Actual source code: ex22f90.F90

slepc-3.20.1 2023-11-27
Report Typos and Errors
  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> ./ex22f90 [-n <n>] [-tau <tau>] [SLEPc opts]
 11: !
 12: !  Description: Delay differential equation. Fortran90 equivalent of ex22.c
 13: !
 14: !  The command line options are:
 15: !    -n <n>, where <n> = number of grid subdivisions
 16: !    -tau <tau>, where <tau> = delay parameter
 17: !
 18: ! ----------------------------------------------------------------------
 19: !  Solve parabolic partial differential equation with time delay tau
 20: !
 21: !           u_t = u_xx + aa*u(t) + bb*u(t-tau)
 22: !           u(0,t) = u(pi,t) = 0
 23: !
 24: !  with aa = 20 and bb(x) = -4.1+x*(1-exp(x-pi)).
 25: !
 26: !  Discretization leads to a DDE of dimension n
 27: !
 28: !           -u' = A*u(t) + B*u(t-tau)
 29: !
 30: !  which results in the nonlinear eigenproblem
 31: !
 32: !           (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
 33: ! ----------------------------------------------------------------------
 34: !
 35:       program main
 36: #include <slepc/finclude/slepcnep.h>
 37:       use slepcnep
 38:       implicit none

 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41: !     Declarations
 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !
 44: !  Variables:
 45: !     nep       nonlinear eigensolver context
 46: !     Id,A,B    problem matrices
 47: !     f1,f2,f3  functions to define the nonlinear operator

 49:       Mat            Id, A, B, mats(3)
 50:       FN             f1, f2, f3, funs(3)
 51:       NEP            nep
 52:       NEPType        tname
 53:       PetscScalar    one, bb, coeffs(2), scal
 54:       PetscReal      tau, h, aa, xi, tol
 55:       PetscInt       n, i, k, nev, Istart, Iend
 56:       PetscMPIInt    rank
 57:       PetscErrorCode ierr
 58:       PetscBool      flg, terse

 60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61: !     Beginning of program
 62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 64:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
 65:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 66:       n = 128
 67:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
 68:       tau = 0.001
 69:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-tau',tau,flg,ierr))
 70:       if (rank .eq. 0) then
 71:         write(*,100) n, tau
 72:       endif
 73:  100  format (/'Delay Eigenproblem, n =',I4,', tau =',F6.3)

 75:       one = 1.0
 76:       aa = 20.0
 77:       h = PETSC_PI/real(n+1)

 79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80: !     Create problem matrices
 81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 83: !     ** Id is the identity matrix
 84:       PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,one,Id,ierr))
 85:       PetscCallA(MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE,ierr))

 87: !     ** A = 1/h^2*tridiag(1,-2,1) + aa*I
 88:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 89:       PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
 90:       PetscCallA(MatSetFromOptions(A,ierr))
 91:       PetscCallA(MatSetUp(A,ierr))
 92:       PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
 93:       coeffs(1) = 1.0/(h*h)
 94:       coeffs(2) = -2.0/(h*h)+aa
 95:       do i=Istart,Iend-1
 96:         if (i .gt. 0) then
 97:           PetscCallA(MatSetValue(A,i,i-1,coeffs(1),INSERT_VALUES,ierr))
 98:         endif
 99:         if (i .lt. n-1) then
100:           PetscCallA(MatSetValue(A,i,i+1,coeffs(1),INSERT_VALUES,ierr))
101:         endif
102:         PetscCallA(MatSetValue(A,i,i,coeffs(2),INSERT_VALUES,ierr))
103:       end do
104:       PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
105:       PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
106:       PetscCallA(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr))

