Actual source code: ex54f.F90

      1: !
      2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
      4: !  Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
      5: !
      6: !  This file is part of SLEPc.
      7: !  SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: !
     10: !  Program usage: mpiexec -n <np> ./ex54f [-help] [-n <n>] [all SLEPc options]
     11: !
     12: !  Description: Illustrates use of shell matrices in callback interface from Fortran.
     13: !  Similar to ex21.c. This one solves a simple diagonal linear eigenproblem as a NEP.
     14: !
     15: !  The command line options are:
     16: !    -n <n>, where <n> = matrix dimension

     18: #include <slepc/finclude/slepcnep.h>

     20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     21: ! Modules needed to pass and get the context to/from the Mat
     22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     24: module shell_ctx
     25:   use petscmat
     26:   implicit none
     27:   type :: MatCtx
     28:     PetscScalar :: lambda
     29:   end type MatCtx
     30: end module shell_ctx

     32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     33: ! Module used to implement the shell matrix operations
     34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     36: module ex54fmodule
     37:   use slepcnep
     38:   implicit none

     40: contains
     41:   ! --------------------------------------------------------------
     42:   ! MyNEPFunction - Computes Function matrix  T(lambda)
     43:   !
     44:   subroutine MyNEPFunction(nep, lambda, T, P, ctx, ierr)
     45:     use slepcnep
     46:     use shell_ctx
     47:     implicit none

     49:     NEP                   :: nep
     50:     PetscScalar           :: lambda
     51:     Mat                   :: T, P
     52:     PetscInt              :: ctx
     53:     type(MatCtx), pointer :: ctxT
     54:     PetscErrorCode, intent(out) :: ierr

     56:     PetscCall(MatShellGetContext(T, ctxT, ierr))
     57:     ctxT%lambda = lambda
     58:   end subroutine MyNEPFunction

     60:   ! --------------------------------------------------------------
     61:   ! MyNEPJacobian - Computes Jacobian matrix  T'(lambda)
     62:   !
     63:   subroutine MyNEPJacobian(nep, lambda, T, ctx, ierr)
     64:     use slepcnep
     65:     use shell_ctx
     66:     implicit none

     68:     NEP                   :: nep
     69:     PetscScalar           :: lambda
     70:     Mat                   :: T
     71:     PetscInt              :: ctx
     72:     type(MatCtx), pointer :: ctxT
     73:     PetscErrorCode, intent(out) :: ierr

     75:     PetscCall(MatShellGetContext(T, ctxT, ierr))
     76:     ctxT%lambda = lambda
     77:   end subroutine MyNEPJacobian

     79:   ! --------------------------------------------------------------
     80:   ! MatMult_A - Shell matrix operation, multiples y=A*x
     81:   ! Here A=(D-lambda*I) where D is a diagonal matrix
     82:   !
     83:   subroutine MatMult_A(A, x, y, ierr)
     84:     use shell_ctx
     85:     implicit none

     87:     Mat                   :: A
     88:     Vec                   :: x, y
     89:     PetscInt              :: i, istart, iend
     90:     PetscScalar           :: val
     91:     type(MatCtx), pointer :: ctxA
     92:     PetscErrorCode, intent(out) :: ierr
     93: !
     94:     PetscCall(VecGetOwnershipRange(x, istart, iend, ierr))
     95:     do i = istart, iend - 1
     96:       val = 1.0/real(i + 1, PETSC_REAL_KIND)
     97:       PetscCall(VecSetValue(y, i, val, INSERT_VALUES, ierr))
     98:     end do
     99:     PetscCall(VecAssemblyBegin(y, ierr))
    100:     PetscCall(VecAssemblyEnd(y, ierr))
    101:     PetscCall(VecPointwiseMult(y, y, x, ierr))
    102:     PetscCall(MatShellGetContext(A, ctxA, ierr))
    103:     PetscCall(VecAXPY(y, -ctxA%lambda, x, ierr))
    104:   end subroutine MatMult_A

    106:   ! --------------------------------------------------------------
    107:   ! MatDuplicate_A - Shell matrix operation, duplicates A
    108:   !
    109:   subroutine MatDuplicate_A(A, opt, M, ierr)
    110:     use shell_ctx
    111:     implicit none

    113:     Mat                   :: A, M
    114:     MatDuplicateOption    :: opt
    115:     PetscInt              :: ml, nl
    116:     type(MatCtx), pointer :: ctxM, ctxA
    117:     PetscErrorCode, intent(out) :: ierr

    119:     PetscCall(MatGetLocalSize(A, ml, nl, ierr))
    120:     PetscCall(MatShellGetContext(A, ctxA, ierr))
    121:     allocate (ctxM)
    122:     ctxM%lambda = ctxA%lambda
    123:     PetscCall(MatCreateShell(PETSC_COMM_WORLD, ml, nl, PETSC_DETERMINE, PETSC_DETERMINE, ctxM, M, ierr))
    124:     PetscCall(MatShellSetOperation(M, MATOP_MULT, MatMult_A, ierr))
    125:     PetscCall(MatShellSetOperation(M, MATOP_DESTROY, MatDestroy_A, ierr))
    126:   end subroutine MatDuplicate_A

    128:   ! --------------------------------------------------------------
    129:   ! MatDestroy_A - Shell matrix operation, destroys A
    130:   !
    131:   subroutine MatDestroy_A(A, ierr)
    132:     use shell_ctx
    133:     implicit none

    135:     Mat                   :: A
    136:     type(MatCtx), pointer :: ctxA
    137:     PetscErrorCode, intent(out) :: ierr

    139:     PetscCall(MatShellGetContext(A, ctxA, ierr))
    140:     deallocate (ctxA)
    141:   end subroutine MatDestroy_A

    143:   ! --------------------------------------------------------------
    144:   ! MatMult_B - Shell matrix operation, multiples y=B*x
    145:   ! Here B=-I
    146:   !
    147:   subroutine MatMult_B(B, x, y, ierr)
    148:     use petscmat
    149:     implicit none

    151:     Mat            :: B
    152:     Vec            :: x
    153:     Vec            :: y
    154:     PetscErrorCode, intent(out) :: ierr
    155:     PetscScalar, parameter :: mone = -1.0

    157:     PetscCall(VecCopy(x, y, ierr))
    158:     PetscCall(VecScale(y, mone, ierr))
    159:   end subroutine MatMult_B

    161: end module ex54fmodule

    163: !=================================================================================================

    165: program ex54f
    166:   use slepcnep
    167:   use shell_ctx
    168:   use ex54fmodule
    169:   implicit none

    171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172: ! Declarations
    173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    175:   NEP            :: nep
    176:   Mat            :: A, B
    177:   PetscInt       :: n = 400, nev = 3, nconv
    178:   PetscErrorCode :: ierr
    179:   PetscScalar    :: sigma
    180:   PetscBool      :: flg, terse
    181:   PetscMPIInt    :: rank
    182:   type(MatCtx)   :: ctxA, ctxB

    184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    185: ! Beginning of program
    186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    188:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
    189:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
    190:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    191:   if (rank == 0) then
    192:     write (*, '(/a,i4)') 'Nonlinear eigenproblem with shell matrices, n =', n
    193:   end if

    195: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    196: ! Create matrix-free operators for A and B corresponding to T and T'
    197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    199:   ctxA%lambda = 0.0
    200:   PetscCallA(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, ctxA, A, ierr))
    201:   PetscCallA(MatShellSetOperation(A, MATOP_MULT, MatMult_A, ierr))
    202:   PetscCallA(MatShellSetOperation(A, MATOP_DUPLICATE, MatDuplicate_A, ierr))

    204:   ctxB%lambda = 0.0   ! unused
    205:   PetscCallA(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, ctxB, B, ierr))
    206:   PetscCallA(MatShellSetOperation(B, MATOP_MULT, MatMult_B, ierr))

    208: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    209: ! Create the eigensolver and set various options
    210: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    212:   PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))
    213:   PetscCallA(NEPSetFunction(nep, A, A, MyNEPFunction, PETSC_NULL_INTEGER, ierr))
    214:   PetscCallA(NEPSetJacobian(nep, B, MyNEPJacobian, PETSC_NULL_INTEGER, ierr))
    215:   PetscCallA(NEPSetDimensions(nep, nev, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
    216:   sigma = 1.05
    217:   PetscCallA(NEPSetTarget(nep, sigma, ierr))
    218:   PetscCallA(NEPSetFromOptions(nep, ierr))

    220: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    221: ! Solve the eigensystem, display solution and clean up
    222: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    224:   PetscCallA(NEPSolve(nep, ierr))

    226:   PetscCallA(NEPGetConverged(nep, nconv, ierr))
    227:   if (rank == 0) then
    228:     write (*, '(a,i2/)') ' Number of converged approximate eigenpairs:', nconv
    229:   end if
    230: !
    231: ! ** show detailed info unless -terse option is given by user
    232:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
    233:   if (terse) then
    234:     PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    235:   else
    236:     PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
    237:     PetscCallA(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_WORLD, ierr))
    238:     PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
    239:     PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
    240:   end if

    242:   PetscCallA(NEPDestroy(nep, ierr))
    243:   PetscCallA(MatDestroy(A, ierr))
    244:   PetscCallA(MatDestroy(B, ierr))
    245:   PetscCallA(SlepcFinalize(ierr))

    247: end program ex54f

    249: !/*TEST
    250: !
    251: !   testset:
    252: !      args: -terse
    253: !      output_file: output/ex54f_1.out
    254: !      filter: grep -v approximate | sed -e "s/[+-]0\.0*i//g"
    255: !      test:
    256: !         suffix: 1_slp
    257: !         args: -nep_type slp -nep_slp_ksp_type gmres -nep_slp_pc_type none
    258: !         requires: double
    259: !      test:
    260: !         suffix: 1_nleigs
    261: !         args: -nep_type nleigs -rg_interval_endpoints 0.2,1.1 -nep_nleigs_ksp_type gmres -nep_nleigs_pc_type none
    262: !         requires: !complex
    263: !      test:
    264: !         suffix: 1_nleigs_complex
    265: !         args: -nep_type nleigs -rg_interval_endpoints 0.2,1.1,-.1,.1 -nep_nleigs_ksp_type gmres -nep_nleigs_pc_type none
    266: !         requires: complex
    267: !
    268: !TEST*/