Actual source code: test3f.F
slepc-3.11.2 2019-07-30
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2019, 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: SlepcConvMonitor ctx
45: PetscViewerAndFormat vf
47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: ! Beginning of program
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
52: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
53: n = 20
54: if (rank .eq. 0) then
55: write(*,100) n
56: endif
57: 100 format (/'Diagonal Quadratic Eigenproblem, n =',I3,' (Fortran)')
59: call MatCreate(PETSC_COMM_WORLD,A(1),ierr)
60: call MatSetSizes(A(1),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
61: call MatSetFromOptions(A(1),ierr)
62: call MatSetUp(A(1),ierr)
63: call MatGetOwnershipRange(A(1),Istart,Iend,ierr)
64: do i=Istart,Iend-1
65: val = i+1
66: call MatSetValue(A(1),i,i,val,INSERT_VALUES,ierr)
67: enddo
68: call MatAssemblyBegin(A(1),MAT_FINAL_ASSEMBLY,ierr)
69: call MatAssemblyEnd(A(1),MAT_FINAL_ASSEMBLY,ierr)
71: call MatCreate(PETSC_COMM_WORLD,A(2),ierr)
72: call MatSetSizes(A(2),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
73: call MatSetFromOptions(A(2),ierr)
74: call MatSetUp(A(2),ierr)
75: call MatGetOwnershipRange(A(2),Istart,Iend,ierr)
76: do i=Istart,Iend-1
77: val = 1
78: call MatSetValue(A(2),i,i,val,INSERT_VALUES,ierr)
79: enddo
80: call MatAssemblyBegin(A(2),MAT_FINAL_ASSEMBLY,ierr)
81: call MatAssemblyEnd(A(2),MAT_FINAL_ASSEMBLY,ierr)
83: call MatCreate(PETSC_COMM_WORLD,A(3),ierr)
84: call MatSetSizes(A(3),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
85: call MatSetFromOptions(A(3),ierr)
86: call MatSetUp(A(3),ierr)
87: call MatGetOwnershipRange(A(3),Istart,Iend,ierr)
88: do i=Istart,Iend-1
89: val = real(n)/real(i+1)
90: call MatSetValue(A(3),i,i,val,INSERT_VALUES,ierr)
91: enddo
92: call MatAssemblyBegin(A(3),MAT_FINAL_ASSEMBLY,ierr)
93: call MatAssemblyEnd(A(3),MAT_FINAL_ASSEMBLY,ierr)
95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: ! Create eigensolver and test interface functions
97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: call PEPCreate(PETSC_COMM_WORLD,pep,ierr)
100: nmat = 3
101: call PEPSetOperators(pep,nmat,A,ierr)
102: call PEPGetNumMatrices(pep,nmat,ierr)
103: if (rank .eq. 0) then
104: write(*,110) nmat-1
105: endif
106: 110 format (' Polynomial of degree ',I2)
107: i = 0
108: call PEPGetOperators(pep,i,B,ierr)
109: call MatView(B,PETSC_NULL_VIEWER,ierr)
111: call PEPSetType(pep,PEPTOAR,ierr)
112: call PEPGetType(pep,tname,ierr)
113: if (rank .eq. 0) then
114: write(*,120) tname
115: endif
116: 120 format (' Type set to ',A)
118: call PEPGetProblemType(pep,ptype,ierr)
119: if (rank .eq. 0) then
120: write(*,130) ptype
121: endif
122: 130 format (' Problem type before changing = ',I2)
123: call PEPSetProblemType(pep,PEP_HERMITIAN,ierr)
124: call PEPGetProblemType(pep,ptype,ierr)
125: if (rank .eq. 0) then
126: write(*,140) ptype
127: endif
128: 140 format (' ... changed to ',I2)
130: call PEPGetExtract(pep,extr,ierr)
131: if (rank .eq. 0) then
132: write(*,150) extr
133: endif
134: 150 format (' Extraction before changing = ',I2)
135: call PEPSetExtract(pep,PEP_EXTRACT_STRUCTURED,ierr)
136: call PEPGetExtract(pep,extr,ierr)
137: if (rank .eq. 0) then
138: write(*,160) extr
139: endif
140: 160 format (' ... changed to ',I2)
142: alpha = .1
143: its = 5
144: lambda = 1.
145: scal = PEP_SCALE_SCALAR
146: call PEPSetScale(pep,scal,alpha,PETSC_NULL_VEC, &
147: & PETSC_NULL_VEC,its,lambda,ierr)
148: call PEPGetScale(pep,scal,alpha,PETSC_NULL_VEC, &
149: & PETSC_NULL_VEC,its,lambda,ierr)
150: if (rank .eq. 0) then
151: write(*,170) scal,alpha,its
152: endif
153: 170 format (' Scaling: ',I2,', alpha=',F7.4,', its=',I2)
155: basis = PEP_BASIS_CHEBYSHEV1
156: call PEPSetBasis(pep,basis,ierr)
157: call PEPGetBasis(pep,basis,ierr)
158: if (rank .eq. 0) then
159: write(*,180) basis
160: endif
161: 180 format (' Polynomial basis: ',I2)
163: np = 1
164: tol = 1e-9
165: its = 2
166: refine = PEP_REFINE_SIMPLE
167: rscheme = PEP_REFINE_SCHEME_SCHUR
168: call PEPSetRefine(pep,refine,np,tol,its,rscheme,ierr)
169: call PEPGetRefine(pep,refine,np,tol,its,rscheme,ierr)
170: if (rank .eq. 0) then
171: write(*,190) refine,tol,its,rscheme
172: endif
173: 190 format (' Refinement: ',I2,', tol=',F7.4,', its=',I2', schem=',I2)
175: tget = 4.8
176: call PEPSetTarget(pep,tget,ierr)
177: call PEPGetTarget(pep,tget,ierr)
178: call PEPSetWhichEigenpairs(pep,PEP_TARGET_MAGNITUDE,ierr)
179: call PEPGetWhichEigenpairs(pep,which,ierr)
180: if (rank .eq. 0) then
181: write(*,200) which,PetscRealPart(tget)
182: endif
183: 200 format (' Which = ',I2,', target = ',F4.1)
185: nev = 4
186: call PEPSetDimensions(pep,nev,PETSC_DEFAULT_INTEGER, &
187: & PETSC_DEFAULT_INTEGER,ierr)
188: call PEPGetDimensions(pep,nev,ncv,mpd,ierr)
189: if (rank .eq. 0) then
190: write(*,210) nev,ncv,mpd
191: endif
192: 210 format (' Dimensions: nev=',I2,', ncv=',I2,', mpd=',I2)
194: tol = 2.2e-4
195: its = 200
196: call PEPSetTolerances(pep,tol,its,ierr)
197: call PEPGetTolerances(pep,tol,its,ierr)
198: if (rank .eq. 0) then
199: write(*,220) tol,its
200: endif
201: 220 format (' Tolerance =',F8.5,', max_its =',I4)
203: call PEPSetConvergenceTest(pep,PEP_CONV_ABS,ierr)
204: call PEPGetConvergenceTest(pep,conv,ierr)
205: call PEPSetStoppingTest(pep,PEP_STOP_BASIC,ierr)
206: call PEPGetStoppingTest(pep,stp,ierr)
207: if (rank .eq. 0) then
208: write(*,230) conv,stp
209: endif
210: 230 format (' Convergence test =',I2,', stopping test =',I2)
212: call PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, &
213: & PETSC_VIEWER_DEFAULT,vf,ierr)
214: call PEPMonitorSet(pep,PEPMONITORFIRST,vf, &
215: & PetscViewerAndFormatDestroy,ierr)
216: call SlepcConvMonitorCreate(PETSC_VIEWER_STDOUT_WORLD, &
217: & PETSC_VIEWER_DEFAULT,ctx,ierr)
218: call PEPMonitorSet(pep,PEPMONITORCONVERGED,ctx, &
219: & SlepcConvMonitorDestroy,ierr)
220: call PEPMonitorCancel(pep,ierr)
222: call PEPGetST(pep,st,ierr)
223: call STGetKSP(st,ksp,ierr)
224: tol = 1.e-8
225: tolabs = 1.e-35
226: call KSPSetTolerances(ksp,tol,tolabs,PETSC_DEFAULT_REAL, &
227: & PETSC_DEFAULT_INTEGER,ierr)
228: call STView(st,PETSC_NULL_VIEWER,ierr)
229: call PEPGetDS(pep,ds,ierr)
230: call DSView(ds,PETSC_NULL_VIEWER,ierr)
232: call PEPSetFromOptions(pep,ierr)
233: call PEPSolve(pep,ierr)
234: call PEPGetConvergedReason(pep,reason,ierr)
235: if (rank .eq. 0) then
236: write(*,240) reason
237: endif
238: 240 format (' Finished - converged reason =',I2)
240: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: ! Display solution and clean up
242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: call PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr)
244: call PEPDestroy(pep,ierr)
245: call MatDestroy(A(1),ierr)
246: call MatDestroy(A(2),ierr)
247: call MatDestroy(A(3),ierr)
249: call SlepcFinalize(ierr)
250: end
252: !/*TEST
253: !
254: ! test:
255: ! suffix: 1
256: ! args: -pep_tol 1e-6 -pep_ncv 22
257: !
258: !TEST*/