c++ - #pragma omp parallel for schedule 使我的程序崩溃

Question

我正在用 C++ 为 Autodesk Maya 2013 构建一个插件。我必须尽快解决一组优化问题。我正在使用 open MP 来完成这项任务。问题是我在并行计算方面没有太多经验。我尝试使用:

#pragma omp parallel for schedule (static)

在我的 for 循环中(没有充分理解它应该如何工作)，它对我的​​一些代码运行良好，但使我的另一部分代码崩溃。

这是一个由于 omp 指令而崩溃的函数示例:

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)));
            }
        }
    }
}

下面是一个代码示例，它与这些指令配合得很好(因此速度更快)

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;
}

我的问题是，是什么使一个函数与 omp 指令配合得如此好，又是什么使另一个函数崩溃？使 omp 指令行为不同的区别是什么？

Answer 1

在使用 openmp 之前，您将一些数据一个一个地推回到 vector elements 中。但是，使用 openmp，将有多个线程并行运行 for 循环中的代码。当多个线程同时将数据推回 vector elements 时，如果没有代码确保一个线程不会在另一个线程完成之前开始推送，就会出现问题。这就是您的代码崩溃的原因。

要解决这个问题，您可以使用局部增益 vector 。每个线程首先将数据推送到其私有(private)的本地缓冲区 vector ，然后您可以将这些缓冲区 vector 连接成一个 vector 。

您会注意到，此方法无法保持 vector elements 中数据元素的原始顺序。如果你想这样做，你可以计算数据元素的每个预期索引，并直接将数据分配到正确的位置。

更新

OpenMP 提供 API 让您知道您使用了多少个线程以及您正在使用哪个线程。参见 omp_get_max_threads()和 omp_get_thread_num()了解更多信息。