# 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?