2

我正在用 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 指令行为不同的区别是什么？

4

1 回答 1

3

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

要解决此问题，您可以使用局部 buff 向量。每个线程首先将数据推送到其私有的本地缓冲区向量，然后您可以将这些缓冲区向量连接在一起成为一个向量。

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

更新

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

于 2013-10-27T16:19:31.080 回答