当前位置: 动力学知识库 > 问答 > 编程问答 >

c - Open MP : Symmetric matrix multiplication for Sparse matrices

问题描述:

I am working on parallelizing the Conjugate gradient method to solve sparse matrices. The CG method calls a subroutine "matrixVectorProduct()" in the algorithm. I am trying to parallelise this subroutine specifically. The following is the code I am working with for SYMMETRIC matrices stored in CSR format.

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 after parallelization is:

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;

}

The first omp pragma loop is working as desired. But there seems to be a problem when i parallelize the second loop similarly. It is not giving the correct output.

Can someone please direct me as to what I am doing wrong in the second pragma loop. I am a newbie to multithreading and open MP.

Thanks.

网友答案:

There is a data race in your second loop since multiple threads can update the same vector element. Use:

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

which ensures atomicity of updates.

网友答案:

Try to use two separate arrays instead of one outVec[], and then add their contents together. The reason is that a parallel algorithm may be writing the results of simultaneously performed operations into one and the same memory cell at the same time, thus "spoiling" the content rather then performing summations one after another. You may also have a look at the following webpage (a parallel CG implementation used for development of that code) and contact the developers for an advice:

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

Hope this helps.

分享给朋友:
您可能感兴趣的文章:
随机阅读: