I am trying to figure out why my sparse (CRS) Gauss Seidel solver is so slow. I tried to find an implementation of the Gauss Seidel method in sparse format online but could only find implementations using dense matrices. Here is my code: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;}Note: I use pre-processor flags to choose how to calculate the error, but I have profiled and this is not the issue. In the above code snippet the arrays r, c, and v represent the sparse matrix A in compressed row storage, b is the right-hand side of Ax=b, n is the dimension of the matrix A, and tol is the error tolerance we want to solve to. My questions are:Is this how you would implement a sparse CRS Gauss Seidel solver?Is there anything obvious that could be changed to speed up this code?Why do I think my implementation is slow? I am using this solver as a smoother for my AMG implementation and the Gauss Seidel smoothing at the highest level takes the majority of the time (significantly more than computing the Galerkin triple matrix product) even though I may only use 72 Gauss Seidel iterations at the highest level. Any help is greatly appreciated.
Efficient compressed row storage Gauss Seidel C/C++
linear solver;c