108: !     ** B = diag(bb(xi))
109:       PetscCallA(MatCreate(PETSC_COMM_WORLD,B,ierr))
110:       PetscCallA(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
111:       PetscCallA(MatSetFromOptions(B,ierr))
112:       PetscCallA(MatSetUp(B,ierr))
113:       PetscCallA(MatGetOwnershipRange(B,Istart,Iend,ierr))
114:       do i=Istart,Iend-1
115:         xi = (i+1)*h
116:         bb  = -4.1+xi*(1.0-exp(xi-PETSC_PI))
117:         PetscCallA(MatSetValue(B,i,i,bb,INSERT_VALUES,ierr))
118:       end do
119:       PetscCallA(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
120:       PetscCallA(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
121:       PetscCallA(MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE,ierr))

123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: !     Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

127:       PetscCallA(FNCreate(PETSC_COMM_WORLD,f1,ierr))
128:       PetscCallA(FNSetType(f1,FNRATIONAL,ierr))
129:       k = 2
130:       coeffs(1) = -1.0
131:       coeffs(2) = 0.0
132:       PetscCallA(FNRationalSetNumerator(f1,k,coeffs,ierr))

134:       PetscCallA(FNCreate(PETSC_COMM_WORLD,f2,ierr))
135:       PetscCallA(FNSetType(f2,FNRATIONAL,ierr))
136:       k = 1
137:       coeffs(1) = 1.0
138:       PetscCallA(FNRationalSetNumerator(f2,k,coeffs,ierr))

140:       PetscCallA(FNCreate(PETSC_COMM_WORLD,f3,ierr))
141:       PetscCallA(FNSetType(f3,FNEXP,ierr))
142:       scal = -tau
143:       PetscCallA(FNSetScale(f3,scal,one,ierr))

145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: !     Create the eigensolver and set various options
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

149: !     ** Create eigensolver context
150:       PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))

152: !     ** Set the split operator. Note that A is passed first so that
153: !     ** SUBSET_NONZERO_PATTERN can be used
154:       k = 3
155:       mats(1) = A
156:       mats(2) = Id
157:       mats(3) = B
158:       funs(1) = f2
159:       funs(2) = f1
160:       funs(3) = f3
161:       PetscCallA(NEPSetSplitOperator(nep,k,mats,funs,SUBSET_NONZERO_PATTERN,ierr))
162:       PetscCallA(NEPSetProblemType(nep,NEP_GENERAL,ierr))

164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: !     Customize nonlinear solver; set runtime options
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

168:       tol = 1e-9
169:       PetscCallA(NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr))
170:       k = 1
171:       PetscCallA(NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr))
172:       k = 0
173:       PetscCallA(NEPRIISetLagPreconditioner(nep,k,ierr))

175: !     ** Set solver parameters at runtime
176:       PetscCallA(NEPSetFromOptions(nep,ierr))

178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: !     Solve the eigensystem
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

182:       PetscCallA(NEPSolve(nep,ierr))

184: !     ** Optional: Get some information from the solver and display it
185:       PetscCallA(NEPGetType(nep,tname,ierr))
186:       if (rank .eq. 0) then
187:         write(*,120) tname
188:       endif
189:  120  format (' Solution method: ',A)
190:       PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
191:       if (rank .eq. 0) then
192:         write(*,130) nev
193:       endif
194:  130  format (' Number of requested eigenvalues:',I4)

196: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: !     Display solution and clean up
198: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

200: !     ** show detailed info unless -terse option is given by user
201:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
202:       if (terse) then
203:         PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
204:       else
205:         PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
206:         PetscCallA(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr))
207:         PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
208:         PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
209:       endif
210:       PetscCallA(NEPDestroy(nep,ierr))
211:       PetscCallA(MatDestroy(Id,ierr))
212:       PetscCallA(MatDestroy(A,ierr))
213:       PetscCallA(MatDestroy(B,ierr))
214:       PetscCallA(FNDestroy(f1,ierr))
215:       PetscCallA(FNDestroy(f2,ierr))
216:       PetscCallA(FNDestroy(f3,ierr))
217:       PetscCallA(SlepcFinalize(ierr))
218:       end

220: !/*TEST
221: !
222: !   test:
223: !      suffix: 1
224: !      args: -terse
225: !      requires: !single
226: !      filter: sed -e "s/[+-]0\.0*i//g"
227: !
228: !TEST*/