Open MP: Symmetrische Matrixmultiplikation für Sparse-Matrizen

Question

Ich arbeite an der Parallelisierung der Conjugate-Gradientenmethode, um dünn besetzte Matrizen zu lösen. Die CG-Methode ruft im Algorithmus eine Subroutine "matrixVectorProduct()" auf. Ich versuche diese Unterroutine gezielt parallel zu machen. Das Folgende ist der Code, mit dem ich für SYMMETRIC-Matrizen arbeite, die im CSR-Format gespeichert sind.

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

My Code nach Parallelisierung ist:

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

Die erste omp pragma Schleife wie gewünscht funktioniert. Aber es scheint ein Problem zu geben, wenn ich die zweite Schleife ähnlich parallelisiere. Es gibt nicht die richtige Ausgabe.

Kann mir bitte jemand sagen, was ich in der zweiten Pragma-Schleife falsch mache. Ich bin ein Neuling Multithreading und MP öffnen.

Danke.

Answer 1

Versuchen Sie, zwei separate Arrays anstelle von einem outVec [] zu verwenden, und fügen Sie dann ihre Inhalte zusammen. Der Grund dafür ist, dass ein paralleler Algorithmus die Ergebnisse von gleichzeitig ausgeführten Operationen in ein und dieselbe Speicherzelle gleichzeitig schreiben kann, wodurch der Inhalt "verdorben" wird und dann nacheinander Summierungen durchgeführt werden. Sie können auch einen Blick auf die folgende Website haben (eine parallele Implementierung CG für die Entwicklung dieses Codes) und kontaktieren Sie die Entwickler für einen Rat:

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

Hoffnung, das hilft.

Answer 2

In Ihrer zweiten Schleife gibt es ein Datenrennen, da mehrere Threads das gleiche Vektorelement aktualisieren können. Verwenden Sie:

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

, die Atomizität von Updates gewährleistet.