# Implementation of Jacobi iteration

I have implemented the Jacobi iteration in C++ using a dense vector and a sparse matrix in CSR format. The code is as follows:

    {
T omega = 2.0 / 3.0;
std::vector<T> temp = result;

for (int niter = 0; niter < maxit; ++niter)
{
for (int i = 0; i < matrix.get_rowptr()->size() - 1; ++i)
{
T rsum = 0.0;
T diag = 0.0;

for (int j = (*matrix.get_rowptr())[i]; j < (*matrix.get_rowptr())[i + 1]; ++j)
{
if ((*matrix.get_columnindex())[j] == i)
diag = (*matrix.get_value())[j];
else
rsum += (*matrix.get_value())[j] * temp[(*matrix.get_columnindex())[j]];
}

if (diag != 0.0)
result[i] = temp[i] + omega * ((b[i] - rsum) / diag);
}

temp = result;
}
}


I have profiled my application and this function is the one with the most time amount. It makes sense as this function is called a couple of times.
Now I am looking for a more efficient way to implement it but was not able to find a solution. By the way, the for loop with i is parallized using OpenMP. I have removed this piece of code for the post.
Any idea how to speed this up? Would it make sense to move only this function to CUDA?

If you really want to speed it up, I'd suggest changing your iteration to a SSOR iteration. This can converge much faster than jacobi, and if you are already OpenMP parallelized, there's not much you can do to improve performance aside from maybe domain decomposition in MPI if your problem size is large enough. The fastest and easiest way to speed this up would be by improving the iterations through either Gauss-Seidel or SSOR iterations, both of which are very easy to implement.

• Thanks for the information. I will take a closer look. – vydesaster Sep 4 '19 at 17:05
• Well, as a matter of fact, use a Krylov space method rather than some fixed point iteration. That's a general rule: If you really want to speed something up, change the algorithm rather than the implementation of the algorithm! – Wolfgang Bangerth Sep 4 '19 at 20:29
• Those do take much more coding to implement and I'd suggest actually understanding the theory behind them before implementing, but from a ideal standpoint, you are of course correct. – EMP Sep 4 '19 at 20:44
• I agree with both the answer and @WolfgangBangerth's comments but I was wondering if you have just tried not storing the matrix itself? Even if you use a sparse representation one of the advantages of using something simple like Jacobi (or SOR) is that the matrix never needs to be explicitly formed often and therefore avoids a lot of need for special caching or looks ups. – Kyle Mandli Sep 4 '19 at 21:19
• Paralellization for Gauss-Seidel is more difficult, at least for general matrices because of the dependencies between the DOFs. For SSOR you need to take care about the update order of the symmetric iterations, as they should be reverse to each other. – dweber Sep 5 '19 at 11:31