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