时间戳

首页 > 解决方案 > openmp - 并行向量矩阵乘积

c - openmp - 并行向量矩阵乘积

问题描述

我正在使用外积方法计算向量矩阵（矩阵是稀疏的并存储在 CSR 中）乘积（不完全是乘积，计算最短距离的轻微变化）。我是并行编程的新手，基本上试图了解使用并行 for 部分与使用任务的更新 VS 的关键部分和进行缩减之间的区别。哪种方法更好，为什么？

注意：此函数调用包含一个 omp parallel 和一个 omp single。

使用并行方法，我这样做：

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq) {
    initialize_T(tReq);
    int index;
    #pragma omp parallel for schedule(static, BLOCK_SIZE)
    for(int i=0;i<N;i++) {
        int num_edges = A->row_ptr[i+1] - A->row_ptr[i];
        index = 0;
        if(num_edges) {
            if(T[i] != INFINITY && tB[i] != INFINITY) {
                for(int j=0;j<num_edges;j++) {
                    index = A->col_ind[A->row_ptr[i] + j];
                    #pragma omp critical 
                    tReq[index] = min(tReq[index], T[i]+A->val[A->row_ptr[i]+j]);      
                }
            }
        }
    }
    return tReq;
}

使用基于任务的减少方法，这本质上是我的想法：

int size = N/BLOCK_SIZE + 1;
double C[size][N];

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq, int size, double C[][N]) {

    int index;

    for(int i=0;i<size;i++) {
        for(int j=0;j<N;j++) {
            C[i][j] = INFINITY;
            tReq[j] = INFINITY;
        }
    }

    int k;

    for(k=0;k<size-1; k++) {
        #pragma omp task firstprivate(k) depend(inout: C[k]) 
        {
            int index;
            int delim;   
            delim = (k==size-1) ? N : BLOCK_SIZE;
            printf("kk is %d\n", k*BLOCK_SIZE);
            // printf("k is %d Delim is %d for thread %d\n", k, delim, omp_get_thread_num());
            for(int i=0;i<delim; i++) {
                int num_edges = A->row_ptr[k*BLOCK_SIZE + i+1] - A->row_ptr[k*BLOCK_SIZE + i];
                index = 0;
                if(num_edges) {
                    if(T[k*BLOCK_SIZE + i] != INFINITY && tB[k*BLOCK_SIZE + i] != INFINITY) {           
                        for(int j=0;j<num_edges;j++) {
                            index = A->col_ind[A->row_ptr[k*BLOCK_SIZE + i] + j];
                            {
                            C[k][index] = min(C[k][index], T[k*BLOCK_SIZE + i]+A->val[A->row_ptr[k*BLOCK_SIZE + i]+j]);                 
                            }
                        }
                    }       
                }   
            }
        }
    }    

    #pragma omp taskwait

    for(int i=0; i<N; i++) {
        {
            double minimum = INFINITY;
            for(int j=0;j<size;j++) {
                if(C[j][i] < minimum) {
                    minimum = C[j][i];
                }
            }
            tReq[i] = minimum;
        }
    }

    return tReq;
}

本质上，与基于任务的方法相比，使用并行有什么缺点吗？

标签： copenmpsparse-matrixmatrix-multiplication

解决方案

没错，您基本上有两个选择：保护数据更新或使用线程特定的副本。但是，您可以为每个选项做得更好：

使用受保护的更新时，您应该只保护绝对必要的内容。大多数情况下，您可以使用初始原子检查来防止关键区域，类似于双重检查锁定模式。

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq) {
    initialize_T(tReq);
    #pragma omp parallel for schedule(static, BLOCK_SIZE)
    for(int i=0;i<N;i++) {
        int num_edges = A->row_ptr[i+1] - A->row_ptr[i];
        if (num_edges) {
            if(T[i] != INFINITY && tB[i] != INFINITY) {
                for(int j=0;j<num_edges;j++) {
                    // !WARNING! You MUST declare index within the parallel region
                    // or explicitly declare it private to avoid data races!
                    int index = A->col_ind[A->row_ptr[i] + j];
                    double tmp = T[i] + A->val[A->row_ptr[i]+j];
                    double old;
                    #pragma omp atomic
                    old = tReq[index];
                    if (tmp < old) {
                        #pragma omp critical
                        {
                            tmp = min(tReq[index], tmp);
                            // Another atomic ensures that the earlier read
                            // outside of critical works correctly
                            #pragma omp atomic
                            tReq[index] = tmp;
                        }
                    }
                }
            }
        }
    }
    return tReq;
}

注意：不幸的是，没有 OpenMP/C 不支持直接原子最小值。

另一种方法是减少，甚至标准本身也支持。所以没有必要重新发明工作共享等。您可以简单地执行以下操作：

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq) {
    initialize_T(tReq);
    #pragma omp parallel for schedule(static, BLOCK_SIZE) reduction(min:tReq)
    for(int i=0;i<N;i++) {
        int num_edges = A->row_ptr[i+1] - A->row_ptr[i];
        if (num_edges) {
            if(T[i] != INFINITY && tB[i] != INFINITY) {
                for(int j=0;j<num_edges;j++) {
                    // !WARNING! You MUST declare index within the parallel region
                    // or explicitly declare it private to avoid data races!
                    int index = A->col_ind[A->row_ptr[i] + j];
                    tReq[index] = min(tReq[index], T[i]+A->val[A->row_ptr[i]+j]);
                }
            }
        }
    }
    return tReq;
}

OpenMP 将神奇地创建线程本地副本tReq并在最后合并（减少）它们。

哪个版本更适合您取决于目标数组的大小和写入速率。critical如果你经常写，减少将是有益的，因为它不会因//atomic糟糕的缓存而减慢。如果您有一个巨大的目标数组或没有那么多更新迭代，第一个解决方案会变得更有趣，因为创建和减少数组的相对开销。

推荐阅读