# Efficient compressed row storage Gauss Seidel C/C++

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:

1. Is this how you would implement a sparse CRS Gauss Seidel solver?
2. 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.

How many multigrid levels are you using? As you go to coarser matrices, they grow progressively more dense, so the cost of doing an individual matrix/vector multiply or Gauss-Seidel smoothing will grow accordingly. Try checking the ratio of the number of non-zero entries to the matrix size at each AMG level. Eventually, you reach a point where any additional gain in the convergence speed of the algorithm is offset by the higher cost of executing the smoother; even though you need fewer outer-level iterations, each one is more expensive. To that end, you can try using fewer multigrid levels.

You may also want to consider using the modified sparse row format. It's a slight variation on CRS which stores all the diagonal entries of the matrix at the beginning of your array v. The advantage of doing so is that you can get rid of the conditional expression

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


inside your loop over k in favor of something more like

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];


Getting rid of conditionals in tight loops can speed things up dramatically because branch prediction is expensive. You can read more about that here. Naturally, you should do some experiments to be sure that it's really helping; these things are famously difficult to predict.