Actual source code: test17f.F90

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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*/