2
$\begingroup$

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?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for the information. I will take a closer look. $\endgroup$
    – vydesaster
    Sep 4, 2019 at 17:05
  • 2
    $\begingroup$ 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! $\endgroup$
    – Wolfgang Bangerth
    Sep 4, 2019 at 20:29
  • 1
    $\begingroup$ 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. $\endgroup$
    – EMP
    Sep 4, 2019 at 20:44
  • $\begingroup$ 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. $\endgroup$
    – Kyle Mandli
    Sep 4, 2019 at 21:19
  • $\begingroup$ 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. $\endgroup$
    – dweber
    Sep 5, 2019 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.