Open MP：稀疏矩阵的对称矩阵乘法答案

【问题标题】：Open MP : Symmetric matrix multiplication for Sparse matricesOpen MP：稀疏矩阵的对称矩阵乘法
【发布时间】：2016-08-25 17:31:21
【问题描述】：

我正在研究并行化共轭梯度法来求解稀疏矩阵。 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;
}

第一个 omp pragma 循环按预期工作。但是当我类似地并行化第二个循环时似乎存在问题。它没有给出正确的输出。

有人可以指导我在第二个编译指示循环中做错了什么。我是多线程和开放 MP 的新手。

谢谢。

【问题讨论】：

标签： c multithreading parallel-processing openmp sparse-matrix

【解决方案1】：

在您的第二个循环中存在数据竞争，因为多个线程可以更新同一个向量元素。使用：

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

确保更新的原子性。

【讨论】：