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>
>
>
