[petsc-users] V-cycle multigrid with matrix shells
Barry Smith
bsmith at mcs.anl.gov
Wed May 18 13:41:54 CDT 2011
This is happening because it is trying to setup the solver on a particular level but the matrix has not yet been provided to that solver with KSPSetOperators()
0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Matrix must be set first!
------------------------------------------------------------------------
[0]PETSC ERROR: PCSetUp() line 775 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 602 in src/ksp/pc/impls/mg/mg.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
In looking at your code it seems you are trying to calling KSPSetOperators(). For some reason in the logic of the code it is not getting to the line that needs to set the operators.
I would run with the option -start_in_debugger noxterm and put a break in KSPSetFromOperators() and PCSetUp() then use "cont" and each time the debugger stops in KSPSetFromOperators() and PCSetUp() use "where" to see the current stack (where the code is). You can then match the KSPSetOperators() to each PCSetUp() and maybe figure out why the correct number of KSPSetOperators() are not being triggered in the code.
I would do it myself except I don't have all the code so cannot build it and run it.
Barry
On May 17, 2011, at 7:39 PM, Sylvain Barbot wrote:
> Barry,
>
> I have a call to KSPSetFromOptions() before KSPSetup() and KSPSolve()
> already. Is there any way I can force the default to be SOR smoother
> in the code instead of from the command line? With the command line
> you suggested, I get a slightly different error:
>
> [0]PETSC ERROR: --------------------- Error Message
> ------------------------------------
> [0]PETSC ERROR: Object is in wrong state!
> [0]PETSC ERROR: Matrix must be set first!
> [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: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: ./statici_petsc_1d_mg on a real4-wit named
> catalina.gps.caltech.edu by sbarbot Tue May 17 17:30:17 2011
> [0]PETSC ERROR: Libraries linked from
> /Users/sbarbot/Documents/work/petsc/real4-with-debug/lib
> [0]PETSC ERROR: Configure run at Thu Oct 7 19:12:09 2010
> [0]PETSC ERROR: Configure options
> --with-blas-lapack-dir=/opt/intel/Compiler/11.1/084/Frameworks/mkl/lib/em64t
> -with-mpi-dir=/shared/bin --with-debugging=1 --with-precision=single
> --with-petsc-arch=real4-with-debug
> [0]PETSC ERROR:
> ------------------------------------------------------------------------
> [0]PETSC ERROR: PCSetUp() line 775 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 602 in src/ksp/pc/impls/mg/mg.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
>
> there is no more reference to PCSetUp_BJacobi(), so it's a half
> success. I attach the code setting up the solver/preconditioner for
> reference (SUBROUTINE green_init_mg).
>
> btw, is there a way to better trace the location of the problem? i am
> using the __FUNCT__ pre-compiling definition, for instance
>
> #undef __FUNCT__
> #define __FUNCT__ "matrixinterp"
> SUBROUTINE matrixinterp(A,x,y,ierr)
>
> but it never shows up in the warning/error message.
>
> Thanks,
> Sylvain
>
>
> 2011/5/17 Barry Smith <bsmith at mcs.anl.gov>:
>>
>> This problem is because it is trying to use block Jacobi on a level as the smoother instead of SOR. You need to tell PCMG to use sor with -mg_levels_pc_type sor and make sure you have a KSPSetFromOptions() call before KSPSolve() so it catches those options.
>>
>> I'm sorry I forgot about the ourgetvecs() you did the correct fix. I will add that to petsc-dev also.
>>
>>
>> Barry
>>
>> The default is block Jacobi with ILU on each process.
>>
>>
>> On May 17, 2011, at 2:23 PM, Sylvain Barbot wrote:
>>
>>> Dear Barry,
>>>
>>> I made your suggested changes and recompiled the library. I also added
>>> the function ourview and ourgetvecs from petsc-dev and changed the
>>> last call of matshellsetoperations_ from
>>>
>>> PetscError(comm,__LINE__,"MatShellSetOperation_Fortran",__FILE__,__SDIR__,PETSC_ERR_ARG_WRONG,PETSC_ERROR_INITIAL,"Cannot
>>> set that matrix operation");
>>>
>>> to
>>>
>>> PetscError(__LINE__,"MatShellSetOperation_Fortran",__FILE__,__SDIR__,0,1,"Cannot
>>> set that matrix operation");
>>>
>>> for compatibility with my current version of Petsc.
>>>
>>> I changed the matrix type of my smoothing matrix to MAT_SOR. Upon
>>> running, i get the following error
>>>
>>> [0]PETSC ERROR: --------------------- Error Message
>>> ------------------------------------
>>> [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: See docs/changes/index.html for recent updates.
>>> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>>> [0]PETSC ERROR: See docs/index.html for manual pages.
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: ./statici_petsc_1d_mg on a real4-wit named
>>> catalina.gps.caltech.edu by sbarbot Tue May 17 12:08:10 2011
>>> [0]PETSC ERROR: Libraries linked from
>>> /Users/sbarbot/Documents/work/petsc/real4-with-debug/lib
>>> [0]PETSC ERROR: Configure run at Thu Oct 7 19:12:09 2010
>>> [0]PETSC ERROR: Configure options
>>> --with-blas-lapack-dir=/opt/intel/Compiler/11.1/084/Frameworks/mkl/lib/em64t
>>> -with-mpi-dir=/shared/bin --with-debugging=1 --with-precision=single
>>> --with-petsc-arch=real4-with-debug
>>> [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
>>> [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
>>>
>>> btw, when I call MatShellSetOperation, what should be the type of the
>>> shell matrix that performs the residuals, MAT_SOR as well, MAT_MULT or
>>> something else?
>>>
>>> thanks again for your help.
>>> Sylvain
>>>
>>> 2011/5/16 Barry Smith <bsmith at mcs.anl.gov>:
>>>>
>>>> On May 16, 2011, at 9:52 PM, Sylvain Barbot wrote:
>>>>
>>>>> Hi Barry,
>>>>>
>>>>> Thanks for your comments.
>>>>>
>>>>>>> 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.
>>>>>
>>>>> My function matrixsmoothfv indeed performs the smoothing operation so
>>>>> it not a multiplication. I can see no documentation about MATOP_SOR,
>>>>> is there an example somewhere online?
>>>>
>>>> For each possible Matrix function there is the name, for matrix vector multiply it is MATOP_MULT and for smoothing it is MATOP_SOR and the calling sequence you use in your function that does the operation is the same as the matrix operation. For example for MatMult() it is Mat,Vec,Vec,PetscErrorCode (see below for the smoother).
>>>>> Assigning MATOP_SOR to my
>>>>> smoothing matrix now generates the following runtime error message:
>>>>>
>>>>> [0]PETSC ERROR: MatShellSetOperation_Fortran() line 107 in
>>>>> src/mat/impls/shell/ftn-custom/zshellf.c
>>>>> application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
>>>>> rank 0 in job 37 catalina.gps.caltech.edu_50093 caused collective
>>>>> abort of all ranks
>>>>> exit status of rank 0: return code 1
>>>>
>>>> This is because we hadn't added support for providing the MATOP_SOR for the Fortran interface. Edit src/mat/imples/shell/ftn-custom/zshellf.c and locate the function matshellsetoperation_() and put instead of it the following code.
>>>>
>>>>
>>>> static PetscErrorCode oursor(Mat mat,Vec b,PetscReal omega,MatSORType flg,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
>>>> {
>>>> PetscErrorCode ierr = 0;
>>>> (*(PetscErrorCode (PETSC_STDCALL *)(Mat*,Vec*,PetscReal*,MatSORType*,PetscReal*,PetscInt*,PetscInt*,Vec*,PetscErrorCode*))(((PetscObject)mat)->fortran_func_pointers[9]))(&mat,&b,&omega,&flg,&shift,&its,&lits,&x,&ierr);
>>>> return ierr;
>>>> }
>>>>
>>>> void PETSC_STDCALL matshellsetoperation_(Mat *mat,MatOperation *op,PetscErrorCode (PETSC_STDCALL *f)(Mat*,Vec*,Vec*,PetscErrorCode*),PetscErrorCode *ierr)
>>>> {
>>>> MPI_Comm comm;
>>>>
>>>> *ierr = PetscObjectGetComm((PetscObject)*mat,&comm);if (*ierr) return;
>>>> PetscObjectAllocateFortranPointers(*mat,10);
>>>> if (*op == MATOP_MULT) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourmult);
>>>> ((PetscObject)*mat)->fortran_func_pointers[0] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_MULT_TRANSPOSE) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourmulttranspose);
>>>> ((PetscObject)*mat)->fortran_func_pointers[2] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_MULT_ADD) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourmultadd);
>>>> ((PetscObject)*mat)->fortran_func_pointers[1] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_MULT_TRANSPOSE_ADD) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourmulttransposeadd);
>>>> ((PetscObject)*mat)->fortran_func_pointers[3] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_GET_DIAGONAL) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourgetdiagonal);
>>>> ((PetscObject)*mat)->fortran_func_pointers[4] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_DIAGONAL_SCALE) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourdiagonalscale);
>>>> ((PetscObject)*mat)->fortran_func_pointers[5] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_DIAGONAL_SET) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourdiagonalset);
>>>> ((PetscObject)*mat)->fortran_func_pointers[6] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_GET_VECS) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourgetvecs);
>>>> ((PetscObject)*mat)->fortran_func_pointers[7] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_VIEW) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)ourview);
>>>> ((PetscObject)*mat)->fortran_func_pointers[8] = (PetscVoidFunction)f;
>>>> } else if (*op == MATOP_SOR) {
>>>> *ierr = MatShellSetOperation(*mat,*op,(PetscVoidFunction)oursor);
>>>> ((PetscObject)*mat)->fortran_func_pointers[9] = (PetscVoidFunction)f;
>>>> } else {
>>>> PetscError(comm,__LINE__,"MatShellSetOperation_Fortran",__FILE__,__SDIR__,PETSC_ERR_ARG_WRONG,PETSC_ERROR_INITIAL,
>>>> "Cannot set that matrix operation");
>>>> *ierr = 1;
>>>> }
>>>> }
>>>>
>>>>
>>>> Then in that directory run make to update the PETSc libraries.
>>>>
>>>> Now note the calling sequence of the relaxation function you need to provide. It is PetscErrorCode MatSOR(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x,PetscErrorCode ierr)
>>>> in your function you can ignore the omega, the flag, and the shift argument. You function should do its*lits smoothing steps.
>>>>
>>>> Once you have this going please let us know what happens.
>>>>
>>>> Barry
>>>>
>>>>
>>>>
>>>>
>>>>>
>>>>>>> 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?
>>>>>
>>>>> I setup the solver with:
>>>>>
>>>>> CALL PCMGSetType(c%pc,PC_MG_MULTIPLICATIVE,ierr);
>>>>>
>>>>> but this error might not be representative of the problem. upstream,
>>>>> my smoothing matrices were not set up properly.
>>>>>
>>>>>>> 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
>>>>>>>
>>>>>> 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.
>>>>
>>>>
>>
>>
> <green.f90>
More information about the petsc-users
mailing list