高效压缩行存储 Gauss Seidel C/C++

计算科学 线性求解器 C
2021-12-22 00:44:23

我想弄清楚为什么我的稀疏(CRS)高斯赛德尔求解器这么慢。我试图在网上找到稀疏格式的 Gauss Seidel 方法的实现,但只能找到使用密集矩阵的实现。这是我的代码:

void gs(int r[], int c[], double v[], double x[], int n, double tol)
{
  //x is initially b in Ax=b
  double *b = new double[n];
  for(int i=0;i<n;i++){b[i] = x[i];}

  int ii = 0, jj = 0;
  double err = 1.0;
  while(err>tol && ii<MAX_ITER){
    //Gauss-Seidel iteration
    double sigma;
    double ajj;
    for(int j=0;j<n;j++){
      sigma = 0.0;
      ajj = 0.0;   //diagonal entry a_jj
      for(int k=r[j];k<r[j+1];k++){
        if(c[k]!=j){
          sigma = sigma + v[k]*x[c[k]];
        }
        else{
          ajj = v[k];
        }
      }
      x[j] = (b[j] - sigma)/ajj;
    }

    if(jj==4){
      //err = error(ar,ac,av,x,b,n);
      err = fast_error(r,c,v,x,b,n,tol);
      jj = 0;
      #if(DEBUG)
        std::cout<<"error: "<<err<<std::endl;
      #endif
    }
    ii++;
    jj++;
  }

  delete[] b;
}

注意:我使用预处理器标志来选择如何计算错误,但我已经分析过,这不是问题。

在上面的代码片段中,数组rcv表示压缩行存储中的稀疏矩阵A,b 是 的右侧Ax=b,n 是矩阵的维数Atol是我们要解决的容错度。我的问题是:

  1. 这就是你将如何实现稀疏 CRS Gauss Seidel 求解器的方法吗?
  2. 有什么明显的东西可以改变来加快这段代码的速度吗?

为什么我认为我的实施很慢?我使用这个求解器作为我的 AMG 实现的平滑器,最高级别的 Gauss Seidel 平滑花费大部分时间(比计算 Galerkin 三重矩阵乘积要多得多),即使我可能只使用 72 次 Gauss Seidel 迭代最高级别。任何帮助是极大的赞赏。

1个回答

您使用了多少个多重网格级别?当您使用更粗糙的矩阵时,它们会逐渐变得更密集,因此进行单个矩阵/向量乘法或高斯-赛德尔平滑的成本将相应增加。尝试检查每个 AMG 级别的非零条目数与矩阵大小的比率。最终,您会达到一个点,即算法收敛速度的任何额外增益都会被执行更平滑的更高成本所抵消;即使您需要较少的外层迭代,但每次迭代都更昂贵。为此,您可以尝试使用更少的多重网格级别。

您可能还想考虑使用修改后的稀疏行格式这是 CRS 的一个细微变化,它将矩阵的所有对角线条目存储在 array 的开头v这样做的好处是可以去掉条件表达式

if(c[k]!=j){
    sigma += v[k] * x[c[k]];
} else {
    ajj = v[k];
}

在你的循环k中支持更喜欢的东西

sigma = 0.0;
for (int k = r[j]; k < r[j + 1]; k++) {
    sigma += v[k] * x[c[k]];
}
x[j] = (b[j] - sigma)/v[j];

摆脱紧循环中的条件可以显着加快速度,因为分支预测很昂贵。您可以在此处阅读更多相关信息。当然,您应该做一些实验以确保它确实有帮助;众所周知,这些事情很难预测。