[petsc-users] V-cycle multigrid with matrix shells
Barry Smith
bsmith at mcs.anl.gov
Thu May 12 08:13:37 CDT 2011
On May 12, 2011, at 1:01 AM, Sylvain Barbot wrote:
> 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".
>
Yes, I understand all this and it sounds fine.
> 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);
You say this is a smoother but you register it as a MATOP_MULT? What is the calling sequence of matrixsmoothfv() and what does it do? Smooth or multiply? MATOP_MULT is for registering a multiple MATOP_SOR (or MATOP_RELAX in older versions of PETSc) is for providing a smoother.
> 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()".
It is trying to use bjacobi because that is the default solver, you need to tell it to use some other solver explicitly if you do not want it to use block Jacobi. So if you registered your matrixsmoothfv() as a MATOP_SOR then you can use -mg_levels_pc_type sor to have it use that smoother. What solver options are you running with?
>
>
> 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.
>
This is not an indication that it is trying to use a direct solver, this is an indication then you never provided ANY kind of matrix at this level. You need to provide a matrix operator on each level for it to get work vectors from etc. It would be help to have the COMPLETE error message here so I can see exactly what is happening instead of a fragment making me guess what the problem is.
Barry
> 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>
>>
>>
More information about the petsc-users
mailing list