[petsc-users] Matrix-free Multigrid

[Timothée Nicolas]

Hello all,

Eventually, I would like to use matrix-free methods in my multigrid solver,
because I have a large problem and matrix memory storage is an issue. Since
all my matrices represent operators, it is quite straightforward to make
them matrix-free, and according to the manual, I can use multigrid on them.
However, depending on the method chosen in KSP, different matrix operations
are used, which the user has to code by himself.

That is, in my understanding, coding the subroutine for MATOP_MULT (which
in principle is sufficient to describe the matrix) is actually not
sufficient, because the KSP may want internally to get, e.g., the diagonal
with MatGetDiagonal, in which case I need a routine to define the operation
MATOP_GET_DIAGONAL (as far as I understand). This means I have to know all
the matrix related routines used internally by the KSP (which will depend
on the ksp method as well as the preconditioner choice), and hard code the
corresponding routines. After all, one could probably do this, but this is
really tedious work and I have the feeling it is not the approach intended
by the developers of PETSc.

Would someone have an advice on how to do this efficiently ? Maybe I am
totally wrong ? But when I tried to do a KSP with multigrid preconditioning
on a matrix-free matrix, I got errors complaining about MatGetDiagonal and
this sort of things.

Best

Timothee

[Dave May]

Hi Timothee,

Your thinking is correct -- however it is not as bad as you imagine.

If you want to define a matrix free definition of your operator within
geometric MG,
the Mat is required to support several methods in order to be used in
conjunction with smoother.

For Krylov methods, you definitely need to define MatMult.
Assuming you are happy to use a Jacobi preconditioner, you will need to
implement MatGetDiagonal.

I use MF operators within geometric multigrid (on all levels) and find only
these
two operations are required to support using smoothers such as
{richardson,cheby,cg,gmres}+jacobi.

If wish to define your restriction and prolongation operators in a
matrix-free manner
you will need to additionally define the operation MatMultAdd which is
called from
MatInterpolateAdd() inside PCMG. If R = P^T, you will also need to define
MatMultTranspose and MatMultTransposeAdd


Cheers
Dave

Post by Timothée Nicolas
Hello all,
Eventually, I would like to use matrix-free methods in my multigrid
solver, because I have a large problem and matrix memory storage is an
issue. Since all my matrices represent operators, it is quite
straightforward to make them matrix-free, and according to the manual, I
can use multigrid on them. However, depending on the method chosen in KSP,
different matrix operations are used, which the user has to code by
himself.
That is, in my understanding, coding the subroutine for MATOP_MULT (which
in principle is sufficient to describe the matrix) is actually not
sufficient, because the KSP may want internally to get, e.g., the diagonal
with MatGetDiagonal, in which case I need a routine to define the operation
MATOP_GET_DIAGONAL (as far as I understand). This means I have to know all
the matrix related routines used internally by the KSP (which will depend
on the ksp method as well as the preconditioner choice), and hard code the
corresponding routines. After all, one could probably do this, but this is
really tedious work and I have the feeling it is not the approach intended
by the developers of PETSc.
Would someone have an advice on how to do this efficiently ? Maybe I am
totally wrong ? But when I tried to do a KSP with multigrid preconditioning
on a matrix-free matrix, I got errors complaining about MatGetDiagonal and
this sort of things.
Best
Timothee