Actual source code: test17f.F90
slepc-3.8.3 2018-04-03
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2017, 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 = ',f6.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 ',f6.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
195: SETERRA(PETSC_COMM_SELF,1,'Too many subintervals')
196: endif
198: subint(1) = int0
199: subint(npart+1) = int1
200: do i=2,npart
201: subint(i) = int0+(i-1)*(int1-int0)/npart
202: enddo
203: call EPSKrylovSchurSetSubintervals(eps,subint,ierr);CHKERRA(ierr)
204: call EPSKrylovSchurGetSubintervals(eps,subint,ierr);CHKERRA(ierr)
205: if (rank .eq. 0) then
206: write(*,*) 'Using sub-interval separations ='
207: do i=2,npart
208: write(*,180) subint(i)
209: enddo
210: endif
211: 180 format (f6.4)
212: endif
214: call EPSSetFromOptions(eps,ierr);CHKERRA(ierr)
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: ! Compute all eigenvalues in interval and display info
218: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: call EPSSetUp(eps,ierr);CHKERRA(ierr)
221: call EPSKrylovSchurGetInertias(eps,k,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
222: if (k>MAXSHI) then
223: SETERRA(PETSC_COMM_SELF,1,'Too many shifts')
224: endif
225: call EPSKrylovSchurGetInertias(eps,k,shifts,inertias,ierr);CHKERRA(ierr)
226: if (rank .eq. 0) then
227: write(*,*) 'Inertias after EPSSetUp:'
228: do i=1,k
229: write(*,185) shifts(i),inertias(i)
230: enddo
231: endif
232: 185 format (' .. ',f3.1,' (',I3,')')
234: call EPSSolve(eps,ierr);CHKERRA(ierr)
235: call EPSGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
236: call EPSGetInterval(eps,int0,int1,ierr);CHKERRA(ierr)
237: if (rank .eq. 0) then
238: write(*,190) nev,int0,int1
239: endif
240: 190 format (' Found ',I2,' eigenvalues in interval [',f6.4,',',f6.4,']')
242: if (size>0) then
243: call EPSKrylovSchurGetSubcommInfo(eps,k,nval,v,ierr);CHKERRA(ierr)
244: if (rank .eq. 0) then
245: write(*,200) rank,k,nval
246: do i=0,nval-1
247: call EPSKrylovSchurGetSubcommPairs(eps,i,eval,v,ierr);CHKERRA(ierr)
248: write(*,210) PetscRealPart(eval)
249: enddo
250: endif
251: 200 format (' Process ',I2,' has worked in sub-interval ',I2,', containing ',I2,' eigenvalues')
252: 210 format (f6.4)
253: call VecDestroy(v,ierr);CHKERRA(ierr)
255: call EPSKrylovSchurGetSubcommMats(eps,As,Bs,ierr);CHKERRA(ierr)
256: call MatGetLocalSize(A,nloc,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
257: call MatGetLocalSize(As,nlocs,mlocs,ierr);CHKERRA(ierr)
258: if (rank .eq. 0) then
259: write(*,220) rank,nloc,nlocs
260: endif
261: 220 format (' Process ',I2,' owns ',I5,', rows of the global',' matrices, and ',I5,' rows in the subcommunicator')
264: ! modify A on subcommunicators
265: call PetscObjectGetComm(As,comm,ierr);CHKERRA(ierr)
266: call MatCreate(comm,Au,ierr);CHKERRA(ierr)
267: call MatSetSizes(Au,nlocs,mlocs,m,m,ierr);CHKERRA(ierr)
268: call MatSetFromOptions(Au,ierr);CHKERRA(ierr)
269: call MatSetUp(Au,ierr);CHKERRA(ierr)
270: call MatGetOwnershipRange(Au,Istart,Iend,ierr);CHKERRA(ierr)
271: do II=Istart,Iend-1
272: value = 0.5
273: call MatSetValue(Au,II,II,value,INSERT_VALUES,ierr);CHKERRA(ierr)
274: end do
275: call MatAssemblyBegin(Au,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
276: call MatAssemblyEnd(Au,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
277: one = 1.0
278: mone = -1.0
279: zero = 0.0
280: call EPSKrylovSchurUpdateSubcommMats(eps,one,mone,Au,zero,zero, &
281: & PETSC_NULL_MAT,DIFFERENT_NONZERO_PATTERN,PETSC_TRUE,ierr);CHKERRA(ierr)
282: call MatDestroy(Au,ierr);CHKERRA(ierr)
283: endif
285: call EPSDestroy(eps,ierr);CHKERRA(ierr)
286: call MatDestroy(A,ierr);CHKERRA(ierr)
287: call MatDestroy(B,ierr);CHKERRA(ierr)
289: call SlepcFinalize(ierr)
290: end