如何解决 OpenMP 内置向量缩减所遇到的这个问题？

Question

我正在为 CSC 存储格式开发自己的稀疏 BLAS 函数实现。为此，我创建了以下数据结构：

typedef struct SparseMatrixCSC {
  int m;            // Number of rows
  int n;            // Number of columns
  int nnz;          // Number of stored values
  double *val;      // Stored values
  int *row_idx;     // Row indices of stored values
  int *col_start;   // Column j contains values with indices in col_start[j]:(col_start[j+1]-1)
} SparseMatrixCSC;

然后，我想使用 OpenMP 来并行化矩阵向量积 (SpMV)。我使用了不同的方法来规避竞争条件。

首先，我使用原子操作如下：

void dcscmv_atomic(SparseMatrixCSC *A, double *x, double *y) {
  for (int i=0; i<A->m; i++) y[i] = 0.;
  #pragma omp parallel for
  for (int j=0; j<A->n; j++) {
    for (int ii=A->col_start[j]; ii<A->col_start[j+1]; ii++) {
      #pragma omp atomic
      y[A->row_idx[ii]] += A->val[ii] * x[j];
    }
  }
}

这工作得很好，但速度非常慢，而且实际上会导致速度减慢而不是加速。

其次，我尝试使用OpenMP的内置向量缩减功能，如下：

void dcscmv_builtin_array_reduction(SparseMatrixCSC *A, double *x, double *y) {
  for (int i=0; i<A->m; i++) y[i] = 0.;
  #pragma omp parallel for reduction(+:y[:A->m])
  for (int j=0; j<A->n; j++) {
    for (int ii=A->col_start[j]; ii<A->col_start[j+1]; ii++) {
      y[A->row_idx[ii]] += A->val[ii] * x[j];
    }
  }
}

虽然此代码可以正确编译，但它可以在单个线程中正常工作，但在使用多线程时会导致分段错误。

第三，由于我无法让 OpenMP 的内置向量缩减工作，我尝试编写自己的缩减代码，如下所示：

void dcscmv_array_reduction_from_scratch(SparseMatrixCSC *A, double *x, double *y) {
  for (int i=0; i<A->m; i++) y[i] = 0.;
  double *YP;
  #pragma omp parallel 
  {
    int P = omp_get_num_threads();
    int p = omp_get_thread_num();
    #pragma omp single
    {
      YP = (double*)mkl_malloc(A->m * P * sizeof(double), sizeof(double));
      for (int i=0; i<A->m*P; i++) YP[i] = 0.;
    }
    #pragma omp for
    for (int j=0; j<A->n; j++) {
      for (int ii=A->col_start[j]; ii<A->col_start[j+1]; ii++) {
        YP[p * A->m + A->row_idx[ii]] += A->val[ii] * x[j];
      }
    }
    #pragma omp for
    for (int i=0; i<A->m; i++) {
      for (int p=0; p<P; p++) {
        y[i] +=  YP[A->m * p + i];
      }
    }
  }
  mkl_free(YP);
}

这个功能有效，但仍然让我速度变慢，尽管不像

dcscmv_atomic

那么糟糕。

我仍然希望，如果我让

dcscmv_builtin_vector_reduction

工作，我也许能够获得一些加速。因此，我的问题是：我在

dcscmv_builtin_vector_reduction

中实现向量缩减的方式有什么问题，以及如何摆脱这种分割错误？

我尝试将 OpenMP 的内置向量缩减功能与

#pragma omp parallel for reduction(+:y[:A->m])

一起应用，但尽管代码已编译，但它会在多线程执行时导致分段错误。

PS：我将 C 语言标识符放在代码块中，stackoverflow 仍然无法正确突出显示代码:-(。

Answer 1

您的问题在于间接寻址，这使得循环不完全并行。正如您所注意到的，简单的解决方案不起作用或者速度很慢。两条出路。

1. 您可以为每个线程提供自己的输出向量，然后最后将它们合并。
2. 使用基于行的算法，而不是基于列的算法。然后输出的每个分量都是可独立计算的，因此是并行的。是的，我知道您是按列存储矩阵，但即使这样也是可能的。您也许可以找到带有 CRS 矩阵的转置乘积的算法。基本上就是这样。