Actual source code: test17f.F90
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: 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 size,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: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
49: if (ierr .ne. 0) then
50: print*,'SlepcInitialize failed'
51: stop
52: endif
53: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr);CHKERRA(ierr)
54: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
55: n = 35
56: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(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: call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
64: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);CHKERRA(ierr)
65: call MatSetFromOptions(A,ierr);CHKERRA(ierr)
66: call MatSetUp(A,ierr);CHKERRA(ierr)
67: call MatCreate(PETSC_COMM_WORLD,B,ierr);CHKERRA(ierr)
68: call MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);CHKERRA(ierr)
69: call MatSetFromOptions(B,ierr);CHKERRA(ierr)
70: call MatSetUp(B,ierr);CHKERRA(ierr)
71: call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
72: do II=Istart,Iend-1
73: i = II/n
74: j = II-i*n
75: value = -1.0
76: row = II
77: if (i>0) then
78: col = II-n
79: call MatSetValue(A,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
80: endif
81: if (i<n-1) then
82: col = II+n
83: call MatSetValue(A,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
84: endif
85: if (j>0) then
86: col = II-1
87: call MatSetValue(A,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
88: endif
89: if (j<n-1) then
90: col = II+1
91: call MatSetValue(A,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
92: endif
93: col = II
94: value = 4.0
95: call MatSetValue(A,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
96: value = 2.0
97: call MatSetValue(B,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
98: enddo
99: if (Istart .eq. 0) then
100: row = 0
101: col = 0
102: value = 6.0
103: call MatSetValue(B,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
104: row = 0
105: col = 1
106: value = -1.0
107: call MatSetValue(B,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
108: row = 1
109: col = 0
110: value = -1.0
111: call MatSetValue(B,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
112: row = 1
113: col = 1
114: value = 1.0
115: call MatSetValue(B,row,col,value,INSERT_VALUES,ierr);CHKERRA(ierr)
116: endif
117: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
118: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
119: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
120: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: ! Create eigensolver and set various options
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: call EPSCreate(PETSC_COMM_WORLD,eps,ierr);CHKERRA(ierr)
127: call EPSSetOperators(eps,A,B,ierr);CHKERRA(ierr)
128: call EPSSetProblemType(eps,EPS_GHEP,ierr);CHKERRA(ierr)
129: call EPSSetType(eps,EPSKRYLOVSCHUR,ierr);CHKERRA(ierr)
131: ! Set interval and other settings for spectrum slicing
133: call EPSSetWhichEigenpairs(eps,EPS_ALL,ierr);CHKERRA(ierr)
134: int0 = 1.1
135: int1 = 1.3
136: call EPSSetInterval(eps,int0,int1,ierr);CHKERRA(ierr)
137: call EPSGetST(eps,st,ierr);CHKERRA(ierr)
138: call STSetType(st,STSINVERT,ierr);CHKERRA(ierr)
139: call STGetKSP(st,ksp,ierr);CHKERRA(ierr)
140: call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
141: call KSPSetType(ksp,KSPPREONLY,ierr);CHKERRA(ierr)
142: call PCSetType(pc,PCCHOLESKY,ierr);CHKERRA(ierr)
144: ! Test interface functions of Krylov-Schur solver
146: call EPSKrylovSchurGetRestart(eps,keep,ierr);CHKERRA(ierr)
147: if (rank .eq. 0) then
148: write(*,110) keep
149: endif
150: 110 format (' Restart parameter before changing = ',f7.4)
151: keep = 0.4
152: call EPSKrylovSchurSetRestart(eps,keep,ierr);CHKERRA(ierr)
153: call EPSKrylovSchurGetRestart(eps,keep,ierr);CHKERRA(ierr)
154: if (rank .eq. 0) then
155: write(*,120) keep
156: endif
157: 120 format (' ... changed to ',f7.4)
159: call EPSKrylovSchurGetLocking(eps,lock,ierr);CHKERRA(ierr)
160: if (rank .eq. 0) then
161: write(*,130) lock
162: endif
163: 130 format (' Locking flag before changing = ',L)
164: call EPSKrylovSchurSetLocking(eps,PETSC_FALSE,ierr);CHKERRA(ierr)
165: call EPSKrylovSchurGetLocking(eps,lock,ierr);CHKERRA(ierr)
166: if (rank .eq. 0) then
167: write(*,140) lock
168: endif
169: 140 format (' ... changed to ',L)
171: call EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
172: if (rank .eq. 0) then
173: write(*,150) nev,ncv,mpd
174: endif
175: 150 format (' Sub-solve dimensions before changing: nev=',I2,', ncv=',I2,', mpd=',I2)
176: nev = 30
177: ncv = 60
178: mpd = 60
179: call EPSKrylovSchurSetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
180: call EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
181: if (rank .eq. 0) then
182: write(*,160) nev,ncv,mpd
183: endif
184: 160 format (' ... changed to: nev=',I2,', ncv=',I2,', mpd=',I2)
186: if (size>0) then
187: npart = size
188: call EPSKrylovSchurSetPartitions(eps,npart,ierr);CHKERRA(ierr)
189: call EPSKrylovSchurGetPartitions(eps,npart,ierr);CHKERRA(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: call EPSKrylovSchurSetSubintervals(eps,subint,ierr);CHKERRA(ierr)
202: call EPSKrylovSchurGetSubintervals(eps,subint,ierr);CHKERRA(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: call EPSSetFromOptions(eps,ierr);CHKERRA(ierr)
214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: ! Compute all eigenvalues in interval and display info
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: call EPSSetUp(eps,ierr);CHKERRA(ierr)
219: call EPSKrylovSchurGetInertias(eps,k,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
220: if (k>MAXSHI) then; SETERRA(PETSC_COMM_SELF,1,'Too many shifts'); endif
221: call EPSKrylovSchurGetInertias(eps,k,shifts,inertias,ierr);CHKERRA(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: call EPSSolve(eps,ierr);CHKERRA(ierr)
231: call EPSGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
232: call EPSGetInterval(eps,int0,int1,ierr);CHKERRA(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 (size>0) then
239: call EPSKrylovSchurGetSubcommInfo(eps,k,nval,v,ierr);CHKERRA(ierr)
240: if (rank .eq. 0) then
241: write(*,200) rank,k,nval
242: do i=0,nval-1
243: call EPSKrylovSchurGetSubcommPairs(eps,i,eval,v,ierr);CHKERRA(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: call VecDestroy(v,ierr);CHKERRA(ierr)
251: call EPSKrylovSchurGetSubcommMats(eps,As,Bs,ierr);CHKERRA(ierr)
252: call MatGetLocalSize(A,nloc,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
253: call MatGetLocalSize(As,nlocs,mlocs,ierr);CHKERRA(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')
260: ! modify A on subcommunicators
261: call PetscObjectGetComm(As,comm,ierr);CHKERRA(ierr)
262: call MatCreate(comm,Au,ierr);CHKERRA(ierr)
263: call MatSetSizes(Au,nlocs,mlocs,m,m,ierr);CHKERRA(ierr)
264: call MatSetFromOptions(Au,ierr);CHKERRA(ierr)
265: call MatSetUp(Au,ierr);CHKERRA(ierr)
266: call MatGetOwnershipRange(Au,Istart,Iend,ierr);CHKERRA(ierr)
267: do II=Istart,Iend-1
268: value = 0.5
269: call MatSetValue(Au,II,II,value,INSERT_VALUES,ierr);CHKERRA(ierr)
270: end do
271: call MatAssemblyBegin(Au,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
272: call MatAssemblyEnd(Au,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
273: one = 1.0
274: mone = -1.0
275: zero = 0.0
276: call EPSKrylovSchurUpdateSubcommMats(eps,one,mone,Au,zero,zero, &
277: & PETSC_NULL_MAT,DIFFERENT_NONZERO_PATTERN,PETSC_TRUE,ierr);CHKERRA(ierr)
278: call MatDestroy(Au,ierr);CHKERRA(ierr)
279: endif
281: call EPSDestroy(eps,ierr);CHKERRA(ierr)
282: call MatDestroy(A,ierr);CHKERRA(ierr)
283: call MatDestroy(B,ierr);CHKERRA(ierr)
285: call SlepcFinalize(ierr)
286: end
288: !/*TEST
289: !
290: ! test:
291: ! suffix: 1
292: ! nsize: 2
293: !
294: !TEST*/