Actual source code: ex22f.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> ./ex22f [-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: #include <slepc/finclude/slepcnep.h>
36: program ex22f
37: use slepcnep
38: implicit none
40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: ! Declarations
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Mat :: Id, A, B, mats(3) ! problem matrices
45: FN :: f1, f2, f3, funs(3) ! functions to define the nonlinear operator
46: NEP :: nep ! nonlinear eigensolver context
47: NEPType :: tname
48: PetscScalar :: bb, coeffs(2), scal
49: PetscReal :: tau, h, aa, xi, tol
50: PetscInt :: n, i, nev, Istart, Iend
51: PetscMPIInt :: rank
52: PetscErrorCode :: ierr
53: PetscBool :: flg, terse
54: PetscScalar, parameter :: one = 1.0
56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: ! Beginning of program
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
61: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
62: n = 128
63: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
64: tau = 0.001
65: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-tau', tau, flg, ierr))
66: if (rank == 0) then
67: write (*, '(/a,i4,a,f6.3)') 'Delay Eigenproblem, n =', n, ', tau =', tau
68: end if
70: aa = 20.0
71: h = PETSC_PI/real(n + 1, PETSC_REAL_KIND)
73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: ! Create problem matrices
75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: ! ** Id is the identity matrix
78: PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, one, Id, ierr))
79: PetscCallA(MatSetOption(Id, MAT_HERMITIAN, PETSC_TRUE, ierr))
81: ! ** A = 1/h^2*tridiag(1,-2,1) + aa*I
82: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
83: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
84: PetscCallA(MatSetFromOptions(A, ierr))
85: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
86: coeffs(1) = 1.0/(h*h)
87: coeffs(2) = -2.0/(h*h) + aa
88: do i = Istart, Iend - 1
89: if (i > 0) then
90: PetscCallA(MatSetValue(A, i, i - 1, coeffs(1), INSERT_VALUES, ierr))
91: end if
92: if (i < n - 1) then
93: PetscCallA(MatSetValue(A, i, i + 1, coeffs(1), INSERT_VALUES, ierr))
94: end if
95: PetscCallA(MatSetValue(A, i, i, coeffs(2), INSERT_VALUES, ierr))
96: end do
97: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
98: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
99: PetscCallA(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE, ierr))
101: ! ** B = diag(bb(xi))
102: PetscCallA(MatCreate(PETSC_COMM_WORLD, B, ierr))
103: PetscCallA(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
104: PetscCallA(MatSetFromOptions(B, ierr))
105: PetscCallA(MatGetOwnershipRange(B, Istart, Iend, ierr))
106: do i = Istart, Iend - 1
107: xi = (i + 1)*h
108: bb = -4.1 + xi*(1.0 - exp(xi - PETSC_PI))
109: PetscCallA(MatSetValue(B, i, i, bb, INSERT_VALUES, ierr))
110: end do
111: PetscCallA(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
112: PetscCallA(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
113: PetscCallA(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE, ierr))
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: PetscCallA(FNCreate(PETSC_COMM_WORLD, f1, ierr))
120: PetscCallA(FNSetType(f1, FNRATIONAL, ierr))
121: coeffs(1) = -1.0
122: coeffs(2) = 0.0
123: PetscCallA(FNRationalSetNumerator(f1, 2_PETSC_INT_KIND, coeffs, ierr))
125: PetscCallA(FNCreate(PETSC_COMM_WORLD, f2, ierr))
126: PetscCallA(FNSetType(f2, FNRATIONAL, ierr))
127: coeffs(1) = 1.0
128: PetscCallA(FNRationalSetNumerator(f2, 1_PETSC_INT_KIND, coeffs, ierr))
130: PetscCallA(FNCreate(PETSC_COMM_WORLD, f3, ierr))
131: PetscCallA(FNSetType(f3, FNEXP, ierr))
132: scal = -tau
133: PetscCallA(FNSetScale(f3, scal, one, ierr))
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! Create the eigensolver and set various options
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: ! ** Create eigensolver context
140: PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))
142: ! ** Set the split operator. Note that A is passed first so that
143: ! ** SUBSET_NONZERO_PATTERN can be used
144: mats(1) = A
145: mats(2) = Id
146: mats(3) = B
147: funs(1) = f2
148: funs(2) = f1
149: funs(3) = f3
150: PetscCallA(NEPSetSplitOperator(nep, 3_PETSC_INT_KIND, mats, funs, SUBSET_NONZERO_PATTERN, ierr))
151: PetscCallA(NEPSetProblemType(nep, NEP_GENERAL, ierr))
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Customize nonlinear solver; set runtime options
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: tol = 1e-9
158: PetscCallA(NEPSetTolerances(nep, tol, PETSC_CURRENT_INTEGER, ierr))
159: PetscCallA(NEPSetDimensions(nep, 1_PETSC_INT_KIND, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
160: PetscCallA(NEPRIISetLagPreconditioner(nep, 0_PETSC_INT_KIND, ierr))
162: ! ** Set solver parameters at runtime
163: PetscCallA(NEPSetFromOptions(nep, ierr))
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: ! Solve the eigensystem
167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: PetscCallA(NEPSolve(nep, ierr))
171: ! ** Optional: Get some information from the solver and display it
172: PetscCallA(NEPGetType(nep, tname, ierr))
173: if (rank == 0) then
174: write (*, '(a,a)') ' Solution method: ', tname
175: end if
176: PetscCallA(NEPGetDimensions(nep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
177: if (rank == 0) then
178: write (*, '(a,i4)') ' Number of requested eigenvalues:', nev
179: end if
181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: ! Display solution and clean up
183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: ! ** show detailed info unless -terse option is given by user
186: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
187: if (terse) then
188: PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
189: else
190: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
191: PetscCallA(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_WORLD, ierr))
192: PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
193: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
194: end if
195: PetscCallA(NEPDestroy(nep, ierr))
196: PetscCallA(MatDestroy(Id, ierr))
197: PetscCallA(MatDestroy(A, ierr))
198: PetscCallA(MatDestroy(B, ierr))
199: PetscCallA(FNDestroy(f1, ierr))
200: PetscCallA(FNDestroy(f2, ierr))
201: PetscCallA(FNDestroy(f3, ierr))
202: PetscCallA(SlepcFinalize(ierr))
203: end program ex22f
205: !/*TEST
206: !
207: ! test:
208: ! suffix: 1
209: ! args: -terse
210: ! requires: !single
211: ! filter: sed -e "s/[+-]0\.0*i//g"
212: !
213: !TEST*/