Open MP：稀疏矩陣的對稱矩陣乘法

我正在並行化共軛梯度法來求解稀疏矩陣。 CG方法在算法中調用子程序「matrixVectorProduct（）」。我正在嘗試專門對此子例程進行並行化。以下是我正在使用的以CSR格式存儲的SYMMETRIC矩陣的代碼。Open MP：稀疏矩陣的對稱矩陣乘法

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; 
}

第一OMP編譯指示循環正在按期望的。但似乎有一個問題，當我平行第二個循環類似。它沒有給出正確的輸出。

有人可以指示我在做什麼我錯了在第二個雜注循環。我是多線程和開放MP的新手。

謝謝。

來源

2016-04-30 a_myth

嘗試使用兩個單獨的數組而不是一個outVec []，然後將它們的內容添加在一起。原因在於並行算法可能會將同時執行的操作的結果同時寫入同一個存儲器單元，從而「破壞」內容而不是一個接一個地執行求和。您也可以看看以下網頁（用於該代碼的開發並行CG實現），並聯系開發者的建議：

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

希望這有助於。

來源

2016-05-01 22:02:36 SparseSolverCodes

由於多個線程可以更新相同的向量元素，因此第二個循環中存在數據競爭。使用：

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

它確保更新的原子性。

來源

2016-05-02 12:26:38

Open MP：稀疏矩陣的對稱矩陣乘法

回答

相關問題