[petsc-users] V-cycle multigrid with matrix shells

[Sylvain Barbot]

Hi Barry,

I have tested the direct solver with the shell matrix and it works
very well. In my code, the matrix-multiply operation is implemented
with the function "matrixantiplanefv". I am now trying to do it with
multi-grid. I am keeping a direct solver solely for the coarse-level.
I will do this coarse solve matrix-free with the function
"matrixantiplanecoarsefv". I am also setting up all the other
operations (residual, smoother, restriction, interpolation) with
matrix-free methods. The setup is done in the subroutine
"green_init_mg".

First I set up a PCMG preconditionner:

CALL KSPCreate(PETSC_COMM_WORLD,c%ksp,ierr);
CALL KSPGetPC(c%ksp,c%pc,ierr);
CALL PCSetType(c%pc,PCMG,ierr);
CALL PCMGSetLevels(c%pc,c%nlevels,PETSC_NULL_OBJECT,ierr);
CALL PCMGSetType(c%pc,PC_MG_MULTIPLICATIVE,ierr);
CALL PCMGSetCycleType(c%pc,PC_MG_CYCLE_V,ierr);
CALL PCMGSetNumberSmoothUp(c%pc,1,ierr);
CALL PCMGSetNumberSmoothDown(c%pc,1,ierr);


Then I set up the coarse solver:

CALL MatCreateShell(PETSC_COMM_WORLD, &

cinfo%zm*cinfo%ym*cinfo%xm*info%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,
&

cinfo%mz*cinfo%my*cinfo%mx*info%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,
&
                        PETSC_NULL,c%cA,ierr);
CALL MatShellSetOperation(c%cA,MATOP_MULT,matrixantiplanecoarsefv,ierr);
CALL MatAssemblyBegin(c%cA,MAT_FINAL_ASSEMBLY,ierr);
CALL MatAssemblyEnd(c%cA,MAT_FINAL_ASSEMBLY,ierr);
CALL KSPSetOperators(c%cksp,c%cA,c%cA,DIFFERENT_NONZERO_PATTERN,ierr);
CALL KSPSetUp(c%cksp,ierr);
CALL PCMGGetCoarseSolve(c%pc,c%cksp,ierr);


Then I set up the interpolation and restriction:

DO l=1,c%nlevels-1
       ...
       ! shell matrix MxN with M>N
       CALL MatCreateShell(PETSC_COMM_WORLD, &

cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,minfo%zm*minfo%ym*minfo%xm*minfo%dof,
&

cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,minfo%mz*minfo%my*minfo%mx*minfo%dof,
&
                           PETSC_NULL,c%interp(l),ierr);

       CALL MatShellSetOperation(c%interp(l),MATOP_MULT,matrixinterp,ierr);
       CALL MatAssemblyBegin(c%interp(l),MAT_FINAL_ASSEMBLY,ierr);
       CALL MatAssemblyEnd(c%interp(l),MAT_FINAL_ASSEMBLY,ierr);

       ! shell matrix MxN with M<N
       CALL MatCreateShell(PETSC_COMM_WORLD, &

minfo%zm*minfo%ym*minfo%xm*minfo%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,
&

minfo%mz*minfo%my*minfo%mx*minfo%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,
&
                           PETSC_NULL,c%restrict(l),ierr);

       CALL MatShellSetOperation(c%restrict(l),MATOP_MULT,matrixrestrict,ierr);
       CALL MatAssemblyBegin(c%restrict(l),MAT_FINAL_ASSEMBLY,ierr);
       CALL MatAssemblyEnd(c%restrict(l),MAT_FINAL_ASSEMBLY,ierr);

       ! define the intergrid transfer operations
       CALL PCMGSetInterpolation(c%pc,l,c%interp(l),ierr);
       CALL PCMGSetRestriction(c%pc,l,c%restrict(l),ierr);

END DO

Then i provide the smoother for each level:

