Actual source code: ex22f90.F90
slepc-3.17.0 2022-03-31
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: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
65: if (ierr .ne. 0) then
66: print*,'SlepcInitialize failed'
67: stop
68: endif
69: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
70: n = 128
71: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(ierr)
72: tau = 0.001
73: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-tau',tau,flg,ierr);CHKERRA(ierr)
74: if (rank .eq. 0) then
75: write(*,100) n, tau
76: endif
77: 100 format (/'Delay Eigenproblem, n =',I4,', tau =',F6.3)
79: one = 1.0
80: aa = 20.0
81: h = PETSC_PI/real(n+1)
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Create problem matrices
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: ! ** Id is the identity matrix
88: call MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,one,Id,ierr);CHKERRA(ierr)
89: call MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)
91: ! ** A = 1/h^2*tridiag(1,-2,1) + aa*I
92: call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
93: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
94: call MatSetFromOptions(A,ierr);CHKERRA(ierr)
95: call MatSetUp(A,ierr);CHKERRA(ierr)
96: call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
97: coeffs(1) = 1.0/(h*h)
98: coeffs(2) = -2.0/(h*h)+aa
99: do i=Istart,Iend-1
100: if (i .gt. 0) then
101: call MatSetValue(A,i,i-1,coeffs(1),INSERT_VALUES,ierr);CHKERRA(ierr)
102: endif
103: if (i .lt. n-1) then
104: call MatSetValue(A,i,i+1,coeffs(1),INSERT_VALUES,ierr);CHKERRA(ierr)
105: endif
106: call MatSetValue(A,i,i,coeffs(2),INSERT_VALUES,ierr);CHKERRA(ierr)
107: end do
108: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
109: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
110: call MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)
112: ! ** B = diag(bb(xi))
113: call MatCreate(PETSC_COMM_WORLD,B,ierr);CHKERRA(ierr)
114: call MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
115: call MatSetFromOptions(B,ierr);CHKERRA(ierr)
116: call MatSetUp(B,ierr);CHKERRA(ierr)
117: call MatGetOwnershipRange(B,Istart,Iend,ierr);CHKERRA(ierr)
118: do i=Istart,Iend-1
119: xi = (i+1)*h
120: bb = -4.1+xi*(1.0-exp(xi-PETSC_PI))
121: call MatSetValue(B,i,i,bb,INSERT_VALUES,ierr);CHKERRA(ierr)
122: end do
123: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
124: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
125: call MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: call FNCreate(PETSC_COMM_WORLD,f1,ierr);CHKERRA(ierr)
132: call FNSetType(f1,FNRATIONAL,ierr);CHKERRA(ierr)
133: k = 2
134: coeffs(1) = -1.0
135: coeffs(2) = 0.0
136: call FNRationalSetNumerator(f1,k,coeffs,ierr);CHKERRA(ierr)
138: call FNCreate(PETSC_COMM_WORLD,f2,ierr);CHKERRA(ierr)
139: call FNSetType(f2,FNRATIONAL,ierr);CHKERRA(ierr)
140: k = 1
141: coeffs(1) = 1.0
142: call FNRationalSetNumerator(f2,k,coeffs,ierr);CHKERRA(ierr)
144: call FNCreate(PETSC_COMM_WORLD,f3,ierr);CHKERRA(ierr)
145: call FNSetType(f3,FNEXP,ierr);CHKERRA(ierr)
146: scal = -tau
147: call FNSetScale(f3,scal,one,ierr);CHKERRA(ierr)
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: ! Create the eigensolver and set various options
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: ! ** Create eigensolver context
154: call NEPCreate(PETSC_COMM_WORLD,nep,ierr);CHKERRA(ierr)
156: ! ** Set the split operator. Note that A is passed first so that
157: ! ** SUBSET_NONZERO_PATTERN can be used
158: k = 3
159: mats(1) = A
160: mats(2) = Id
161: mats(3) = B
162: funs(1) = f2
163: funs(2) = f1
164: funs(3) = f3
165: call NEPSetSplitOperator(nep,k,mats,funs,SUBSET_NONZERO_PATTERN,ierr);CHKERRA(ierr)
166: call NEPSetProblemType(nep,NEP_GENERAL,ierr);CHKERRA(ierr)
168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: ! Customize nonlinear solver; set runtime options
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: tol = 1e-9
173: call NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
174: k = 1
175: call NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
176: k = 0
177: call NEPRIISetLagPreconditioner(nep,k,ierr);CHKERRA(ierr)
179: ! ** Set solver parameters at runtime
180: call NEPSetFromOptions(nep,ierr);CHKERRA(ierr)
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: ! Solve the eigensystem
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: call NEPSolve(nep,ierr);CHKERRA(ierr)
188: ! ** Optional: Get some information from the solver and display it
189: call NEPGetType(nep,tname,ierr);CHKERRA(ierr)
190: if (rank .eq. 0) then
191: write(*,120) tname
192: endif
193: 120 format (' Solution method: ',A)
194: call NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
195: if (rank .eq. 0) then
196: write(*,130) nev
197: endif
198: 130 format (' Number of requested eigenvalues:',I4)
200: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: ! Display solution and clean up
202: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: ! ** show detailed info unless -terse option is given by user
205: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
206: if (terse) then
207: call NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
208: else
209: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
210: call NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
211: call NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
212: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
213: endif
214: call NEPDestroy(nep,ierr);CHKERRA(ierr)
215: call MatDestroy(Id,ierr);CHKERRA(ierr)
216: call MatDestroy(A,ierr);CHKERRA(ierr)
217: call MatDestroy(B,ierr);CHKERRA(ierr)
218: call FNDestroy(f1,ierr);CHKERRA(ierr)
219: call FNDestroy(f2,ierr);CHKERRA(ierr)
220: call FNDestroy(f3,ierr);CHKERRA(ierr)
221: call SlepcFinalize(ierr)
222: end
224: !/*TEST
225: !
226: ! test:
227: ! suffix: 1
228: ! args: -terse
229: ! requires: !single
230: ! filter: sed -e "s/[+-]0\.0*i//g"
231: !
232: !TEST*/