使用 openmp 使用 eigen c++ 并行化 Jacobi 算法

Question

我已经根据《数值食谱》一书中描述的例程实现了 Jacobi 算法，但由于我计划处理非常大的矩阵，因此我尝试使用 openmp 对其进行并行化。

void ROTATE(MatrixXd &a, int i, int j, int k, int l, double s, double tau)
{
double g,h;
g=a(i,j);
h=a(k,l);
a(i,j)=g-s*(h+g*tau);
a(k,l)=h+s*(g-h*tau);

}

void jacobi(int n, MatrixXd &a, MatrixXd &v, VectorXd &d )
{
int j,iq,ip,i;
double tresh,theta,tau,t,sm,s,h,g,c;

VectorXd b(n);
VectorXd z(n);

v.setIdentity();    
z.setZero();

#pragma omp parallel for 
for (ip=0;ip<n;ip++)
{   
    d(ip)=a(ip,ip);
    b(ip)=d(ip);
}

for (i=0;i<50;i++) 
{
    sm=0.0;
    for (ip=0;ip<n-1;ip++) 
    {
        #pragma omp parallel for reduction (+:sm)
        for (iq=ip+1;iq<n;iq++)
            sm += fabs(a(ip,iq));
    }
    if (sm == 0.0) {
        break;
    }

    if (i < 3)
    tresh=0.2*sm/(n*n); 
    else
    tresh=0.0;  

    #pragma omp parallel for private (ip,g,h,t,theta,c,s,tau)
    for (ip=0;ip<n-1;ip++)
    {
    //#pragma omp parallel for private (g,h,t,theta,c,s,tau)
        for (iq=ip+1;iq<n;iq++)
        {
            g=100.0*fabs(a(ip,iq));
            if (i > 3 && (fabs(d(ip))+g) == fabs(d[ip]) && (fabs(d[iq])+g) == fabs(d[iq]))
            a(ip,iq)=0.0;
            else if (fabs(a(ip,iq)) > tresh)
            {
                h=d(iq)-d(ip);
                if ((fabs(h)+g) == fabs(h))
                {
                    t=(a(ip,iq))/h;
                }   
                else 
                {
                    theta=0.5*h/(a(ip,iq));
                    t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
                    if (theta < 0.0)
                    {
                        t = -t;
                    }
                    c=1.0/sqrt(1+t*t);
                    s=t*c;
                    tau=s/(1.0+c);
                    h=t*a(ip,iq);

                   #pragma omp critical
                    {
                    z(ip)=z(ip)-h;
                    z(iq)=z(iq)+h;
                    d(ip)=d(ip)-h;
                    d(iq)=d(iq)+h;
                    a(ip,iq)=0.0;


                    for (j=0;j<ip;j++)
                        ROTATE(a,j,ip,j,iq,s,tau);
                    for (j=ip+1;j<iq;j++)
                        ROTATE(a,ip,j,j,iq,s,tau);
                    for (j=iq+1;j<n;j++)
                        ROTATE(a,ip,j,iq,j,s,tau);
                    for (j=0;j<n;j++)
                        ROTATE(v,j,ip,j,iq,s,tau);
                    }

                }
            } 
        }
    }


}

}

我想并行化执行大部分计算的循环以及代码中插入的两个 cmets：

 //#pragma omp parallel for private (ip,g,h,t,theta,c,s,tau)
 //#pragma omp parallel for private (g,h,t,theta,c,s,tau)

是我的尝试。不幸的是，它们最终都产生了不正确的结果。我怀疑问题可能出在这个块：

z(ip)=z(ip)-h;
z(iq)=z(iq)+h;
d(ip)=d(ip)-h;
d(iq)=d(iq)+h;

因为通常这种累积需要减少，但由于每个线程访问数组的不同部分，我不确定这一点。

我不确定我是否以正确的方式进行并行化，因为我最近才开始使用 openmp，因此也欢迎任何建议或建议。

旁注：我知道有更快的特征值和特征向量确定算法，包括 Eigen 中的 SelfAdjointEigenSolver，但这些算法并没有给我在特征向量中所需的精度，而这个算法是。

提前致谢。

编辑：我认为正确的答案是 The Quantum Physicist 提供的答案，因为我所做的并没有减少最大为 4096x4096 的系统的计算时间。无论如何，我更正了代码以使其正常工作，也许对于足够大的系统，它可能会有一些用处。我建议使用计时器来测试

#pragma omp for

实际上减少了计算时间。

Answer 1

我会尽力提供帮助，但我不确定这是否是您问题的答案。

您的代码存在大量问题。我对你的友好建议是：如果你不了解你所做的事情的含义，就不要做平行的事情。

出于某种原因，您似乎认为将所有内容并行处理#pragma for 会使其更快。这是非常错误的。因为产生线程是一件昂贵的事情，并且（相对）花费大量的内存和时间。因此，如果您为每个循环重做 #pragma for，您将为每个循环重新生成线程，这将显着降低程序的速度...除非：您的矩阵非常庞大并且计算时间 >> 比成本产生它们。

当我想以元素方式乘以巨大的矩阵时，我遇到了类似的问题（然后我需要求和以获得量子力学中的一些期望值）。为了使用 OpenMP，我必须将矩阵展平为线性数组，然后将每个数组块分配给一个线程，然后运行一个 for 循环，其中每个循环迭代肯定使用独立于其他元素的元素，我使它们都独立发展。这是相当快的。为什么？因为我从来不用两次重生线程。

为什么你得到错误的结果？我相信原因是因为您不尊重共享内存规则。您有 一些 变量同时被多个线程修改。它藏在某个地方，你必须找到它！比如z这个函数有什么作用？它是否通过引用来获取东西？我在这里看到的：

z(ip)=z(ip)-h;
z(iq)=z(iq)+h;
d(ip)=d(ip)-h;
d(iq)=d(iq)+h;

看起来非常多线程不安全，我不明白你在做什么。您是否返回必须修改的参考？这是线程不安全的秘诀。你为什么不创建干净的数组并处理它们而不是这个？

如何调试：从一个小例子开始（可能是 2x2 矩阵），只使用 2 个线程，并尝试了解发生了什么。使用调试器并定义断点，并检查线程之间共享哪些信息。

还可以考虑使用互斥锁来检查哪些数据在共享时会被破坏。 Here 是怎么做的。

我的建议：除非您打算只生成一次线程，否则不要使用 OpenMP。实际上，我相信 OpenMP 很快就会因为 C++11 而消亡。当 C++ 没有任何本机多线程实现时，OpenMP 很漂亮。所以学习如何使用std::thread，并使用它，如果你需要在线程中运行很多东西，那么学习如何使用std::thread创建一个线程池。 This 是一本学习多线程的好书。