Actual source code: ex10f.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> ./ex10f [-help] [-n <n>] [all SLEPc options]
     11: !
     12: !  Description: Illustrates the use of shell spectral transformations.
     13: !  The problem to be solved is the same as ex1.c and corresponds to the
     14: !  Laplacian operator in 1 dimension
     15: !
     16: !  The command line options are:
     17: !    nm <n>, where <n> is the number of grid subdivisions = matrix dimension
     18: ! ----------------------------------------------------------------------

     20: ! Module contains data needed by shell ST
     21: !
     22: #include <slepc/finclude/slepceps.h>
     23: module ex10fmodule
     24:   use slepceps
     25:   implicit none

     27:   KSP myksp

     29: contains
     30:   ! -------------------------------------------------------------------
     31:   ! STApply_User - This routine demonstrates the use of a user-provided spectral
     32:   ! transformation. The transformation implemented in this code is just OP=A^-1.
     33:   !
     34:   ! Input Parameters:
     35:   !   st - spectral transformation context
     36:   !   x - input vector
     37:   !
     38:   ! Output Parameter:
     39:   !   y - output vector
     40:   !
     41:   subroutine STApply_User(st, x, y, ierr)
     42:     use slepceps
     43:     implicit none

     45:     ST             :: st
     46:     Vec            :: x, y
     47:     PetscErrorCode, intent(out) :: ierr

     49:     PetscCall(KSPSolve(myksp, x, y, ierr))
     50:   end subroutine

     52:   ! -------------------------------------------------------------------
     53:   ! STApplyTranspose_User - This is not required unless using a two-sided eigensolver
     54:   !
     55:   ! Input Parameters:
     56:   !   st - spectral transformation context
     57:   !   x - input vector
     58:   !
     59:   ! Output Parameter:
     60:   !   y - output vector
     61:   !
     62:   subroutine STApplyTranspose_User(st, x, y, ierr)
     63:     use slepceps
     64:     implicit none

     66:     ST             :: st
     67:     Vec            :: x, y
     68:     PetscErrorCode, intent(out) :: ierr

     70:     PetscCall(KSPSolveTranspose(myksp, x, y, ierr))
     71:   end subroutine

     73: #if defined(PETSC_USE_COMPLEX)
     74:   ! -------------------------------------------------------------------
     75:   ! STApplyHermitianTranspose_User - This is not required unless using a two-sided eigensolver
     76:   ! in complex scalars
     77:   !
     78:   ! Input Parameters:
     79:   !   st - spectral transformation context
     80:   !   x - input vector
     81:   !
     82:   ! Output Parameter:
     83:   !   y - output vector
     84:   !
     85:   subroutine STApplyHermitianTranspose_User(st, x, y, ierr)
     86:     use slepceps
     87:     implicit none

     89:     ST             :: st
     90:     Vec            :: x, y, w
     91:     PetscErrorCode, intent(out) :: ierr

     93:     PetscCall(VecDuplicate(x, w, ierr))
     94:     PetscCall(VecCopy(x, w, ierr))
     95:     PetscCall(VecConjugate(w, ierr))
     96:     PetscCall(KSPSolveTranspose(myksp, w, y, ierr))
     97:     PetscCall(VecConjugate(y, ierr))
     98:     PetscCall(VecDestroy(w, ierr))
     99:   end subroutine
    100: #endif

    102:   ! -------------------------------------------------------------------
    103:   ! STBackTransform_User - This routine demonstrates the use of a user-provided spectral
    104:   ! transformation
    105:   !
    106:   ! Input Parameters:
    107:   !   st - spectral transformation context
    108:   !   n  - number of eigenvalues to transform
    109:   !
    110:   ! Output Parameters:
    111:   !   eigr - real part of eigenvalues
    112:   !   eigi - imaginary part of eigenvalues
    113:   !
    114:   subroutine STBackTransform_User(st, n, eigr, eigi, ierr)
    115:     use slepceps
    116:     implicit none

    118:     ST             :: st
    119:     PetscInt       :: n, j
    120:     PetscScalar    :: eigr(*), eigi(*)
    121:     PetscErrorCode, intent(out) :: ierr

    123:     do j = 1, n
    124:       eigr(j) = 1.0/eigr(j)
    125:     end do
    126:     ierr = 0
    127:   end subroutine

    129: end module ex10fmodule

    131: ! ----------------------------------------------------------------------

    133: program ex10f
    134:   use slepceps
    135:   use ex10fmodule
    136:   implicit none

    138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    139: ! Declarations
    140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    141: !
    142:   Mat            :: A    ! operator matrix
    143:   EPS            :: eps  ! eigenproblem solver context
    144:   ST             :: st
    145:   EPSType        :: tname
    146:   PetscInt       :: n, i, Istart, Iend
    147:   PetscInt       :: nev, row(1), col(3)
    148:   PetscScalar    :: val(3)
    149:   PetscBool      :: flg, isShell, terse
    150:   PetscMPIInt    :: rank
    151:   PetscErrorCode :: ierr

    153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    154: ! Beginning of program
    155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    157:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
    158:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    159:   n = 30
    160:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))

    162:   if (rank == 0) then
    163:     write (*, '(/a,i6/)') '1-D Laplacian Eigenproblem (shell-enabled), n=', n
    164:   end if

    166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    167: ! Compute the operator matrix that defines the eigensystem, Ax=kx
    168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    170:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
    171:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    172:   PetscCallA(MatSetFromOptions(A, ierr))

    174:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
    175:   if (Istart == 0) then
    176:     row(1) = 0
    177:     col(1) = 0
    178:     col(2) = 1
    179:     val(1) = 2.0
    180:     val(2) = -1.0
    181:     PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, row, 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
    182:     Istart = Istart + 1
    183:   end if
    184:   if (Iend == n) then
    185:     row(1) = n - 1
    186:     col(1) = n - 2
    187:     col(2) = n - 1
    188:     val(1) = -1.0
    189:     val(2) = 2.0
    190:     PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, row, 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
    191:     Iend = Iend - 1
    192:   end if
    193:   val(1) = -1.0
    194:   val(2) = 2.0
    195:   val(3) = -1.0
    196:   do i = Istart, Iend - 1
    197:     row(1) = i
    198:     col(1) = i - 1
    199:     col(2) = i
    200:     col(3) = i + 1
    201:     PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, row, 3_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
    202:   end do

    204:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
    205:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

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

    211: ! ** Create eigensolver context
    212:   PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))

    214: ! ** Set operators. In this case, it is a standard eigenvalue problem
    215:   PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
    216:   PetscCallA(EPSSetProblemType(eps, EPS_NHEP, ierr))

    218: ! ** Set solver parameters at runtime
    219:   PetscCallA(EPSSetFromOptions(eps, ierr))

    221: ! ** Initialize shell spectral transformation if selected by user
    222:   PetscCallA(EPSGetST(eps, st, ierr))
    223:   PetscCallA(PetscObjectTypeCompare(st, STSHELL, isShell, ierr))

    225:   if (isShell) then
    226: !   ** Change sorting criterion since this ST example computes values
    227: !   ** closest to 0
    228:     PetscCallA(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL, ierr))

    230: !   ** In Fortran, instead of a context for the user-defined spectral transform
    231: !   ** we use a module containing any application-specific data, initialized here
    232:     PetscCallA(KSPCreate(PETSC_COMM_WORLD, myksp, ierr))
    233:     PetscCallA(KSPAppendOptionsPrefix(myksp, "st_", ierr))

    235: !   ** (Required) Set the user-defined routine for applying the operator
    236:     PetscCallA(STShellSetApply(st, STApply_User, ierr))

    238: !   ** (Optional) Set the user-defined routine for applying the transposed operator
    239:     PetscCallA(STShellSetApplyTranspose(st, STApplyTranspose_User, ierr))

    241: #if defined(PETSC_USE_COMPLEX)
    242: !   ** (Optional) Set the user-defined routine for applying the conjugate-transposed operator
    243:     PetscCallA(STShellSetApplyHermitianTranspose(st, STApplyHermitianTranspose_User, ierr))
    244: #endif

    246: !   ** (Optional) Set the user-defined routine for back-transformation
    247:     PetscCallA(STShellSetBackTransform(st, STBackTransform_User, ierr))

    249: !   ** (Optional) Set a name for the transformation, used for STView()
    250:     PetscCallA(PetscObjectSetName(st, 'MyTransformation', ierr))

    252: !   ** (Optional) Do any setup required for the new transformation
    253:     PetscCallA(KSPSetOperators(myksp, A, A, ierr))
    254:     PetscCallA(KSPSetFromOptions(myksp, ierr))
    255:   end if

    257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    258: ! Solve the eigensystem
    259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    261:   PetscCallA(EPSSolve(eps, ierr))

    263: ! ** Optional: Get some information from the solver and display it
    264:   PetscCallA(EPSGetType(eps, tname, ierr))
    265:   if (rank == 0) then
    266:     write (*, '(a,a/)') ' Solution method: ', tname
    267:   end if
    268:   PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    269:   if (rank == 0) then
    270:     write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
    271:   end if

    273: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    274: ! Display solution and clean up
    275: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    277: ! ** show detailed info unless -terse option is given by user
    278:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
    279:   if (terse) then
    280:     PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    281:   else
    282:     PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
    283:     PetscCallA(EPSConvergedReasonView(eps, PETSC_VIEWER_STDOUT_WORLD, ierr))
    284:     PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
    285:     PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
    286:   end if
    287:   if (isShell) then
    288:     PetscCallA(KSPDestroy(myksp, ierr))
    289:   end if
    290:   PetscCallA(EPSDestroy(eps, ierr))
    291:   PetscCallA(MatDestroy(A, ierr))
    292:   PetscCallA(SlepcFinalize(ierr))
    293: end program ex10f

    295: !/*TEST
    296: !
    297: !   testset:
    298: !      args: -eps_nev 5 -eps_non_hermitian -terse
    299: !      output_file: output/ex10_1.out
    300: !      requires: !single
    301: !      test:
    302: !         suffix: 1_sinvert
    303: !         args: -st_type sinvert
    304: !      test:
    305: !         suffix: 1_shell
    306: !         args: -st_type shell -eps_two_sided {{0 1}}
    307: !
    308: !TEST*/