Actual source code: test2f.F90
slepc-3.21.2 2024-09-25
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: ! Description: Simple example to test the NEP Fortran interface.
11: !
12: ! ----------------------------------------------------------------------
13: !
14: program main
15: #include <slepc/finclude/slepcnep.h>
16: use slepcnep
17: implicit none
19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: ! Declarations
21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
22: Mat A(3),B
23: FN f(3),g
24: NEP nep
25: DS ds
26: RG rg
27: PetscReal tol
28: PetscScalar coeffs(2),tget,val
29: PetscInt n,i,its,Istart,Iend
30: PetscInt nev,ncv,mpd,nterm
31: PetscInt nc,np
32: NEPWhich which
33: NEPConvergedReason reason
34: NEPType tname
35: NEPRefine refine
36: NEPRefineScheme rscheme
37: NEPConv conv
38: NEPStop stp
39: NEPProblemType ptype
40: MatStructure mstr
41: PetscMPIInt rank
42: PetscErrorCode ierr
43: PetscViewerAndFormat vf
45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: ! Beginning of program
47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
50: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
51: n = 20
52: if (rank .eq. 0) then
53: write(*,100) n
54: endif
55: 100 format (/'Diagonal Nonlinear Eigenproblem, n =',I3,' (Fortran)')
57: ! Matrices
58: PetscCallA(MatCreate(PETSC_COMM_WORLD,A(1),ierr))
59: PetscCallA(MatSetSizes(A(1),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
60: PetscCallA(MatSetFromOptions(A(1),ierr))
61: PetscCallA(MatGetOwnershipRange(A(1),Istart,Iend,ierr))
62: do i=Istart,Iend-1
63: val = i+1
64: PetscCallA(MatSetValue(A(1),i,i,val,INSERT_VALUES,ierr))
65: enddo
66: PetscCallA(MatAssemblyBegin(A(1),MAT_FINAL_ASSEMBLY,ierr))
67: PetscCallA(MatAssemblyEnd(A(1),MAT_FINAL_ASSEMBLY,ierr))
69: PetscCallA(MatCreate(PETSC_COMM_WORLD,A(2),ierr))
70: PetscCallA(MatSetSizes(A(2),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
71: PetscCallA(MatSetFromOptions(A(2),ierr))
72: PetscCallA(MatGetOwnershipRange(A(2),Istart,Iend,ierr))
73: do i=Istart,Iend-1
74: val = 1
75: PetscCallA(MatSetValue(A(2),i,i,val,INSERT_VALUES,ierr))
76: enddo
77: PetscCallA(MatAssemblyBegin(A(2),MAT_FINAL_ASSEMBLY,ierr))
78: PetscCallA(MatAssemblyEnd(A(2),MAT_FINAL_ASSEMBLY,ierr))
80: PetscCallA(MatCreate(PETSC_COMM_WORLD,A(3),ierr))
81: PetscCallA(MatSetSizes(A(3),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
82: PetscCallA(MatSetFromOptions(A(3),ierr))
83: PetscCallA(MatGetOwnershipRange(A(3),Istart,Iend,ierr))
84: do i=Istart,Iend-1
85: val = real(n)/real(i+1)
86: PetscCallA(MatSetValue(A(3),i,i,val,INSERT_VALUES,ierr))
87: enddo
88: PetscCallA(MatAssemblyBegin(A(3),MAT_FINAL_ASSEMBLY,ierr))
89: PetscCallA(MatAssemblyEnd(A(3),MAT_FINAL_ASSEMBLY,ierr))
91: ! Functions: f0=-lambda, f1=1.0, f2=sqrt(lambda)
92: PetscCallA(FNCreate(PETSC_COMM_WORLD,f(1),ierr))
93: PetscCallA(FNSetType(f(1),FNRATIONAL,ierr))
94: nc = 2
95: coeffs(1) = -1.0
96: coeffs(2) = 0.0
97: PetscCallA(FNRationalSetNumerator(f(1),nc,coeffs,ierr))
99: PetscCallA(FNCreate(PETSC_COMM_WORLD,f(2),ierr))
100: PetscCallA(FNSetType(f(2),FNRATIONAL,ierr))
101: nc = 1
102: coeffs(1) = 1.0
103: PetscCallA(FNRationalSetNumerator(f(2),nc,coeffs,ierr))
105: PetscCallA(FNCreate(PETSC_COMM_WORLD,f(3),ierr))
106: PetscCallA(FNSetType(f(3),FNSQRT,ierr))
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: ! Create eigensolver and test interface functions
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
113: nterm = 3
114: mstr = SAME_NONZERO_PATTERN
115: PetscCallA(NEPSetSplitOperator(nep,nterm,A,f,mstr,ierr))
116: PetscCallA(NEPGetSplitOperatorInfo(nep,nterm,mstr,ierr))
117: if (rank .eq. 0) then
118: write(*,110) nterm
119: endif
120: 110 format (' Nonlinear function with ',I2,' terms')
121: i = 0
122: PetscCallA(NEPGetSplitOperatorTerm(nep,i,B,g,ierr))
123: PetscCallA(MatView(B,PETSC_NULL_VIEWER,ierr))
124: PetscCallA(FNView(g,PETSC_NULL_VIEWER,ierr))
126: PetscCallA(NEPSetType(nep,NEPRII,ierr))
127: PetscCallA(NEPGetType(nep,tname,ierr))
128: if (rank .eq. 0) then
129: write(*,120) tname
130: endif
131: 120 format (' Type set to ',A)
133: PetscCallA(NEPGetProblemType(nep,ptype,ierr))
134: if (rank .eq. 0) then
135: write(*,130) ptype
136: endif
137: 130 format (' Problem type before changing = ',I2)
138: PetscCallA(NEPSetProblemType(nep,NEP_RATIONAL,ierr))
139: PetscCallA(NEPGetProblemType(nep,ptype,ierr))
140: if (rank .eq. 0) then
141: write(*,140) ptype
142: endif
143: 140 format (' ... changed to ',I2)
145: np = 1
146: tol = 1e-9
147: its = 2
148: refine = NEP_REFINE_SIMPLE
149: rscheme = NEP_REFINE_SCHEME_EXPLICIT
150: PetscCallA(NEPSetRefine(nep,refine,np,tol,its,rscheme,ierr))
151: PetscCallA(NEPGetRefine(nep,refine,np,tol,its,rscheme,ierr))
152: if (rank .eq. 0) then
153: write(*,190) refine,tol,its,rscheme
154: endif
155: 190 format (' Refinement: ',I2,', tol=',F12.9,', its=',I2,', scheme=',I2)
157: tget = 1.1
158: PetscCallA(NEPSetTarget(nep,tget,ierr))
159: PetscCallA(NEPGetTarget(nep,tget,ierr))
160: PetscCallA(NEPSetWhichEigenpairs(nep,NEP_TARGET_MAGNITUDE,ierr))
161: PetscCallA(NEPGetWhichEigenpairs(nep,which,ierr))
162: if (rank .eq. 0) then
163: write(*,200) which,PetscRealPart(tget)
164: endif
165: 200 format (' Which = ',I2,', target = ',F4.1)
167: nev = 1
168: ncv = 12
169: PetscCallA(NEPSetDimensions(nep,nev,ncv,PETSC_DEFAULT_INTEGER,ierr))
170: PetscCallA(NEPGetDimensions(nep,nev,ncv,mpd,ierr))
171: if (rank .eq. 0) then
172: write(*,210) nev,ncv,mpd
173: endif
174: 210 format (' Dimensions: nev=',I2,', ncv=',I2,', mpd=',I2)
176: tol = 1.0e-6
177: its = 200
178: PetscCallA(NEPSetTolerances(nep,tol,its,ierr))
179: PetscCallA(NEPGetTolerances(nep,tol,its,ierr))
180: if (rank .eq. 0) then
181: write(*,220) tol,its
182: endif
183: 220 format (' Tolerance =',F9.6,', max_its =',I4)
185: PetscCallA(NEPSetConvergenceTest(nep,NEP_CONV_ABS,ierr))
186: PetscCallA(NEPGetConvergenceTest(nep,conv,ierr))
187: PetscCallA(NEPSetStoppingTest(nep,NEP_STOP_BASIC,ierr))
188: PetscCallA(NEPGetStoppingTest(nep,stp,ierr))
189: if (rank .eq. 0) then
190: write(*,230) conv,stp
191: endif
192: 230 format (' Convergence test =',I2,', stopping test =',I2)
194: PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr))
195: PetscCallA(NEPMonitorSet(nep,NEPMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr))
196: PetscCallA(NEPMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,PETSC_NULL_VEC,vf,ierr))
197: PetscCallA(NEPMonitorSet(nep,NEPMONITORCONVERGED,vf,NEPMonitorConvergedDestroy,ierr))
198: PetscCallA(NEPMonitorCancel(nep,ierr))
200: PetscCallA(NEPGetDS(nep,ds,ierr))
201: PetscCallA(DSView(ds,PETSC_NULL_VIEWER,ierr))
202: PetscCallA(NEPSetFromOptions(nep,ierr))
204: PetscCallA(NEPGetRG(nep,rg,ierr))
205: PetscCallA(RGView(rg,PETSC_NULL_VIEWER,ierr))
207: PetscCallA(NEPSolve(nep,ierr))
208: PetscCallA(NEPGetConvergedReason(nep,reason,ierr))
209: if (rank .eq. 0) then
210: write(*,240) reason
211: endif
212: 240 format (' Finished - converged reason =',I2)
214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: ! Display solution and clean up
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: PetscCallA(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
218: PetscCallA(NEPDestroy(nep,ierr))
219: PetscCallA(MatDestroy(A(1),ierr))
220: PetscCallA(MatDestroy(A(2),ierr))
221: PetscCallA(MatDestroy(A(3),ierr))
222: PetscCallA(FNDestroy(f(1),ierr))
223: PetscCallA(FNDestroy(f(2),ierr))
224: PetscCallA(FNDestroy(f(3),ierr))
226: PetscCallA(SlepcFinalize(ierr))
227: end
229: !/*TEST
230: !
231: ! test:
232: ! suffix: 1
233: ! requires: !single
234: !
235: !TEST*/