Actual source code: test17f.F90

slepc-3.16.0 2021-09-30
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2021, 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:       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,nprc,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:       if (nprc>0) then
140:         npart = nprc
141:         call EPSKrylovSchurSetPartitions(eps,npart,ierr);CHKERRA(ierr)
142:       endif
143:       call EPSKrylovSchurGetKSP(eps,ksp,ierr);CHKERRA(ierr)
144:       call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
145:       call KSPSetType(ksp,KSPPREONLY,ierr);CHKERRA(ierr)
146:       call PCSetType(pc,PCCHOLESKY,ierr);CHKERRA(ierr)

148: !     Test interface functions of Krylov-Schur solver

150:       call EPSKrylovSchurGetRestart(eps,keep,ierr);CHKERRA(ierr)
151:       if (rank .eq. 0) then
152:         write(*,110) keep
153:       endif
154:  110  format (' Restart parameter before changing = ',f7.4)
155:       keep = 0.4
156:       call EPSKrylovSchurSetRestart(eps,keep,ierr);CHKERRA(ierr)
157:       call EPSKrylovSchurGetRestart(eps,keep,ierr);CHKERRA(ierr)
158:       if (rank .eq. 0) then
159:         write(*,120) keep
160:       endif
161:  120  format (' ... changed to ',f7.4)

163:       call EPSKrylovSchurGetLocking(eps,lock,ierr);CHKERRA(ierr)
164:       if (rank .eq. 0) then
165:         write(*,130) lock
166:       endif
167:  130  format (' Locking flag before changing = ',L4)
168:       call EPSKrylovSchurSetLocking(eps,PETSC_FALSE,ierr);CHKERRA(ierr)
169:       call EPSKrylovSchurGetLocking(eps,lock,ierr);CHKERRA(ierr)
170:       if (rank .eq. 0) then
171:         write(*,140) lock
172:       endif
173:  140  format (' ... changed to ',L4)

175:       call EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
176:       if (rank .eq. 0) then
177:         write(*,150) nev,ncv,mpd
178:       endif
179:  150  format (' Sub-solve dimensions before changing: nev=',I2,', ncv=',I2,', mpd=',I2)
180:       nev = 30
181:       ncv = 60
182:       mpd = 60
183:       call EPSKrylovSchurSetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
184:       call EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
185:       if (rank .eq. 0) then
186:         write(*,160) nev,ncv,mpd
187:       endif
188:  160  format (' ... changed to: nev=',I2,', ncv=',I2,', mpd=',I2)

190:       if (nprc>0) then
191:         call EPSKrylovSchurGetPartitions(eps,npart,ierr);CHKERRA(ierr)
192:         if (rank .eq. 0) then
193:           write(*,170) npart
194:         endif
195:  170    format (' Using ',I2,' partitions')
196:         if (npart>MAXSUB) then; SETERRA(PETSC_COMM_SELF,1,'Too many subintervals'); 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 (f7.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; SETERRA(PETSC_COMM_SELF,1,'Too many shifts'); endif
223:       call EPSKrylovSchurGetInertias(eps,k,shifts,inertias,ierr);CHKERRA(ierr)
224:       if (rank .eq. 0) then
225:         write(*,*) 'Inertias after EPSSetUp:'
226:         do i=1,k
227:           write(*,185) shifts(i),inertias(i)
228:         enddo
229:       endif
230:  185  format (' .. ',f4.1,' (',I3,')')

232:       call EPSSolve(eps,ierr);CHKERRA(ierr)
233:       call EPSGetDimensions(eps,nev,ncv,mpd,ierr);CHKERRA(ierr)
234:       call EPSGetInterval(eps,int0,int1,ierr);CHKERRA(ierr)
235:       if (rank .eq. 0) then
236:         write(*,190) nev,int0,int1
237:       endif
238:  190  format (' Found ',I2,' eigenvalues in interval [',f7.4,',',f7.4,']')

240:       if (nprc>0) then
241:         call EPSKrylovSchurGetSubcommInfo(eps,k,nval,v,ierr);CHKERRA(ierr)
242:         if (rank .eq. 0) then
243:           write(*,200) rank,k,nval
244:           do i=0,nval-1
245:             call EPSKrylovSchurGetSubcommPairs(eps,i,eval,v,ierr);CHKERRA(ierr)
246:             write(*,210) PetscRealPart(eval)
247:           enddo
248:         endif
249:  200    format (' Process ',I2,' has worked in sub-interval ',I2,', containing ',I2,' eigenvalues')
250:  210    format (f7.4)
251:         call VecDestroy(v,ierr);CHKERRA(ierr)

253:         call EPSKrylovSchurGetSubcommMats(eps,As,Bs,ierr);CHKERRA(ierr)
254:         call MatGetLocalSize(A,nloc,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
255:         call MatGetLocalSize(As,nlocs,mlocs,ierr);CHKERRA(ierr)
256:         if (rank .eq. 0) then
257:           write(*,220) rank,nloc,nlocs
258:         endif
259:  220    format (' Process ',I2,' owns ',I5,', rows of the global',' matrices, and ',I5,' rows in the subcommunicator')

261: !       modify A on subcommunicators
262:         call PetscObjectGetComm(As,comm,ierr);CHKERRA(ierr)
263:         call MatCreate(comm,Au,ierr);CHKERRA(ierr)
264:         call MatSetSizes(Au,nlocs,mlocs,m,m,ierr);CHKERRA(ierr)
265:         call MatSetFromOptions(Au,ierr);CHKERRA(ierr)
266:         call MatSetUp(Au,ierr);CHKERRA(ierr)
267:         call MatGetOwnershipRange(Au,Istart,Iend,ierr);CHKERRA(ierr)
268:         do II=Istart,Iend-1
269:           value = 0.5
270:           call MatSetValue(Au,II,II,value,INSERT_VALUES,ierr);CHKERRA(ierr)
271:         end do
272:         call MatAssemblyBegin(Au,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
273:         call MatAssemblyEnd(Au,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
274:         one = 1.0
275:         mone = -1.0
276:         zero = 0.0
277:         call EPSKrylovSchurUpdateSubcommMats(eps,one,mone,Au,zero,zero, &
278:      & PETSC_NULL_MAT,DIFFERENT_NONZERO_PATTERN,PETSC_TRUE,ierr);CHKERRA(ierr)
279:         call MatDestroy(Au,ierr);CHKERRA(ierr)
280:       endif

282:       call EPSDestroy(eps,ierr);CHKERRA(ierr)
283:       call MatDestroy(A,ierr);CHKERRA(ierr)
284:       call MatDestroy(B,ierr);CHKERRA(ierr)

286:       call SlepcFinalize(ierr)
287:       end

289: !/*TEST
290: !
291: !   test:
292: !      suffix: 1
293: !      nsize: 2
294: !
295: !TEST*/