CODEWITHSUNDEEP
c++openmp
DontCareBear
DontCareBear

Reputation: 834

#pragma omp parallel for schedule crashes my program

I am building a plugin for autodesk maya 2013 in c++. I have to solve a set of optimization problems as fast as i can. I am using open MP for this task. the problem is I don't have very much experience with parallel computing. I tried to use:

#pragma omp parallel for schedule (static)

on my for loops (without enough understanding of how it's supposed to work) and it worked very well for some of my code, but crashed another portion of my code.

Here is an example of a function that crashes because of the omp directive:

void PlanarizationConstraint::fillSparseMatrix(const Optimizer& opt, vector<T>& elements, double mu)
{
    int size = 3;
    #pragma omp parallel for schedule (static)
    for(int i = 0; i < opt.FVIc.outerSize(); i++)
    {
        int index = 3*i;
        Eigen::Matrix<double,3,3> Qxyz = Eigen::Matrix<double,3,3>::Zero();
        for(SpMat::InnerIterator it(opt.FVIc,i); it; ++it)
        {
            int face = it.row();
            for(int n = 0; n < size; n++)
            {
                Qxyz.row(n) += N(face,n)*N.row(face);
                elements.push_back(T(index+n,offset+face,(1 - mu)*N(face,n)));
            }
        }

        for(int n = 0; n < size; n++)
        {
            for(int k = 0; k < size; k++)
            {
                elements.push_back(T(index+n,index+k,(1-mu)*Qxyz(n,k)));
            }
        }
    }

    #pragma omp parallel for schedule (static)
    for(int j = 0; j < opt.VFIc.outerSize(); j++)
    {
        elements.push_back(T(offset+j,offset+j,opt.fvi[j]));
        for(SpMat::InnerIterator it(opt.VFIc,j); it; ++it)
        {
            int index = 3*it.row();
            for(int n = 0; n < size; n++)
            {
                elements.push_back(T(offset+j,index+n,N(j,n)));
            }
        }
    }
}

And here is an example of code that works very well with those directives (and is faster because of it)

Eigen::MatrixXd Optimizer::OptimizeLLGeneral()
{
    ConstraintsManager manager;
    SurfaceConstraint surface(1,true);
    PlanarizationConstraint planarization(1,true,3^Nv,Nf);
    manager.addConstraint(&surface);
    manager.addConstraint(&planarization);
    double mu = mu0;
    for(int k = 0; k < iterations; k++)
    {
        #pragma omp parallel for schedule (static)
        for(int j = 0; j < VFIc.outerSize(); j++)
        {
            manager.calcVariableMatrix(*this,j);
        }
        #pragma omp parallel for schedule (static)
        for(int i = 0; i < FVIc.outerSize(); i++)
        {
            Eigen::MatrixXd A = Eigen::Matrix<double, 3, 3>::Zero();
            Eigen::MatrixXd b = Eigen::Matrix<double, 1, 3>::Zero();
            manager.addLocalMatrixComponent(*this,i,A,b,mu);
            Eigen::VectorXd temp = b.transpose();
            Q.row(i) = A.colPivHouseholderQr().solve(temp);
        }
        mu = r*mu;
    }
    return Q;
}

My question is what makes one function work so well with the omp directive and what makes the other function crash? what is the difference that makes the omp directive act differently?

Upvotes: 2

Views: 3415

Answers (1)

kangshiyin
kangshiyin

Reputation: 9781

Before using openmp, you pushed back some data to the vector elements one by one. However, with openmp, there will be several threads running the code in the for loop in parallel. When more than one thread are pushing back data to the vector elements at the same time, and when there's no code to ensure that one thread will not start pushing before another one finishes, problem will happen. That's why your code crashes.

To solve this problem, you could use local buff vectors. Each thread first push data to its private local buffer vector, then you can concatenate these buffer vectors together into a single vector.

You will notice that this method can not maintain the original order of the data elements in the vector elements. If you want to do that, you could calculate each expected index of the data element and assign the data to the right position directly.

update

OpenMP provides APIs to let you know how many threads you use and which thread you are using. See omp_get_max_threads() and omp_get_thread_num() for more info.

Upvotes: 3

Related Questions