#pragma omp parallel for schedule crashes my program

Question

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?

Would it be fair to assume the crash occurs on one of the `elements.push_back()` calls? — Hasturkun, Oct 27 '13 at 16:10
it's very fair to assume that. the problem is there but i don't understand open mp enough to know what in that calls creates the problem :P — DontCareBear, Oct 27 '13 at 16:12
AFAIK, `std::vector` might not be threadsafe, it's likely being modified by two (or more) threads in parallel. You'll need to establish a barrier to ensure that it's only modified by one thread at a time. — Hasturkun, Oct 27 '13 at 16:15
Related: http://stackoverflow.com/questions/18669296/c-openmp-parallel-for-loop-alternatives-to-stdvector — Hasturkun, Oct 27 '13 at 16:17
Only `std::vector::operator[]` is thread-safe since it does not change the memory organisation of the internal storage. — Hristo Iliev, Oct 28 '13 at 09:59

kangshiyin · Accepted Answer · 2013-10-27T16:38:44.483

3

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.

edited Oct 27 '13 at 16:38

answered Oct 27 '13 at 16:19

kangshiyin

9,681
1
17
29

its no problem if the data would not maintain its original order. first of all thank you! second: how do you concatenate the buff vectors you created for each thread? you need to know how many threads you have before or something like that? and then create an array of this buff vectors, or you need to just declare an inner variable in the loop as private in the opn mp directive? – DontCareBear Oct 27 '13 at 16:29
see [`omp_get_max_threads()`](http://msdn.microsoft.com/en-us/library/62akecca.aspx) and [`omp_get_thread_num()`](http://msdn.microsoft.com/en-us/library/4xs4beeb.aspx) for more info about how to get max possible number of threads and which thread are you in. – kangshiyin Oct 27 '13 at 16:36
i tried this link it worked! thanks!! http://stackoverflow.com/questions/18669296/c-openmp-parallel-for-loop-alternatives-to-stdvector – DontCareBear Oct 27 '13 at 17:59
Yes, that's a more elegant one. – kangshiyin Oct 27 '13 at 18:12

#pragma omp parallel for schedule crashes my program

1 Answers1

update