Dear petsc-/slepc-users,
I have been trying to understand matrix-free/shell matrices in PETSc for
eventual use in solving a non-linear eigenvalue problem using SLEPC. But
I seem to be having trouble with calls to MatShellGetContext. As far as
I understand, this function should initialize a pointer (second
argument) so that the subroutine output will point to the context
associated with my shell-matrix (let's say of TYPE(MatCtx))? When
calling that subroutine, I get the following error message:
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 2
In my code, the second input argument to the routine is a null-pointer
of TYPE(MatCtx),POINTER :: arg2. Which the error message appears to be
unhappy with. I noticed that there is no error message if I instead pass
an object TYPE(MatCtx) :: arg2 to the routine... which doesn't really
make sense to me? Could someone maybe explain to me what is going on, here?
Just in case, let me also attach my concrete example code (it is
supposed to be a Fortran version of the slepc-example in
slepc-3.8.1/src/nep/examples/tutorials/ex21.c). I have added an extra
call to MatShellGetContext on line 138, after the function and jacobian
should supposedly have been set up.
Thanks a lot for your help!
Cheers,
Samuel
! ----------------------------------------------------
! user-defined context
! ----------------------------------------------------
MODULE solver_context
#include "slepc/finclude/slepc.h"
! ------
USE slepcsys
USE slepceps
USE slepcnep
! ------
IMPLICIT NONE
EXTERNAL :: MatMult_Fun, MatGetDiagonal_Fun, MatDestroy_Fun, MatDuplicate_Fun
EXTERNAL :: MatMult_Jac, MatGetDiagonal_Jac, MatDestroy_Jac
TYPE :: ApplicationCtx
PetscScalar :: kappa
PetscReal :: h
END TYPE ApplicationCtx
TYPE :: MatCtx
PetscScalar :: lambda,kappa
PetscReal :: h
END TYPE MatCtx
END MODULE solver_context
MODULE solver_context_interfaces
USE solver_context
IMPLICIT NONE
! ----------------------------------------------------
INTERFACE MatShellSetContext
! --------------
SUBROUTINE MatShellSetContext(mat,ctx,ierr)
USE solver_context
Mat :: mat
TYPE(MatCtx) :: ctx
PetscErrorCode :: ierr
END SUBROUTINE MatShellSetContext
END INTERFACE MatShellSetContext
! ----------------------------------------------------
! ----------------------------------------------------
INTERFACE MatShellGetContext
! --------------
SUBROUTINE MatShellGetContext(mat,ctx,ierr)
USE solver_context
Mat :: mat
TYPE(MatCtx), POINTER :: ctx
PetscErrorCode :: ierr
END SUBROUTINE MatShellGetContext
END INTERFACE MatShellGetContext
! ----------------------------------------------------
END MODULE solver_context_interfaces
! ----------------------------------------------------
! main program
! ----------------------------------------------------
PROGRAM main
#include "slepc/finclude/slepc.h"
! ------
USE slepcsys
USE slepceps
USE slepcnep
! ------
USE solver_context
IMPLICIT NONE
CHARACTER(len=128) :: str
! --- matrix-free matrix evp (non-linear)
NEP :: nep
Mat :: F,J
TYPE(ApplicationCtx) :: ctx
TYPE(MatCtx) :: ctxF,ctxJ
TYPE(MatCtx), POINTER :: ctxF_pt
PetscInt :: n=128,nev
KSP :: ksp
PC :: pc
PetscMPIInt :: size
PetscBool :: terse, flg
PetscErrorCode :: ierr
! ---
EXTERNAL :: FormFunction
EXTERNAL :: FormJacobian
! initialize SLEPc
CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,size,ierr)
IF(size.NE.1) THEN
CALL MPI_ABORT(PETSC_COMM_WORLD,-1,ierr)
END IF
CALL PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-n",n,flg,ierr)
CALL PetscPrintf(PETSC_COMM_WORLD,"**************************\n",ierr)
CALL PetscPrintf(PETSC_COMM_WORLD,"1-D Nonlinear Eigenproblem\n",ierr)
WRITE(str,*) "n = ",n,"\n"
CALL PetscPrintf(PETSC_COMM_WORLD,TRIM(str),ierr)
CALL PetscPrintf(PETSC_COMM_WORLD,"**************************\n",ierr)
ctx%h = 1.0d0/n
ctx%kappa = 1.0;
! Create nonlinear eigensolver context
CALL NEPCreate(PETSC_COMM_WORLD,nep,ierr)
! ------------------------------------------------------------------------------
! Create matrix data structure; set Function evaluation routine
ctxF%h = ctx%h
ctxF%kappa = ctx%kappa
CALL PetscPrintf(PETSC_COMM_WORLD,"Function evaluation routine\n",ierr)
CALL MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,ctxF,F,ierr)
CALL MatShellSetOperation(F,MATOP_MULT,MatMult_Fun,ierr)
CALL MatShellSetOperation(F,MATOP_GET_DIAGONAL,MatGetDiagonal_Fun,ierr)
CALL MatShellSetOperation(F,MATOP_DESTROY,MatDestroy_Fun,ierr)
! Set function matrix data structure and default function evaluation routine
CALL NEPSetFunction(nep,F,F,FormFunction,ctx,ierr)
! ------------------------------------------------------------------------------
! ------------------------------------------------------------------------------
! Create Jacobian matrix data structure; set Jacobian evaluation routine
ctxJ%h = ctx%h
ctxJ%kappa = ctx%kappa
CALL PetscPrintf(PETSC_COMM_WORLD,"Jacobian evaluation routine\n",ierr)
CALL MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,ctxJ,J,ierr)
CALL MatShellSetOperation(J,MATOP_MULT,MatMult_Jac,ierr)
CALL MatShellSetOperation(J,MATOP_GET_DIAGONAL,MatGetDiagonal_Jac,ierr)
CALL MatShellSetOperation(J,MATOP_DESTROY,MatDestroy_Jac,ierr)
! Set function matrix data structure and default function evaluation routine
CALL NEPSetJacobian(nep,J,FormJacobian,ctx,ierr)
! ------------------------------------------------------------------------------
CALL MatShellGetContext(F,ctxF_pt,ierr)
PRINT*,'ctxF_pt%lambda: ',ctxF_pt%lambda
ctxF_pt%lambda = 11.
PRINT*,'ctxF_pt%lambda: ',ctxF_pt%lambda
CALL MatShellSetContext(F,ctxF_pt,ierr)
CALL MatShellGetContext(F,ctxF_pt,ierr)
PRINT*,'ctxF_pt%lambda: ',ctxF_pt%lambda
STOP
IF(ierr.NE.0) THEN
PRINT*,' '
PRINT*,'--------------------------------'
PRINT*,'********************************'
PRINT*,'Problem with MatShellGetContext.'
PRINT*,'Abort.'
PRINT*,'********************************'
PRINT*,'--------------------------------'
PRINT*,' '
CALL MPI_ABORT(MPI_COMM_WORLD,-1,ierr)
END IF
! ------------------------------------------------------------------------------
! Customize nonlinear solver; set runtime options
CALL PetscPrintf(PETSC_COMM_WORLD,"Customize nonlinear solver\n",ierr)
CALL NEPSetType(nep,NEPRII,ierr)
CALL NEPRIISetLagPreconditioner(nep,0,ierr)
CALL NEPRIIGetKSP(nep,ksp,ierr)
CALL KSPSetType(ksp,KSPBCGS,ierr)
CALL KSPGetPC(ksp,pc,ierr)
CALL PCSetType(pc,PCJACOBI,ierr)
! set solver parameters at runtime
CALL NEPSetFromOptions(nep,ierr)
! ------------------------------------------------------------------------------
! ------------------------------------------------------------------------------
! Solve the eigensystem
CALL PetscPrintf(PETSC_COMM_WORLD,"Solve eigensystem\n",ierr)
CALL NEPSolve(nep,ierr)
CALL NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
WRITE(str,*) "Number of requested eigenvalues: ", nev, "\n"
CALL PetscPrintf(PETSC_COMM_WORLD,TRIM(str),ierr)
! ------------------------------------------------------------------------------
! ------------------------------------------------------------------------------
! Display solution and clean up
! ------------------------------------------------------------------------------
! show detailed info unless -terse option is given by user
CALL PetscPrintf(PETSC_COMM_WORLD,"Display solution and clean up\n",ierr)
CALL PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr)
CALL NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr)
CALL NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr)
CALL PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr)
CALL NEPDestroy(nep,ierr)
CALL MatDestroy(F,ierr)
CALL MatDestroy(J,ierr)
! finalize SLEPc
CALL SlepcFinalize(ierr)
END PROGRAM main
SUBROUTINE FormInitialGuess(x,ierr)
USE solver_context
IMPLICIT NONE
Vec :: x
PetscErrorCode :: ierr
!
PetscScalar :: one
!
one = 1.0d0
CALL VecSet(x,one,ierr)
!
END SUBROUTINE FormInitialGuess
SUBROUTINE FormFunction(nep,lambda,fun,B,ctx)
USE solver_context
IMPLICIT NONE
NEP :: nep
PetscScalar :: lambda
Mat :: fun,B
TYPE(MatCtx) :: ctx
! ----
TYPE(MatCtx),POINTER :: ctxF
PetscErrorCode :: ierr
!
CALL MatShellGetContext(fun,ctxF,ierr)
ctxF%lambda = lambda
!
END SUBROUTINE FormFunction
SUBROUTINE FormJacobian(nep,lambda,jac,ctx)
USE solver_context
IMPLICIT NONE
NEP :: nep
PetscScalar :: lambda
Mat :: jac
TYPE(MatCtx) :: ctx
! ----
TYPE(MatCtx),POINTER :: ctxJ
PetscErrorCode :: ierr
!
CALL MatShellGetContext(jac,ctxJ,ierr)
ctxJ%lambda = lambda
!
END SUBROUTINE FormJacobian
! ----------------------------------------------
SUBROUTINE MatMult_Fun(A,x,y,ierr)
USE solver_context
IMPLICIT NONE
Mat :: A
Vec :: x,y
! ----
PetscErrorCode :: ierr
TYPE(MatCtx),POINTER :: ctx
PetscInt :: i,n
PetscScalar, POINTER :: px(:)
PetscScalar, POINTER :: py(:)
PetscScalar :: c,d,de,oe
PetscReal :: h
PRINT*,'MatMult_Fun'
!
CALL MatShellGetContext(A,ctx,ierr)
CALL VecGetArrayReadF90(x,px,ierr)
CALL VecGetArrayF90(y,py,ierr)
!
CALL VecGetSize(x,n,ierr)
h = ctx%h
c = ctx%kappa/(ctx%lambda-ctx%kappa)
d = n
de = 2.0*(d-ctx%lambda*h/6.0)
oe = -d-ctx%lambda*h/6.0
py(1) = de*px(1) + oe*px(2)
DO i=2,n-1
py(i) = oe*px(i-1) + de*px(i) + oe*px(i+1)
END DO
de = d-ctx%lambda*h/3.0+ctx%lambda*c
py(n) = oe*px(n-1) + de*px(n)
!
CALL VecRestoreArrayReadF90(x,px,ierr)
CALL VecRestoreArray(y,py,ierr)
END SUBROUTINE MatMult_Fun
! ----------------------------------------------
SUBROUTINE MatGetDiagonal_Fun(A,diag,ierr)
USE solver_context
IMPLICIT NONE
Mat :: A
Vec :: diag
PetscErrorCode :: ierr
!
TYPE(MatCtx),POINTER :: ctx
PetscInt :: n
PetscScalar,POINTER :: pd(:)
PetscScalar :: c,d
PetscReal :: h
PRINT*,'MatGetDiagonal_Fun'
!
CALL MatShellGetContext(A,ctx,ierr)
CALL VecGetSize(diag,n,ierr)
h = ctx%h;
c = ctx%kappa/(ctx%lambda-ctx%kappa);
d = n;
CALL VecSet(diag,2.0*(d-ctx%lambda*h/3.0),ierr)
CALL VecGetArrayF90(diag,pd,ierr)
pd(n) = d-ctx%lambda*h/3.0+ctx%lambda*c;
CALL VecRestoreArrayF90(diag,pd,ierr)
END SUBROUTINE MatGetDiagonal_Fun
! ----------------------------------------------
SUBROUTINE MatDestroy_Fun(A,ierr)
USE solver_context
IMPLICIT NONE
Mat :: A
PetscErrorCode :: ierr
!
PRINT*,'MatDestroy_Fun'
!
END SUBROUTINE MatDestroy_Fun
! ----------------------------------------------
SUBROUTINE MatMult_Jac(A,x,y,ierr)
USE solver_context
IMPLICIT NONE
Mat :: A
Vec :: x,y
PetscErrorCode :: ierr
! ----
TYPE(MatCtx),POINTER :: ctx
PetscInt :: i,n
PetscScalar, POINTER :: px(:)
PetscScalar, POINTER :: py(:)
PetscScalar :: c,d,de,oe
PetscReal :: h
PRINT*,'MatMul_Jac'
!
CALL MatShellGetContext(A,ctx,ierr)
CALL VecGetArrayReadF90(x,px,ierr)
CALL VecGetArrayF90(y,py,ierr)
!
CALL VecGetSize(x,n,ierr)
h = ctx%h;
c = ctx%kappa/(ctx%lambda-ctx%kappa);
de = -2.0*h/3.0; ! /* diagonal entry */
oe = -h/6.0; ! /* offdiagonal entry */
py(1) = de*px(1) + oe*px(1);
DO i=2,n-1
py(i) = oe*px(i-1) + de*px(i) + oe*px(i+1)
END DO
de = -h/3.0-c*c; !/* diagonal entry of last row */
py(n) = oe*px(n-1) + de*px(n);
!
CALL VecRestoreArrayReadF90(x,px,ierr)
CALL VecRestoreArrayF90(y,py,ierr)
!
END SUBROUTINE MatMult_Jac
! ----------------------------------------------
SUBROUTINE MatGetDiagonal_Jac(A,diag,ierr)
USE solver_context
IMPLICIT NONE
Mat :: A
Vec :: diag
PetscErrorCode :: ierr
!
TYPE(MatCtx),POINTER :: ctx
PetscInt :: n
PetscScalar,POINTER :: pd(:)
PetscScalar :: c,petsc_val
PetscReal :: h
PRINT*,'Get diagonal_Jac'
CALL MatShellGetContext(A,ctx,ierr)
CALL VecGetSize(diag,n,ierr)
h = ctx%h;
c = ctx%kappa/(ctx%lambda-ctx%kappa);
petsc_val = -2.0*h/3.0
CALL VecSet(diag,petsc_val,ierr)
CALL VecGetArray(diag,pd,ierr)
pd(n) = -h/3.0-c*c;
CALL VecRestoreArray(diag,pd,ierr)
END SUBROUTINE MatGetDiagonal_Jac
! ----------------------------------------------
SUBROUTINE MatDestroy_Jac(A,ierr)
USE solver_context
IMPLICIT NONE
Mat :: A
PetscErrorCode :: ierr
END SUBROUTINE MatDestroy_Jac