DO l=0,c%nlevels-1
       ! global size at current level
       cinfo%mx=(info%mx/2**(c%nlevels-l))*2+1
       cinfo%my=(info%my/2**(c%nlevels-l))*2+1
       cinfo%mz=1
       ! locally-owned size at current level
       CALL getmxl(x,info%xm,c%nlevels-l-1,cinfo%xm)
       CALL getmxl(y,info%ym,c%nlevels-l-1,cinfo%ym)
       cinfo%zm=1
       cinfo%dof=info%dof

       !-----------------------------------------------------------------------------
       ! smoothing matrix
       !-----------------------------------------------------------------------------

       ! square smoothing matrix
       CALL MatCreateShell(PETSC_COMM_WORLD, &

cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,
&

cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,
&
                           PETSC_NULL,c%sA(l+1),ierr);

       CALL MatShellSetOperation(c%sA(l+1),MATOP_MULT,matrixsmoothfv,ierr);
       CALL MatAssemblyBegin(c%sA(l+1),MAT_FINAL_ASSEMBLY,ierr);
       CALL MatAssemblyEnd(c%sA(l+1),MAT_FINAL_ASSEMBLY,ierr);

       CALL PCMGGetSmoother(c%pc,l,sksp,ierr);
       CALL KSPCreate(PETSC_COMM_WORLD,sksp,ierr);
       CALL KSPSetOperators(sksp,c%sA(l+1),c%sA(l+1),DIFFERENT_NONZERO_PATTERN,ierr);
       CALL KSPSetUp(sksp,ierr);

       !-----------------------------------------------------------------------------
       ! residual operator
       !-----------------------------------------------------------------------------

       ! square forward-model matrix
       CALL MatCreateShell(PETSC_COMM_WORLD, &

cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,
&

cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,
&
                           PETSC_NULL,c%rA(l+1),ierr);

       CALL MatShellSetOperation(c%rA(l+1),MATOP_MULT,matrixresidualfv,ierr);
       CALL MatAssemblyBegin(c%rA(l+1),MAT_FINAL_ASSEMBLY,ierr);
       CALL MatAssemblyEnd(c%rA(l+1),MAT_FINAL_ASSEMBLY,ierr);

       ! compute formula b - Ax with matrix-free matrix
       CALL PCMGSetResidual(c%pc,l,matrixresidualfv,c%rA(l+1),ierr);

END DO

I also provide work space...
Despite that, I have a crash upon running KSPSetup() with "[0]PETSC
ERROR: PCSetUp_BJacobi()".


Also, I do no wish to provide a shell matrix for a direct solver at
the finest level but if I don't I get the error:

MatGetVecs() line 7265 in src/mat/interface/matrix.c
[0]PETSC ERROR: KSPGetVecs() line 806 in src/ksp/ksp/interface/iterativ.c
[0]PETSC ERROR: KSPSetUp_GMRES() line 94 in src/ksp/ksp/impls/gmres/gmres.c
[0]PETSC ERROR: KSPSetUp() line

So despite all my efforts, it looks that the solver still attempts to
do a direct solve.

I would appreciate your comments.
Best wishes,
Sylvain

