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.
>>>
>>>
>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: green.f90
Type: application/octet-stream
Size: 45011 bytes
Desc: not available
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110517/8aeab216/attachment-0001.obj>
