Skip to content

Ifpack2,Tpetra: Optimize BlockCrs Jacobi & (S)GS per Issue #96

James Willenbring requested to merge mhoemmen:master into master

Created by: mhoemmen

@trilinos/ifpack2 @trilinos/tpetra @amklinv

This is a pull request and not just a commit, because I'm not convinced that the unit tests fully exercise correctness of the algorithms in question. The algorithms do have unit tests, but the test problems might not be as hard as the ones in the applications. I would like interested parties to try out the changes first, before I commit them to master.

Per Issue #96 (closed), I optimized Ifpack2's implementation of block Jacobi and block (Symmetric) Gauss-Seidel, in Ifpack2::Relaxation. I did the following:

  1. Ifpack2 now stores and applies the inverse of the matrix's block diagonal explicitly.
  2. Use the new BlockMultiVector::blockWiseMultiply method to optimize the first sweep of block Jacobi
  3. Use the new BlockMultiVector::blockJacobiUpdate method for subsequent sweeps, rather than hard-coding the low-level kernel in Ifpack2::Relaxation

In (1), "explicitly" means that each block's inverse is computed in place (LU factorization (GETF2) + GETRI) and applied using a matrix-vector product (GEMV). Before, Ifpack2::Relaxation would store the LU factorization of each block (result of GETF2), and apply it using a linear solve (GETRS).

Ifpack2::Relaxation uses the same storage format for the block diagonal, for both (block) Jacobi and (block) (Symmetric) Gauss-Seidel. Thus, I had to change both at the same time. I also had to change Tpetra::Experimental::BlockCrsMatrix's localGaussSeidel interface, because it no longer needs the array of pivots from the LU factorizations (GETRS needs the pivots, but GEMV does not). This makes BlockCrsMatrix's localGaussSeidel interface more consistent with CrsMatrix's interface.

Merge request reports