Actual source code: ex27f.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> ./ex27f [-help] [-n <n>] [all SLEPc options]
     11: !
     12: !  Description: Simple NLEIGS example. Fortran90 equivalent of ex27.c
     13: !
     14: !  The command line options are:
     15: !    -n <n>, where <n> = matrix dimension
     16: !
     17: ! ----------------------------------------------------------------------
     18: !   Solve T(lambda)x=0 using NLEIGS solver
     19: !      with T(lambda) = -D+sqrt(lambda)*I
     20: !      where D is the Laplacian operator in 1 dimension
     21: !      and with the interpolation interval [.01,16]
     22: ! ----------------------------------------------------------------------
     23: !
     24: #include <slepc/finclude/slepcnep.h>

     26: module ex27fmodule
     27:   use slepcnep
     28:   implicit none

     30: contains
     31:   ! --------------------------------------------------------------
     32:   ! FormFunction - Computes Function matrix  T(lambda)
     33:   !
     34:   subroutine FormFunction(nep, lambda, fun, B, ctx, ierr)
     35:     use slepcnep
     36:     implicit none

     38:     NEP            :: nep
     39:     PetscScalar    :: lambda, val(0:2), t
     40:     Mat            :: fun, B
     41:     PetscInt       :: ctx, i, n, col(0:2), Istart, Iend, Istart0, Iend0
     42:     PetscBool      :: FirstBlock = PETSC_FALSE, LastBlock = PETSC_FALSE
     43:     PetscErrorCode, intent(out) :: ierr

     45:     ! ** Compute Function entries and insert into matrix
     46:     t = sqrt(lambda)
     47:     PetscCall(MatGetSize(fun, n, PETSC_NULL_INTEGER, ierr))
     48:     PetscCall(MatGetOwnershipRange(fun, Istart, Iend, ierr))
     49:     if (Istart == 0) FirstBlock = PETSC_TRUE
     50:     if (Iend == n) LastBlock = PETSC_TRUE
     51:     val(0) = 1.0
     52:     val(1) = t - 2.0
     53:     val(2) = 1.0

     55:     Istart0 = Istart
     56:     if (FirstBlock) Istart0 = Istart + 1
     57:     Iend0 = Iend
     58:     if (LastBlock) Iend0 = Iend - 1

     60:     do i = Istart0, Iend0 - 1
     61:       col(0) = i - 1
     62:       col(1) = i
     63:       col(2) = i + 1
     64:       PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
     65:     end do

     67:     if (LastBlock) then
     68:       i = n - 1
     69:       col(0) = n - 2
     70:       col(1) = n - 1
     71:       val(0) = 1.0
     72:       val(1) = t - 2.0
     73:       PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
     74:     end if

     76:     if (FirstBlock) then
     77:       i = 0
     78:       col(0) = 0
     79:       col(1) = 1
     80:       val(0) = t - 2.0
     81:       val(1) = 1.0
     82:       PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
     83:     end if

     85:     ! ** Assemble matrix
     86:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
     87:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
     88:     PetscCall(MatAssemblyBegin(fun, MAT_FINAL_ASSEMBLY, ierr))
     89:     PetscCall(MatAssemblyEnd(fun, MAT_FINAL_ASSEMBLY, ierr))

     91:   end subroutine FormFunction

     93:   ! --------------------------------------------------------------
     94:   ! FormJacobian - Computes Jacobian matrix  T'(lambda)
     95:   !
     96:   subroutine FormJacobian(nep, lambda, jac, ctx, ierr)
     97:     use slepcnep
     98:     implicit none

    100:     NEP            :: nep
    101:     PetscScalar    :: lambda, t
    102:     Mat            :: jac
    103:     PetscInt       :: ctx
    104:     Vec            :: d
    105:     PetscErrorCode, intent(out) :: ierr

    107:     PetscCall(MatCreateVecs(jac, d, PETSC_NULL_VEC, ierr))
    108:     t = 0.5/sqrt(lambda)
    109:     PetscCall(VecSet(d, t, ierr))
    110:     PetscCall(MatDiagonalSet(jac, d, INSERT_VALUES, ierr))
    111:     PetscCall(VecDestroy(d, ierr))

    113:   end subroutine FormJacobian

    115:   ! --------------------------------------------------------------
    116:   !  ComputeSingularities - This is a user-defined routine to compute maxnp
    117:   !  points (at most) in the complex plane where the function T(.) is not analytic.
    118:   !
    119:   !  In this case, we discretize the singularity region (-inf,0)~(-1e+6,-1e-5)
    120:   !
    121:   !  Input Parameters:
    122:   !    nep   - nonlinear eigensolver context
    123:   !    maxnp - on input number of requested points in the discretization (can be set)
    124:   !    xi    - computed values of the discretization
    125:   !    dummy - optional user-defined monitor context (unused here)
    126:   !
    127:   subroutine ComputeSingularities(nep, maxnp, xi, dummy, ierr)
    128:     use slepcnep
    129:     implicit none

    131:     NEP            :: nep
    132:     PetscInt       :: maxnp, dummy
    133:     PetscScalar    :: xi(0:maxnp - 1)
    134:     PetscReal      :: h
    135:     PetscInt       :: i
    136:     PetscErrorCode, intent(out) :: ierr

    138:     h = 11.0/real(maxnp - 1, PETSC_REAL_KIND)
    139:     xi(0) = -1e-5
    140:     xi(maxnp - 1) = -1e+6
    141:     do i = 1, maxnp - 2
    142:       xi(i) = -10**(-5 + h*i)
    143:     end do
    144:     ierr = 0

    146:   end subroutine ComputeSingularities

    148: end module ex27fmodule

    150: program ex27f
    151:   use slepcnep
    152:   use ex27fmodule
    153:   implicit none

    155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    156: ! Declarations
    157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    159:   NEP                :: nep
    160:   Mat                :: A(2), F, J
    161:   NEPType            :: ntype
    162:   PetscInt           :: n = 100, nev, Istart, Iend, i, col
    163:   PetscErrorCode     :: ierr
    164:   PetscBool          :: terse, flg, split = PETSC_TRUE
    165:   PetscReal          :: ia, ib, ic, id
    166:   RG                 :: rg
    167:   FN                 :: fn(2)
    168:   PetscScalar        :: coeffs, sigma
    169:   character(len=128) :: string
    170:   PetscScalar, parameter :: one = 1.0

    172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    173: ! Beginning of program
    174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    176:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
    177:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-n", n, flg, ierr))
    178:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-split", split, flg, ierr))
    179:   if (split) then
    180:     write (string, *) 'Square root eigenproblem, n=', n, ' (split-form)\n'
    181:   else
    182:     write (string, *) 'Square root eigenproblem, n=', n, '\n'
    183:   end if
    184:   PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(string), ierr))

    186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    187: ! Create nonlinear eigensolver context and set options
    188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    190:   PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))
    191:   PetscCallA(NEPSetType(nep, NEPNLEIGS, ierr))
    192:   PetscCallA(NEPNLEIGSSetSingularitiesFunction(nep, ComputeSingularities, 0, ierr))
    193:   PetscCallA(NEPGetRG(nep, rg, ierr))
    194:   PetscCallA(RGSetType(rg, RGINTERVAL, ierr))
    195:   ia = 0.01
    196:   ib = 16.0
    197: #if defined(PETSC_USE_COMPLEX)
    198:   ic = -0.001
    199:   id = 0.001
    200: #else
    201:   ic = 0.0
    202:   id = 0.0
    203: #endif
    204:   PetscCallA(RGIntervalSetEndpoints(rg, ia, ib, ic, id, ierr))
    205:   sigma = 1.1
    206:   PetscCallA(NEPSetTarget(nep, sigma, ierr))

    208: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    209: ! Define the nonlinear problem
    210: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    212:   if (split) then
    213:     ! ** Create matrices for the split form
    214:     PetscCallA(MatCreate(PETSC_COMM_WORLD, A(1), ierr))
    215:     PetscCallA(MatSetSizes(A(1), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    216:     PetscCallA(MatSetFromOptions(A(1), ierr))
    217:     PetscCallA(MatGetOwnershipRange(A(1), Istart, Iend, ierr))
    218:     coeffs = -2.0
    219:     do i = Istart, Iend - 1
    220:       if (i > 0) then
    221:         col = i - 1
    222:         PetscCallA(MatSetValue(A(1), i, col, one, INSERT_VALUES, ierr))
    223:       end if
    224:       if (i < n - 1) then
    225:         col = i + 1
    226:         PetscCallA(MatSetValue(A(1), i, col, one, INSERT_VALUES, ierr))
    227:       end if
    228:       PetscCallA(MatSetValue(A(1), i, i, coeffs, INSERT_VALUES, ierr))
    229:     end do
    230:     PetscCallA(MatAssemblyBegin(A(1), MAT_FINAL_ASSEMBLY, ierr))
    231:     PetscCallA(MatAssemblyEnd(A(1), MAT_FINAL_ASSEMBLY, ierr))

    233:     PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, one, A(2), ierr))

    235:     ! ** Define functions for the split form
    236:     PetscCallA(FNCreate(PETSC_COMM_WORLD, fn(1), ierr))
    237:     PetscCallA(FNSetType(fn(1), FNRATIONAL, ierr))
    238:     PetscCallA(FNRationalSetNumerator(fn(1), 1_PETSC_INT_KIND, [one], ierr))
    239:     PetscCallA(FNCreate(PETSC_COMM_WORLD, fn(2), ierr))
    240:     PetscCallA(FNSetType(fn(2), FNSQRT, ierr))
    241:     PetscCallA(NEPSetSplitOperator(nep, 2_PETSC_INT_KIND, A, fn, SUBSET_NONZERO_PATTERN, ierr))
    242:   else
    243:     ! ** Callback form: create matrix and set Function evaluation routine
    244:     PetscCallA(MatCreate(PETSC_COMM_WORLD, F, ierr))
    245:     PetscCallA(MatSetSizes(F, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    246:     PetscCallA(MatSetFromOptions(F, ierr))
    247:     PetscCallA(MatSeqAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    248:     PetscCallA(MatMPIAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    249:     PetscCallA(NEPSetFunction(nep, F, F, FormFunction, PETSC_NULL_INTEGER, ierr))

    251:     PetscCallA(MatCreate(PETSC_COMM_WORLD, J, ierr))
    252:     PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    253:     PetscCallA(MatSetFromOptions(J, ierr))
    254:     PetscCallA(MatSeqAIJSetPreallocation(J, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    255:     PetscCallA(MatMPIAIJSetPreallocation(J, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    256:     PetscCallA(NEPSetJacobian(nep, J, FormJacobian, PETSC_NULL_INTEGER, ierr))
    257:   end if

    259:   PetscCallA(NEPSetFromOptions(nep, ierr))

    261: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    262: ! Solve the eigensystem
    263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    265:   PetscCallA(NEPSolve(nep, ierr))
    266:   PetscCallA(NEPGetType(nep, ntype, ierr))
    267:   write (string, *) 'Solution method: ', ntype, '\n'
    268:   PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(string), ierr))
    269:   PetscCallA(NEPGetDimensions(nep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    270:   write (string, *) 'Number of requested eigenvalues:', nev, '\n'
    271:   PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(string), ierr))

    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(NEPErrorView(nep, NEP_ERROR_BACKWARD, PETSC_NULL_VIEWER, ierr))
    281:   else
    282:     PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
    283:     PetscCallA(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_WORLD, ierr))
    284:     PetscCallA(NEPErrorView(nep, NEP_ERROR_BACKWARD, PETSC_VIEWER_STDOUT_WORLD, ierr))
    285:     PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
    286:   end if

    288:   if (split) then
    289:     PetscCallA(MatDestroy(A(1), ierr))
    290:     PetscCallA(MatDestroy(A(2), ierr))
    291:     PetscCallA(FNDestroy(fn(1), ierr))
    292:     PetscCallA(FNDestroy(fn(2), ierr))
    293:   else
    294:     PetscCallA(MatDestroy(F, ierr))
    295:     PetscCallA(MatDestroy(J, ierr))
    296:   end if
    297:   PetscCallA(NEPDestroy(nep, ierr))
    298:   PetscCallA(SlepcFinalize(ierr))

    300: end program ex27f

    302: !/*TEST
    303: !
    304: !   test:
    305: !      suffix: 1
    306: !      args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
    307: !      requires: !single
    308: !      filter: sed -e "s/[+-]0\.0*i//g"
    309: !
    310: !   test:
    311: !      suffix: 2
    312: !      args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
    313: !      requires: !single
    314: !      filter: sed -e "s/[+-]0\.0*i//g"
    315: !
    316: !TEST*/