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