Actual source code: test3f.F90
slepc-3.21.1 2024-04-26
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 PEP Fortran interface.
11: !
12: ! ----------------------------------------------------------------------
13: !
14: program main
15: #include <slepc/finclude/slepcpep.h>
16: use slepcpep
17: implicit none
19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: ! Declarations
21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
22: Mat A(3),B
23: PEP pep
24: ST st
25: KSP ksp
26: DS ds
27: PetscReal tol,tolabs,alpha,lambda
28: PetscScalar tget,val
29: PetscInt n,i,its,Istart,Iend
30: PetscInt nev,ncv,mpd,nmat,np
31: PEPWhich which
32: PEPConvergedReason reason
33: PEPType tname
34: PEPExtract extr
35: PEPBasis basis
36: PEPScale scal
37: PEPRefine refine
38: PEPRefineScheme rscheme
39: PEPConv conv
40: PEPStop stp
41: PEPProblemType ptype
42: PetscMPIInt rank
43: PetscErrorCode ierr
44: PetscViewerAndFormat vf
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: ! Beginning of program
48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
51: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
52: n = 20
53: if (rank .eq. 0) then
54: write(*,100) n
55: endif
56: 100 format (/'Diagonal Quadratic Eigenproblem, n =',I3,' (Fortran)')
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: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: ! Create eigensolver and test interface functions
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: PetscCallA(PEPCreate(PETSC_COMM_WORLD,pep,ierr))
96: nmat = 3
97: PetscCallA(PEPSetOperators(pep,nmat,A,ierr))
98: PetscCallA(PEPGetNumMatrices(pep,nmat,ierr))
99: if (rank .eq. 0) then
100: write(*,110) nmat-1
101: endif
102: 110 format (' Polynomial of degree ',I2)
103: i = 0
104: PetscCallA(PEPGetOperators(pep,i,B,ierr))
105: PetscCallA(MatView(B,PETSC_NULL_VIEWER,ierr))
107: PetscCallA(PEPSetType(pep,PEPTOAR,ierr))
108: PetscCallA(PEPGetType(pep,tname,ierr))
109: if (rank .eq. 0) then
110: write(*,120) tname
111: endif
112: 120 format (' Type set to ',A)
114: PetscCallA(PEPGetProblemType(pep,ptype,ierr))
115: if (rank .eq. 0) then
116: write(*,130) ptype
117: endif
118: 130 format (' Problem type before changing = ',I2)
119: PetscCallA(PEPSetProblemType(pep,PEP_HERMITIAN,ierr))
120: PetscCallA(PEPGetProblemType(pep,ptype,ierr))
121: if (rank .eq. 0) then
122: write(*,140) ptype
123: endif
124: 140 format (' ... changed to ',I2)
126: PetscCallA(PEPGetExtract(pep,extr,ierr))
127: if (rank .eq. 0) then
128: write(*,150) extr
129: endif
130: 150 format (' Extraction before changing = ',I2)
131: PetscCallA(PEPSetExtract(pep,PEP_EXTRACT_STRUCTURED,ierr))
132: PetscCallA(PEPGetExtract(pep,extr,ierr))
133: if (rank .eq. 0) then
134: write(*,160) extr
135: endif
136: 160 format (' ... changed to ',I2)
138: alpha = .1
139: its = 5
140: lambda = 1.
141: scal = PEP_SCALE_SCALAR
142: PetscCallA(PEPSetScale(pep,scal,alpha,PETSC_NULL_VEC,PETSC_NULL_VEC,its,lambda,ierr))
143: PetscCallA(PEPGetScale(pep,scal,alpha,PETSC_NULL_VEC,PETSC_NULL_VEC,its,lambda,ierr))
144: if (rank .eq. 0) then
145: write(*,170) scal,alpha,its
146: endif
147: 170 format (' Scaling: ',I2,', alpha=',F7.4,', its=',I2)
149: basis = PEP_BASIS_CHEBYSHEV1
150: PetscCallA(PEPSetBasis(pep,basis,ierr))
151: PetscCallA(PEPGetBasis(pep,basis,ierr))
152: if (rank .eq. 0) then
153: write(*,180) basis
154: endif
155: 180 format (' Polynomial basis: ',I2)
157: np = 1
158: tol = 1e-9
159: its = 2
160: refine = PEP_REFINE_SIMPLE
161: rscheme = PEP_REFINE_SCHEME_SCHUR
162: PetscCallA(PEPSetRefine(pep,refine,np,tol,its,rscheme,ierr))
163: PetscCallA(PEPGetRefine(pep,refine,np,tol,its,rscheme,ierr))
164: if (rank .eq. 0) then
165: write(*,190) refine,tol,its,rscheme
166: endif
167: 190 format (' Refinement: ',I2,', tol=',F7.4,', its=',I2,', schem=',I2)
169: tget = 4.8
170: PetscCallA(PEPSetTarget(pep,tget,ierr))
171: PetscCallA(PEPGetTarget(pep,tget,ierr))
172: PetscCallA(PEPSetWhichEigenpairs(pep,PEP_TARGET_MAGNITUDE,ierr))
173: PetscCallA(PEPGetWhichEigenpairs(pep,which,ierr))
174: if (rank .eq. 0) then
175: write(*,200) which,PetscRealPart(tget)
176: endif
177: 200 format (' Which = ',I2,', target = ',F4.1)
179: nev = 4
180: PetscCallA(PEPSetDimensions(pep,nev,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr))
181: PetscCallA(PEPGetDimensions(pep,nev,ncv,mpd,ierr))
182: if (rank .eq. 0) then
183: write(*,210) nev,ncv,mpd
184: endif
185: 210 format (' Dimensions: nev=',I2,', ncv=',I2,', mpd=',I2)
187: tol = 2.2e-4
188: its = 200
189: PetscCallA(PEPSetTolerances(pep,tol,its,ierr))
190: PetscCallA(PEPGetTolerances(pep,tol,its,ierr))
191: if (rank .eq. 0) then
192: write(*,220) tol,its
193: endif
194: 220 format (' Tolerance =',F8.5,', max_its =',I4)
196: PetscCallA(PEPSetConvergenceTest(pep,PEP_CONV_ABS,ierr))
197: PetscCallA(PEPGetConvergenceTest(pep,conv,ierr))
198: PetscCallA(PEPSetStoppingTest(pep,PEP_STOP_BASIC,ierr))
199: PetscCallA(PEPGetStoppingTest(pep,stp,ierr))
200: if (rank .eq. 0) then
201: write(*,230) conv,stp
202: endif
203: 230 format (' Convergence test =',I2,', stopping test =',I2)
205: PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr))
206: PetscCallA(PEPMonitorSet(pep,PEPMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr))
207: PetscCallA(PEPMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,PETSC_NULL_VEC,vf,ierr))
208: PetscCallA(PEPMonitorSet(pep,PEPMONITORCONVERGED,vf,PEPMonitorConvergedDestroy,ierr))
209: PetscCallA(PEPMonitorCancel(pep,ierr))
211: PetscCallA(PEPGetST(pep,st,ierr))
212: PetscCallA(STGetKSP(st,ksp,ierr))
213: tol = 1.e-8
214: tolabs = 1.e-35
215: PetscCallA(KSPSetTolerances(ksp,tol,tolabs,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr))
216: PetscCallA(STView(st,PETSC_NULL_VIEWER,ierr))
217: PetscCallA(PEPGetDS(pep,ds,ierr))
218: PetscCallA(DSView(ds,PETSC_NULL_VIEWER,ierr))
220: PetscCallA(PEPSetFromOptions(pep,ierr))
221: PetscCallA(PEPSolve(pep,ierr))
222: PetscCallA(PEPGetConvergedReason(pep,reason,ierr))
223: if (rank .eq. 0) then
224: write(*,240) reason
225: endif
226: 240 format (' Finished - converged reason =',I2)
228: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: ! Display solution and clean up
230: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: PetscCallA(PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
232: PetscCallA(PEPDestroy(pep,ierr))
233: PetscCallA(MatDestroy(A(1),ierr))
234: PetscCallA(MatDestroy(A(2),ierr))
235: PetscCallA(MatDestroy(A(3),ierr))
237: PetscCallA(SlepcFinalize(ierr))
238: end
240: !/*TEST
241: !
242: ! test:
243: ! suffix: 1
244: ! args: -pep_tol 1e-6 -pep_ncv 22
245: ! filter: sed -e "s/[+-]0\.0*i//g"
246: !
247: !TEST*/