2011/5/10 Barry Smith <bsmith at mcs.anl.gov>:
>
>
>   Barry
>
>
> On May 10, 2011, at 9:20 PM, Sylvain Barbot wrote:
>
>> Dear Jed,
>>
>> Thank you for your previous comments. I am still trying to design a
>> multigrid preconditionner for the Navier's equation of elasticity. I
>> have defined shell matrices for the interpolation, restriction,
>> smoothing, residual calculation and the fine-level solver (see
>> attached my_green.f90 file with stub functions). I thought that all
>> these shell matrices would be enough to do the job, but upon running
>> KSPSetup(), I get the following error:
>>
>> [0]PETSC ERROR: No support for this operation for this object type!
>> [0]PETSC ERROR: This matrix type, shell, does not support getting
>> diagonal block!
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 4, Fri Jul 30
>> 14:42:02 CDT 2010
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: PCSetUp_BJacobi() line 162 in src/ksp/pc/impls/bjacobi/bjacobi.c
>> [0]PETSC ERROR: PCSetUp() line 795 in src/ksp/pc/interface/precon.c
>> [0]PETSC ERROR: KSPSetUp() line 237 in src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: PCSetUp_MG() line 574 in src/ksp/pc/impls/mg/mg.c
>>
>> regarding your previous comments: "For Jacobi, you can have your
>> MatShell implement MatGetDiagonal and it will work. (...) If you want
>> to use a multiplicative relaxation like SOR, you would have to
>> implement it yourself.", I actually want to use SOR with my smoother
>> shell matrix. What am I doing wrong to make petsc still believe I want
>> to use Jacobi?
>>
>
>  By default, it is using block Jacobi on each level (except the coarsest). You need to tell it to use the relaxation you have coded. For example if your shell matrix provides a SOR/relaxation method then you need to use -mg_levels_pc_type sor  and it will call that on each level instead of block Jacobi. In looking at your code I see the only operation your provide your shell matrices are matrix-vector multiply. Therefore you cannot use SOR/relaxation. The only think you can use is GMRES (or some other Krylov method) with no preconditioner as the smoother. For this you can use -mg_levels_pc_type none -mg_levels_ksp_type gmres
>
>   In your message  you say  "I actually want to use SOR with my smoother" to do this you need to provide the smoother function and set it with MatShellSetOperations().
>
>
>> another related question: with a direct solver, one needs to assign a
>> shell matrix with a KSP object (KSPSetOperators); is it still required
>> when using multigrid? if I don't do so, I get the following error:
>>
>  Where are you using the direct solver for, on the entire linear system or on the coarse mesh? In either case you NEED to provide an explicit (not shell) sparse matrix otherwise a direct solver cannot be used. From the message below it looks like you did not provide a matrix somewhere. If you don't provide a matrix a direct solver cannot be used.
>
>
>   Barry
>
>> [0]PETSC ERROR: Null argument, when expecting valid pointer!
>> [0]PETSC ERROR: Null Object: Parameter # 1!
>> [0]PETSC ERROR: MatGetVecs() line 7265 in src/mat/interface/matrix.c
>> [0]PETSC ERROR: KSPGetVecs() line 806 in src/ksp/ksp/interface/iterativ.c
>> [0]PETSC ERROR: KSPSetUp_GMRES() line 94 in src/ksp/ksp/impls/gmres/gmres.c
>> [0]PETSC ERROR: KSPSetUp() line 199 in src/ksp/ksp/interface/itfunc.c
>>
>> however, I don't understand why I would need this matrix in a
>> multigrid framework.
>>
>> Is there any example out there of using PCMG with shell matrices?
>>
>> Thanks in advance for your response.
>> Best wishes,
>> Sylvain Barbot
>>
>>
>> 2010/10/7 Jed Brown <jed at 59a2.org>:
>>> On Thu, Oct 7, 2010 at 11:45, Sylvain Barbot <sylbar.vainbot at gmail.com>
>>> wrote:
>>>>
>>>> questions:
>>>> 1) is it at all possible to specify this mode of operation, from the
>>>> finest to the coarser level, and back? Any examples out there?
>>>
>>> -pc_mg_type multiplicative, or
>>>   PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
>>>
>>>>
>>>> 2) is it readily possible to use matrix shells with DMMG? I imagine
>>>> the Jacobian matrix may simply be provided as a matrix shell. Is there
>>>> any examples of multi-grid methods with shell matrices online?
>>>
>>> You can do this, but you will have to define a custom smoother.  For Jacobi,
>>> you can have your MatShell implement MatGetDiagonal and it will work.  I
>>> thought you could implement MatGetDiagonalBlock for PBJacobi, but it's not
>>> currently written that way (though it should be and that would be an easy
>>> change to make).  If you want to use a multiplicative relaxation like SOR,
>>> you would have to implement it yourself.  If you need something like ILU for
>>> a smoother, then you will have to pay for a matrix.  Note that one
>>> possibility is to assemble a matrix on all but the finest level so you can
>>> use stronger smoothers there, then make do with Jacobi on the finest level.
>>>
>>>>
>>>> 3) to deal with non-uniform sampling: can I provide the coordinates of
>>>> the finest grid with DASetCoordinates, then expect DMMG to provide the
>>>> subsampled coordinates at the coarser levels?
>>>
>>> Currently no, you have to set them on each level.  Perhaps you could do this
>>> by rolling a loop over levels and applying MatRestrict using the restriction
>>> matrix from DMMG (this might not be the sampling that you want).
>>> Jed
>> <my_green.f90>
>
>