Actual source code: ex20f.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> ./ex20f [-n <n>] [SLEPc opts]
     11: !
     12: !  Description: Simple 1-D nonlinear eigenproblem. Fortran90 equivalent of ex20.c
     13: !
     14: !  The command line options are:
     15: !    -n <n>, where <n> = number of grid subdivisions
     16: !
     17: ! ----------------------------------------------------------------------
     18: !  Solve 1-D PDE
     19: !           -u'' = lambda*u
     20: !  on [0,1] subject to
     21: !           u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
     22: ! ----------------------------------------------------------------------
     23: !

     25: #include <slepc/finclude/slepcnep.h>
     26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     27: ! User-defined module with application context and callback functions
     28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     29: module ex20fmodule
     30:   use slepcnep
     31:   type User
     32:     PetscScalar :: kappa
     33:     PetscReal   :: h
     34:   end type User

     36: contains
     37:   ! ---------------  Evaluate Function matrix  T(lambda)  ----------------

     39:   subroutine FormFunction(nep, lambda, fun, B, ctx, ierr)
     40:     implicit none
     41:     NEP            :: nep
     42:     PetscScalar    :: lambda, A(3), c, d
     43:     Mat            :: fun, B
     44:     type(User)     :: ctx
     45:     PetscReal      :: h
     46:     PetscInt       :: i, n, j(3), Istart, Iend
     47:     PetscErrorCode, intent(out) :: ierr

     49: !   ** Compute Function entries and insert into matrix
     50:     PetscCall(MatGetSize(fun, n, PETSC_NULL_INTEGER, ierr))
     51:     PetscCall(MatGetOwnershipRange(fun, Istart, Iend, ierr))
     52:     h = ctx%h
     53:     c = ctx%kappa/(lambda - ctx%kappa)
     54:     d = n

     56: !   ** Boundary points
     57:     if (Istart == 0) then
     58:       i = 0
     59:       j(1) = 0
     60:       j(2) = 1
     61:       A(1) = 2.0*(d - lambda*h/3.0)
     62:       A(2) = -d - lambda*h/6.0
     63:       PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
     64:       Istart = Istart + 1
     65:     end if

     67:     if (Iend == n) then
     68:       i = n - 1
     69:       j(1) = n - 2
     70:       j(2) = n - 1
     71:       A(1) = -d - lambda*h/6.0
     72:       A(2) = d - lambda*h/3.0 + lambda*c
     73:       PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
     74:       Iend = Iend - 1
     75:     end if

     77: !   ** Interior grid points
     78:     do i = Istart, Iend - 1
     79:       j(1) = i - 1
     80:       j(2) = i
     81:       j(3) = i + 1
     82:       A(1) = -d - lambda*h/6.0
     83:       A(2) = 2.0*(d - lambda*h/3.0)
     84:       A(3) = -d - lambda*h/6.0
     85:       PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
     86:     end do

     88: !   ** Assemble matrix
     89:     PetscCall(MatAssemblyBegin(fun, MAT_FINAL_ASSEMBLY, ierr))
     90:     PetscCall(MatAssemblyEnd(fun, MAT_FINAL_ASSEMBLY, ierr))

     92:   end subroutine

     94:   ! ---------------  Evaluate Jacobian matrix  T'(lambda)  ---------------

     96:   subroutine FormJacobian(nep, lambda, jac, ctx, ierr)
     97:     implicit none
     98:     NEP            :: nep
     99:     PetscScalar    :: lambda, A(3), c
    100:     Mat            :: jac
    101:     type(User)     :: ctx
    102:     PetscReal      :: h
    103:     PetscInt       :: i, n, j(3), Istart, Iend
    104:     PetscErrorCode, intent(out) :: ierr

    106: !   ** Compute Jacobian entries and insert into matrix
    107:     PetscCall(MatGetSize(jac, n, PETSC_NULL_INTEGER, ierr))
    108:     PetscCall(MatGetOwnershipRange(jac, Istart, Iend, ierr))
    109:     h = ctx%h
    110:     c = ctx%kappa/(lambda - ctx%kappa)

    112: !   ** Boundary points
    113:     if (Istart == 0) then
    114:       i = 0
    115:       j(1) = 0
    116:       j(2) = 1
    117:       A(1) = -2.0*h/3.0
    118:       A(2) = -h/6.0
    119:       PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
    120:       Istart = Istart + 1
    121:     end if

    123:     if (Iend == n) then
    124:       i = n - 1
    125:       j(1) = n - 2
    126:       j(2) = n - 1
    127:       A(1) = -h/6.0
    128:       A(2) = -h/3.0 - c*c
    129:       PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
    130:       Iend = Iend - 1
    131:     end if

    133: !   ** Interior grid points
    134:     do i = Istart, Iend - 1
    135:       j(1) = i - 1
    136:       j(2) = i
    137:       j(3) = i + 1
    138:       A(1) = -h/6.0
    139:       A(2) = -2.0*h/3.0
    140:       A(3) = -h/6.0
    141:       PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
    142:     end do

    144: !   ** Assemble matrix
    145:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
    146:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))

    148:   end subroutine

    150: end module ex20fmodule

    152: program ex20f
    153:   use ex20fmodule
    154:   implicit none

    156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    157: ! Declarations
    158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    160:   NEP            :: nep      ! nonlinear eigensolver
    161:   Vec            :: x, v     ! eigenvector, auxiliary vector
    162:   PetscScalar    :: lambda   ! eigenvalue
    163:   Mat            :: F, J     ! Function and Jacobian matrices
    164:   type(User)     :: ctx      ! user-defined context
    165:   NEPType        :: tname
    166:   PetscInt       :: n, i, nev, its, maxit, nconv
    167:   PetscReal      :: tol, norm
    168:   PetscMPIInt    :: rank
    169:   PetscBool      :: flg
    170:   PetscErrorCode :: ierr
    171:   PetscScalar, parameter :: one = 1.0

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

    177:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
    178:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    179:   n = 128
    180:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
    181:   if (rank == 0) then
    182:     write (*, '(/a,i4)') 'Nonlinear Eigenproblem, n =', n
    183:   end if

    185:   ctx%h = 1.0/real(n, PETSC_REAL_KIND)
    186:   ctx%kappa = 1.0

    188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    189: ! Create matrix data structure to hold the Function and the Jacobian
    190: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    192:   PetscCallA(MatCreate(PETSC_COMM_WORLD, F, ierr))
    193:   PetscCallA(MatSetSizes(F, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    194:   PetscCallA(MatSetFromOptions(F, ierr))
    195:   PetscCallA(MatSeqAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    196:   PetscCallA(MatMPIAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))

    198:   PetscCallA(MatCreate(PETSC_COMM_WORLD, J, ierr))
    199:   PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    200:   PetscCallA(MatSetFromOptions(J, ierr))
    201:   PetscCallA(MatSeqAIJSetPreallocation(J, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
    202:   PetscCallA(MatMPIAIJSetPreallocation(J, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))

    204: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    205: ! Create the eigensolver and set various options
    206: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    208: ! ** Create eigensolver context
    209:   PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))

    211: ! ** Set routines for evaluation of Function and Jacobian
    212:   PetscCallA(NEPSetFunction(nep, F, F, FormFunction, ctx, ierr))
    213:   PetscCallA(NEPSetJacobian(nep, J, FormJacobian, ctx, ierr))

    215: ! ** Customize nonlinear solver
    216:   tol = 1e-9
    217:   PetscCallA(NEPSetTolerances(nep, tol, PETSC_CURRENT_INTEGER, ierr))
    218:   PetscCallA(NEPSetDimensions(nep, 1_PETSC_INT_KIND, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))

    220: ! ** Set solver parameters at runtime
    221:   PetscCallA(NEPSetFromOptions(nep, ierr))

    223: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    224: ! Solve the eigensystem
    225: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    227: ! ** Evaluate initial guess
    228:   PetscCallA(MatCreateVecs(F, x, PETSC_NULL_VEC, ierr))
    229:   PetscCallA(VecDuplicate(x, v, ierr))
    230:   PetscCallA(VecSet(v, one, ierr))
    231:   PetscCallA(NEPSetInitialSpace(nep, 1_PETSC_INT_KIND, [v], ierr))
    232:   PetscCallA(VecDestroy(v, ierr))

    234: ! ** Call the solver
    235:   PetscCallA(NEPSolve(nep, ierr))
    236:   PetscCallA(NEPGetIterationNumber(nep, its, ierr))
    237:   if (rank == 0) then
    238:     write (*, '(a,i3)') ' Number of NEP iterations =', its
    239:   end if

    241: ! ** Optional: Get some information from the solver and display it
    242:   PetscCallA(NEPGetType(nep, tname, ierr))
    243:   if (rank == 0) then
    244:     write (*, '(a,a10)') ' Solution method: ', tname
    245:   end if
    246:   PetscCallA(NEPGetDimensions(nep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    247:   if (rank == 0) then
    248:     write (*, '(a,i4)') ' Number of requested eigenvalues:', nev
    249:   end if
    250:   PetscCallA(NEPGetTolerances(nep, tol, maxit, ierr))
    251:   if (rank == 0) then
    252:     write (*, '(a,f12.9,a,i5)') ' Stopping condition: tol=', tol, ', maxit=', maxit
    253:   end if

    255: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    256: ! Display solution and clean up
    257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    259:   PetscCallA(NEPGetConverged(nep, nconv, ierr))
    260:   if (rank == 0) then
    261:     write (*, '(a,i2/)') ' Number of converged approximate eigenpairs:', nconv
    262:   end if

    264: ! ** Display eigenvalues and relative errors
    265:   if (nconv > 0) then
    266:     if (rank == 0) then
    267:       write (*, *) '        k              ||T(k)x||'
    268:       write (*, *) '----------------- ------------------'
    269:     end if
    270:     do i = 0, nconv - 1
    271: !     ** Get converged eigenpairs: (in this example they are always real)
    272:       PetscCallA(NEPGetEigenpair(nep, i, lambda, PETSC_NULL_SCALAR, x, PETSC_NULL_VEC, ierr))

    274: !     ** Compute residual norm and error
    275:       PetscCallA(NEPComputeError(nep, i, NEP_ERROR_RELATIVE, norm, ierr))
    276:       if (rank == 0) then
    277:         write (*, '(1p,e15.4,e18.4)') PetscRealPart(lambda), norm
    278:       end if
    279:     end do
    280:     if (rank == 0) then
    281:       write (*, *)
    282:     end if
    283:   end if

    285:   PetscCallA(NEPDestroy(nep, ierr))
    286:   PetscCallA(MatDestroy(F, ierr))
    287:   PetscCallA(MatDestroy(J, ierr))
    288:   PetscCallA(VecDestroy(x, ierr))
    289:   PetscCallA(SlepcFinalize(ierr))
    290: end program ex20f

    292: !/*TEST
    293: !
    294: !   test:
    295: !      suffix: 1
    296: !      args: -nep_target 4
    297: !      filter: sed -e "s/[0-9]\.[0-9]*E-[0-9]*/removed/g" -e "s/ Number of NEP iterations = [ 0-9]*/ Number of NEP iterations = /"
    298: !      requires: !single
    299: !
    300: !TEST*/