Actual source code: test17f.F90
slepc-3.22.1 2024-10-28
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: Test Fortran interface of spectrum-slicing Krylov-Schur.
11: !
12: ! ----------------------------------------------------------------------
13: !
14: program main
15: #include <slepc/finclude/slepceps.h>
16: use slepceps
17: implicit none
19: #define MAXSUB 16
20: #define MAXSHI 16
22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
23: ! Declarations
24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
25: Mat A,B,As,Bs,Au
26: EPS eps
27: ST st
28: KSP ksp
29: PC pc
30: Vec v
31: PetscScalar value
32: PetscInt n,m,i,j,k,Istart,Iend
33: PetscInt nev,ncv,mpd,nval
34: PetscInt row,col,nloc,nlocs,mlocs
35: PetscInt II,npart,inertias(MAXSHI)
36: PetscBool flg,lock
37: PetscMPIInt nprc,rank
38: PetscReal int0,int1,keep,subint(MAXSUB)
39: PetscReal shifts(MAXSHI)
40: PetscScalar eval,one,mone,zero
41: PetscErrorCode ierr
42: MPI_Comm comm
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Beginning of program
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
49: if (ierr .ne. 0) then
50: print*,'SlepcInitialize failed'
51: stop
52: endif
53: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,nprc,ierr))
54: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
55: n = 35
56: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
57: m = n*n
58: if (rank .eq. 0) then
59: write(*,100) n
60: endif
61: 100 format (/'Spectrum-slicing test, n =',I3,' (Fortran)'/)
63: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
64: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr))
65: PetscCallA(MatSetFromOptions(A,ierr))
66: PetscCallA(MatCreate(PETSC_COMM_WORLD,B,ierr))
67: PetscCallA(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr))
68: PetscCallA(MatSetFromOptions(B,ierr))
69: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
70: do II=Istart,Iend-1
71: i = II/n
72: j = II-i*n
73: value = -1.0
74: row = II
75: if (i>0) then
76: col = II-n
77: PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
78: endif
79: if (i<n-1) then
80: col = II+n
81: PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
82: endif
83: if (j>0) then
84: col = II-1
85: PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
86: endif
87: if (j<n-1) then
88: col = II+1
89: PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
90: endif
91: col = II
92: value = 4.0
93: PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
94: value = 2.0
95: PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
96: enddo
97: if (Istart .eq. 0) then
98: row = 0
99: col = 0
100: value = 6.0
101: PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
102: row = 0
103: col = 1
104: value = -1.0
105: PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
106: row = 1
107: col = 0
108: value = -1.0
109: PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
110: row = 1
111: col = 1
112: value = 1.0
113: PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
114: endif
115: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
116: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
117: PetscCallA(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
118: PetscCallA(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Create eigensolver and set various options
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: PetscCallA(EPSCreate(PETSC_COMM_WORLD,eps,ierr))
125: PetscCallA(EPSSetOperators(eps,A,B,ierr))
126: PetscCallA(EPSSetProblemType(eps,EPS_GHEP,ierr))
127: PetscCallA(EPSSetType(eps,EPSKRYLOVSCHUR,ierr))
129: ! Set interval and other settings for spectrum slicing
131: PetscCallA(EPSSetWhichEigenpairs(eps,EPS_ALL,ierr))
132: int0 = 1.1
133: int1 = 1.3
134: PetscCallA(EPSSetInterval(eps,int0,int1,ierr))
135: PetscCallA(EPSGetST(eps,st,ierr))
136: PetscCallA(STSetType(st,STSINVERT,ierr))
137: if (nprc>0) then
138: npart = nprc
139: PetscCallA(EPSKrylovSchurSetPartitions(eps,npart,ierr))
140: endif
141: PetscCallA(EPSKrylovSchurGetKSP(eps,ksp,ierr))
142: PetscCallA(KSPGetPC(ksp,pc,ierr))
143: PetscCallA(KSPSetType(ksp,KSPPREONLY,ierr))
144: PetscCallA(PCSetType(pc,PCCHOLESKY,ierr))
146: ! Test interface functions of Krylov-Schur solver
148: PetscCallA(EPSKrylovSchurGetRestart(eps,keep,ierr))
149: if (rank .eq. 0) then
150: write(*,110) keep
151: endif
152: 110 format (' Restart parameter before changing = ',f7.4)
153: keep = 0.4
154: PetscCallA(EPSKrylovSchurSetRestart(eps,keep,ierr))
155: PetscCallA(EPSKrylovSchurGetRestart(eps,keep,ierr))
156: if (rank .eq. 0) then
157: write(*,120) keep
158: endif
159: 120 format (' ... changed to ',f7.4)
161: PetscCallA(EPSKrylovSchurGetLocking(eps,lock,ierr))
162: if (rank .eq. 0) then
163: write(*,130) lock
164: endif
165: 130 format (' Locking flag before changing = ',L4)
166: PetscCallA(EPSKrylovSchurSetLocking(eps,PETSC_FALSE,ierr))
167: PetscCallA(EPSKrylovSchurGetLocking(eps,lock,ierr))
168: if (rank .eq. 0) then
169: write(*,140) lock
170: endif
171: 140 format (' ... changed to ',L4)
173: PetscCallA(EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr))
174: if (rank .eq. 0) then
175: write(*,150) nev,ncv,mpd
176: endif
177: 150 format (' Sub-solve dimensions before changing: nev=',I2,', ncv=',I2,', mpd=',I2)
178: nev = 30
179: ncv = 60
180: mpd = 60
181: PetscCallA(EPSKrylovSchurSetDimensions(eps,nev,ncv,mpd,ierr))
182: PetscCallA(EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr))
183: if (rank .eq. 0) then
184: write(*,160) nev,ncv,mpd
185: endif
186: 160 format (' ... changed to: nev=',I2,', ncv=',I2,', mpd=',I2)
188: if (nprc>0) then
189: PetscCallA(EPSKrylovSchurGetPartitions(eps,npart,ierr))
190: if (rank .eq. 0) then
191: write(*,170) npart
192: endif
193: 170 format (' Using ',I2,' partitions')
194: if (npart>MAXSUB) then; SETERRA(PETSC_COMM_SELF,1,'Too many subintervals'); endif
196: subint(1) = int0
197: subint(npart+1) = int1
198: do i=2,npart
199: subint(i) = int0+(i-1)*(int1-int0)/npart
200: enddo
201: PetscCallA(EPSKrylovSchurSetSubintervals(eps,subint,ierr))
202: PetscCallA(EPSKrylovSchurGetSubintervals(eps,subint,ierr))
203: if (rank .eq. 0) then
204: write(*,*) 'Using sub-interval separations ='
205: do i=2,npart
206: write(*,180) subint(i)
207: enddo
208: endif
209: 180 format (f7.4)
210: endif
212: PetscCallA(EPSSetFromOptions(eps,ierr))
214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: ! Compute all eigenvalues in interval and display info
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: PetscCallA(EPSSetUp(eps,ierr))
219: PetscCallA(EPSKrylovSchurGetInertias(eps,k,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr))
220: if (k>MAXSHI) then; SETERRA(PETSC_COMM_SELF,1,'Too many shifts'); endif
221: PetscCallA(EPSKrylovSchurGetInertias(eps,k,shifts,inertias,ierr))
222: if (rank .eq. 0) then
223: write(*,*) 'Inertias after EPSSetUp:'
224: do i=1,k
225: write(*,185) shifts(i),inertias(i)
226: enddo
227: endif
228: 185 format (' .. ',f4.1,' (',I3,')')
230: PetscCallA(EPSSolve(eps,ierr))
231: PetscCallA(EPSGetDimensions(eps,nev,ncv,mpd,ierr))
232: PetscCallA(EPSGetInterval(eps,int0,int1,ierr))
233: if (rank .eq. 0) then
234: write(*,190) nev,int0,int1
235: endif
236: 190 format (' Found ',I2,' eigenvalues in interval [',f7.4,',',f7.4,']')
238: if (nprc>0) then
239: PetscCallA(EPSKrylovSchurGetSubcommInfo(eps,k,nval,v,ierr))
240: if (rank .eq. 0) then
241: write(*,200) rank,k,nval
242: do i=0,nval-1
243: PetscCallA(EPSKrylovSchurGetSubcommPairs(eps,i,eval,v,ierr))
244: write(*,210) PetscRealPart(eval)
245: enddo
246: endif
247: 200 format (' Process ',I2,' has worked in sub-interval ',I2,', containing ',I2,' eigenvalues')
248: 210 format (f7.4)
249: PetscCallA(VecDestroy(v,ierr))
251: PetscCallA(EPSKrylovSchurGetSubcommMats(eps,As,Bs,ierr))
252: PetscCallA(MatGetLocalSize(A,nloc,PETSC_NULL_INTEGER,ierr))
253: PetscCallA(MatGetLocalSize(As,nlocs,mlocs,ierr))
254: if (rank .eq. 0) then
255: write(*,220) rank,nloc,nlocs
256: endif
257: 220 format (' Process ',I2,' owns ',I5,', rows of the global',' matrices, and ',I5,' rows in the subcommunicator')
259: ! modify A on subcommunicators
260: PetscCallA(PetscObjectGetComm(As,comm,ierr))
261: PetscCallA(MatCreate(comm,Au,ierr))
262: PetscCallA(MatSetSizes(Au,nlocs,mlocs,m,m,ierr))
263: PetscCallA(MatSetFromOptions(Au,ierr))
264: PetscCallA(MatGetOwnershipRange(Au,Istart,Iend,ierr))
265: do II=Istart,Iend-1
266: value = 0.5
267: PetscCallA(MatSetValue(Au,II,II,value,INSERT_VALUES,ierr))
268: end do
269: PetscCallA(MatAssemblyBegin(Au,MAT_FINAL_ASSEMBLY,ierr))
270: PetscCallA(MatAssemblyEnd(Au,MAT_FINAL_ASSEMBLY,ierr))
271: one = 1.0
272: mone = -1.0
273: zero = 0.0
274: PetscCallA(EPSKrylovSchurUpdateSubcommMats(eps,one,mone,Au,zero,zero,PETSC_NULL_MAT,DIFFERENT_NONZERO_PATTERN,PETSC_TRUE,ierr))
275: PetscCallA(MatDestroy(Au,ierr))
276: endif
278: PetscCallA(EPSDestroy(eps,ierr))
279: PetscCallA(MatDestroy(A,ierr))
280: PetscCallA(MatDestroy(B,ierr))
282: PetscCallA(SlepcFinalize(ierr))
283: end
285: !/*TEST
286: !
287: ! test:
288: ! suffix: 1
289: ! nsize: 2
290: !
291: !TEST*/