用特征c++实现Jacobi算法的并行化
Parallelization of Jacobi algorithm using eigen c++ using openmp
我已经根据《数值食谱》一书中描述的例程实现了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);
}
}
}
}
}
}
}
我想并行化执行大部分计算的循环和插入代码中的注释:
//#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,所以任何建议或建议也将受到欢迎。
旁注:我知道有更快的算法用于特征值和特征向量的确定,包括特征中的SelfAdjointEigenSolver,但这些并没有给我在特征向量和这个算法中需要的精度。
提前表示感谢。
编辑:我认为正确的答案是量子物理学家提供的答案,因为我所做的并没有减少大小为4096x4096的系统的计算时间。在任何情况下,我纠正了代码,以使其工作,也许对于足够大的系统,它可能是一些用途。我建议使用计时器来测试
#pragma omp for
实际上减少了计算时间。
我会尽力帮忙,但我不确定这就是你问题的答案。
你的代码有很多问题。我给你的友好建议是:如果你不明白你所做的事情的含义,就不要做平行的事情。
出于某种原因,看起来你认为把所有的东西放在并行的#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个线程,并尝试理解发生了什么。使用调试器并定义断点,并检查线程之间共享哪些信息。
还可以考虑使用互斥锁来检查哪些数据在被共享时被破坏了。以下是如何做到这一点。
我的建议:不要使用OpenMP,除非您计划只生成一次线程。实际上,我相信OpenMP很快就会因为c++ 11而消亡。在c++还没有任何本地多线程实现的时候,OpenMP非常漂亮。因此,学习如何使用std::thread
并使用它,如果您需要在线程中运行许多事情,那么学习如何使用std::thread
创建线程池。这是一本学习多线程的好书。
- 如何使用OpenMP并行化此矩阵时间矢量运算
- 如何使用 MPI 的远程内存访问 (RMA) 功能并行化数据聚合?
- 在C++中使用并行化的预期速度是多少(不是 OpenMp,而是 <thread>)
- 如何使用 OpenMP 并行化最近邻搜索
- Malloc 在使用线程并行化 SSH 调用时存在问题
- 如何使用 OpenMP 正确并行化 for 循环?
- 如何将矩阵的行随机复制到内存中的另一个矩阵的过程并行化?
- 如何使用 Pthreads 并行化图像翻转?
- MPI:反复并行化缓冲区
- 是否可以使用OpenMP并行化一个列表,该列表可以在每次迭代中添加新元素
- 如何在Visual Studio中并行化armadillo
- 嵌套循环 OpenMP 并行化、私有索引还是公共索引?
- 如何并行化增加循环的大小
- 迭代卡拉苏巴算法在C++中使用OpenACC并行化和矢量化
- 使用 c++17 算法并行化简单循环
- 使用两个语句并行化 while 循环(弗洛伊德循环检测算法)
- 并行化埃拉托色尼筛 用于查找素数的算法
- 用特征c++实现Jacobi算法的并行化
- 如何有效地并行化分治算法
- c++中带线程的蛮力搜索算法并行化