代码之家 › 专栏 › 技术社区 › Nicolas Venkovic

我的CSR的稀疏矩阵多向量(SpMM)函数有什么问题?

blas sparse-matrix linear-algebra

Nicolas Venkovic · 技术社区 · 2 年前

假设CSR存储格式,我有以下C中稀疏矩阵向量(SpMV)乘积的代码:

void dcsrmv(SparseMatrixCSR *A, double *x, double *y) {
  for (int i=0; i<A->m; i++) {
    double sum = 0.;
    for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
      sum += A->val[jj] * x[A->col_idx[jj]];
    }
    y[i] = sum;
  }
}

当我使用数组运行此函数100次时 ecology1 从SuiteParse矩阵集合中,平均运行时间为5.9毫秒,而如果我运行MKL mkl_sparse_d_mv 函数的平均运行时间为8.4ms。

现在,我想为相同的存储格式提供一个稀疏矩阵多矢量(SpMM)乘积函数,我用以下3种方法对其进行了编码:

void dcsrmm0(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  for (int vec=0; vec<nvec; vec++) {
    double *x = X + vec * A->n;
    double *y = Y + vec * A->m;
    for (int i=0; i<A->m; i++) {
      double sum = 0.;
      for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
        sum += A->val[jj] * x[A->col_idx[jj]];
      }
      y[i] = sum;
    }
  }
}

void dcsrmm1(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  for (int i=0; i<A->m; i++) {
    for (int vec=0; vec<nvec; vec++) {
      double sum = 0.;
      double *x = X + vec * A->n;
      double *y = Y + vec * A->m;
      for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
        sum += A->val[jj] * x[A->col_idx[jj]];
      }
      y[i] = sum;
    }
  }
}

和

void dcsrmm2(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  double *x[nvec], *y[nvec];
  for (int vec=0; vec<nvec; vec++) x[vec] = X + vec * A->n;
  for (int vec=0; vec<nvec; vec++) y[vec] = Y + vec * A->m;
  double sum[nvec];
  for (int i=0; i<A->m; i++) {
    for (int vec=0; vec<nvec; vec++) sum[vec] = 0.;
    for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
      for (int vec=0; vec<nvec; vec++) {
        sum[vec] += A->val[jj] * x[vec][A->col_idx[jj]];
      }
    }
    for (int vec=0; vec<nvec; vec++) y[vec][i] = sum[vec];
  }
}

哪里 X 和 Y 按列存储。

然后我用 nvec=10 .在这种情况下 dcsrmm0 平均耗时49.8ms, dcsrmm1 平均耗时38.3ms,并且 dcsrmm2 平均耗时34.8毫秒,而 mkl_sparse_d_mm 只需26.7ms。

所以我的问题是,为什么我的实施 dcsrmv 做得比 mkl_sparse_d_mv ,但是我的 dcsrmm0-2 功能落后太多 mkl_sparse_d_mm ? 我可以做些什么来改进这些功能?

0 回复 | 直到 2 年前

Nicolas Venkovic 2 年前

我发现以下备选方案给了我最好的结果:

void dcsrmm2_improved(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  double *x[nvec];
  for (int vec=0; vec<nvec; vec++) x[vec] = X + vec * A->n;
  for (int i=0; i<A->m; i++) {
    double *y = Y + i;
    for (int vec=0; vec<nvec; vec++) y[vec * A->m] = 0.;
    for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
      double val = A->val[jj];
      int j = A->col_idx[jj];
      for (int vec=0; vec<nvec; vec++) {
        y[vec * A->m] += val * x[vec][j];
      }
    }
  }
}

当我使用数组运行此函数100次时 ecology1 从SuiteParse矩阵集合中,我得到了27.6毫秒的平均运行时间,这与使用 mkl_sparse_d_mm 使用CSR格式。

如果有人知道如何进一步改进,请告诉我。