Open MP：Sparse行列の対称行列乗算

Question

Sparse行列を解決するためにConjugate勾配法を並列化する作業を進めています。 CGメソッドはアルゴリズムの中で "matrixVectorProduct（）"というサブルーチンを呼び出します。私はこのサブルーチンを具体的に並列化しようとしています。以下は、CSRフォーマットで格納されたSYMMETRIC行列のためのコードです。

void matrixVectorProduct(MTX *MAT, double* inVec, double* outVec){ 

int i,j, ckey; 

if((matcode[1] == 'X')&&(matcode[3] == 'S')) 
{ 
    //Initialize outVec to zeros 
    for(int j=0;j<MAT->nrows;j++) 
      outVec[j] = 0.0; 

     for(i=0;i<MAT->nrows;i++) 
      for(ckey=MAT->row_ptr[i];ckey<MAT->row_ptr[i+1];ckey++) { 
       j = MAT->JA[ckey]; 
       outVec[i] = outVec[i] + MAT->val[ckey] * inVec[j]; 
      } 

     for(int i=0;i<MAT->nrows;i++) 
      for(int ckey=MAT->row_ptr[i];ckey<MAT->row_ptr[i+1];ckey++) { 
       j = MAT->JA[ckey]; 
       if(j!=i) 
        outVec[j] += MAT->val[ckey] * inVec[i];; 
      } 
} 
else 
{ 
    fprintf(stderr,"\n Not a symmetric Matrix. CG method not applicable\n"); 
    exit(1); 
}  
return;} 

並列後、マイコードである：

void matrixVectorProduct(MTX *MAT, double* inVec, double* outVec){ 

int i,j, ckey; 

if((matcode[1] == 'X')&&(matcode[3] == 'S')) 
{ 
    //Initialize outVec to zeros 
    for(int j=0;j<MAT->nrows;j++) 
      outVec[j] = 0.0; 

    #pragma omp parallel 
    { 
     #pragma omp for private(ckey,j) schedule(static) 
     for(i=0;i<MAT->nrows;i++) { 
      double zi = 0.0; 
      for(ckey=MAT->row_ptr[i];ckey<MAT->row_ptr[i+1];ckey++) { 
       j = MAT->JA[ckey]; 
       zi = zi + MAT->val[ckey] * inVec[j]; 
      } 
     outVec[i] += zi; 
     } 
    } 

    #pragma omp parallel 
    { 
     #pragma omp for private(ckey,j) schedule(static) 
     for(int i=0;i<MAT->nrows;i++) 
      for(int ckey=MAT->row_ptr[i];ckey<MAT->row_ptr[i+1];ckey++) { 
       j = MAT->JA[ckey]; 
       if(j!=i) 
        outVec[j] += MAT->val[ckey] * inVec[i];; 
      } 
    } 

} 
else 
{ 
    fprintf(stderr,"\n Not a symmetric Matrix. CG method not applicable\n"); 
    exit(1); 
} 

return; 
} 

所望のように、第1 OMPプラグマループが動作しています。しかし、私は同様に2番目のループを並列化するときに問題があるようです。正しい出力が得られていません。

2番目のプラグマループで何が間違っているのか教えてもらえますか？私はマルチスレッドとオープンMPへの初心者です。

ありがとうございました。

Answer 1

1つのoutVec []の代わりに2つの別々の配列を使用し、それらの内容を一緒に追加してみてください。その理由は、並列アルゴリズムが、同時に実行される演算の結果を同一のメモリセルに同時に書き込む可能性があり、その結果、次々に加算を実行するのではなく、内容を「損なわせる」ことである。 また、次のWebページ（そのコードの開発のために使用されるパラレルCGの実装）を見ていると助言のための開発者に連絡することがあります。

http://members.ozemail.com.au/~comecau/CMA_Sparse.htm

・ホープ、このことができます。

Answer 2

複数のスレッドが同じベクトル要素を更新できるため、2番目のループにデータ競合があります。アップデートのアトミック性を保証するために、

if (j != i) { 
    #pragma omp atomic 
    outVec[j] += MAT->val[ckey] * inVec[i]; 
} 

を使用してください。