• 如果您觉得本站非常有看点,那么赶紧使用Ctrl+D 收藏吧

并行使用openMP的SVD分解不如预期的那样执行

openmp 来源:Teol11 5次浏览

我最近编写了一个基于“单侧雅可比旋转”算法的并行SVD分解例程。代码正常工作,但速度非常慢。事实上,它应该利用内部循环的并行性for(int g=0;g<n;g++),但在注释掉#pragma omp paralell for指令时,我可以看到性能略有下降。换句话说,并行处理没有明显的速度(代码与4个线程并行运行)。注1:几乎所有的工作都集中在涉及矩阵A和V的相对较大的三个后续循环中。并行使用openMP的SVD分解不如预期的那样执行

for(h=0;h<N;h++) 
{ 
    p+=A[h+N*i]*A[h+N*j];//columns dot product:Ai * Aj 
    qi+=A[h+N*i]*A[h+N*i];// ||Ai||^2 
    qj+=A[h+N*j]*A[h+N*j];// ||Aj||^2 
} 

double Ahi,Vhi; 
for(h=0;h<N;h++)//...rotate Ai & Aj (only columns i & j are changend) 
{ 
    Ahi=A[h+N*i]; 
    A[h+N*i]=cs*A[h+N*i]+sn*A[h+N*j]; 
    A[h+N*j]=-sn*Ahi+cs*A[h+N*j]; 
} 
//store & update rotation matrix V (only columns i & j are updated) 
for(h=0;h<N;h++) 
{ 
    Vhi=V[h+N*i]; 
    V[h+N*i]=cs*V[h+N*i]+sn*V[h+N*j]; 
    V[h+N*j]=-sn*Vhi+cs*V[h+N*j]; 
} 

所有并行应该有剥削,但并非如此。我不明白为什么。

注2:在Windows(cygWin编译器)和Linux(GCC)平台上都会发生同样的情况。

注3:矩阵按列重大阵列

所以我希望在查明原因的平行不被剥削一些帮助表示。我错过了什么?平行中有一些隐藏的开销,因为我看不到?

非常感谢您的任何建议

int sweep(double* A,double*V,int N,double tol) 
{ 
static int*I=new int[(int)ceil(0.5*(N-1))]; 
static int*J=new int[(int)ceil(0.5*(N-1))]; 
int ntol=0; 
for(int r=0;r<N;r++) //fill in i,j indexes of parallel rotations in vectors  I & J 
{ 
    int k=r+1; 
    if (k==N) 
    { 
    for(int i=2;i<=(int)ceil(0.5*N);i++){ 
     I[i-2]=i-1; 
     J[i-2]=N+2-i-1; 
     } 
    } 
    else 
    { 
     for(int i=1;i<=(int)ceil(0.5*(N-k));i++)I[i-1]=i-1; 
     for(int i=1;i<=(int)ceil(0.5*(N-k));i++)J[i-1]=N-k+2-i-1; 
     if(k>2) 
     { 
      int j=(int)ceil(0.5*(N-k)); 
      for(int i=N-k+2;i<=N-(int)floor(0.5*k);i++){ 
       I[j]=i-1; 
       J[j]=2*N-k+2-i-1; 
       j++; 
      } 
     } 
    } 

    int n=(k%2==0)?(int)floor(0.5*(N-1)):(int)floor(0.5*N); 

    #pragma omp parallel for schedule(dynamic,5) reduction(+:ntol) default(none) shared(std::cout,I,J,A,V,N,n,tol) 
    for(int g=0;g<n;g++) 
    { 
     int i=I[g]; 
     int j=J[g]; 
     double p=0; 
     double qi=0; 
     double qj=0; 
     double cs,sn,q,c; 
     int h; 
     for(h=0;h<N;h++) 
     { 
      p+=A[h+N*i]*A[h+N*j];//columns dot product:Ai * Aj 
      qi+=A[h+N*i]*A[h+N*i];// ||Ai||^2 
      qj+=A[h+N*j]*A[h+N*j];// ||Aj||^2 
     } 
     q=qi-qj; 
     if(p*p/(qi*qj)<tol) ntol++; //if Ai & Aj are orthogonal enough... 
     else      //if Ai & Aj are not orthogonal enough then... rotate them 
     { 
      c=sqrt(4*p*p+q*q); 
      if(q>=0){ 
       cs=sqrt((c+q)/(2*c)); 
       sn=p/(c*cs); 
      } 
      else{ 
       sn=(p>=0)?sqrt((c-q)/2/c):-sqrt((c-q)/2/c); 
       cs=p/(c*sn); 
      } 
      //...rotate Ai & Aj (only columns i & j are changend) 
      double Ahi,Vhi; 
      for(h=0;h<N;h++) 
      { 
       Ahi=A[h+N*i]; 
       A[h+N*i]=cs*A[h+N*i]+sn*A[h+N*j]; 
       A[h+N*j]=-sn*Ahi+cs*A[h+N*j]; 
      } 
      //store & update rotation matrix V (only columns i & j are updated) 
      for(h=0;h<N;h++) 
      { 
       Vhi=V[h+N*i]; 
       V[h+N*i]=cs*V[h+N*i]+sn*V[h+N*j]; 
       V[h+N*j]=-sn*Vhi+cs*V[h+N*j]; 
      } 
     } 
    } 
} 
if(2*ntol==(N*(N-1)))return(1);//if each columns of A is orthogonal enough  to each other stop sweep 

return(0); 
} 


===========解决方案如下:

的FLOPS的奇异值分解(SVD)应该为O(N^3)而读为O(N^2)。这意味着可以对算法进行并行化并使其与核心数量保持一致。

但是,您的SVD实现是内存带宽限制,这意味着它无法与单个套接字系统的内核数量成正比。目前,您的三个嵌套循环遍历整个范围N,这会导致大量数据在下次需要重用时从缓存中逐出。

为了计算绑定你将不得不改变你的算法,这样,而不是在长条工作/尺寸N的点产品,它采用循环平铺最大化n^3(其中n是大小的块)计算在大小为n^2的缓存中。

这里是paper for parallel SVD computation来自谷歌搜索其说,在我们提出的块JRS算法的抽象

的关键点是通过在矩阵的块执行计算(B行)重新使用加载的数据到高速缓冲存储器中,而不是如JRS迭代算法中的矢量条带。

I used loop tiling/blocks to achieve good scaling with the number of cores for cholesky decomposition。

请注意,使用tiles/blocks也会提高单线程代码的性能。


版权声明:本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系管理员进行删除。
喜欢 (0)