Actual source code: test17f.F90
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: #include <slepc/finclude/slepceps.h>
15: program test17f
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 :: val, eval
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: PetscErrorCode :: ierr
41: MPIU_Comm :: comm
42: PetscScalar, parameter :: one = 1.0, mone = -1.0, zero = 0.0
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Beginning of program
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
49: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, nprc, ierr))
50: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
51: n = 35
52: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
53: m = n*n
54: if (rank == 0) then
55: write (*, '(/a,i3,a/)') 'Spectrum-slicing test, n =', n, ' (Fortran)'
56: end if
58: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
59: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m, ierr))
60: PetscCallA(MatSetFromOptions(A, ierr))
61: PetscCallA(MatCreate(PETSC_COMM_WORLD, B, ierr))
62: PetscCallA(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, m, ierr))
63: PetscCallA(MatSetFromOptions(B, ierr))
64: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
65: do II = Istart, Iend - 1
66: i = II/n
67: j = II - i*n
68: val = -1.0
69: row = II
70: if (i > 0) then
71: col = II - n
72: PetscCallA(MatSetValue(A, row, col, val, INSERT_VALUES, ierr))
73: end if
74: if (i < n - 1) then
75: col = II + n
76: PetscCallA(MatSetValue(A, row, col, val, INSERT_VALUES, ierr))
77: end if
78: if (j > 0) then
79: col = II - 1
80: PetscCallA(MatSetValue(A, row, col, val, INSERT_VALUES, ierr))
81: end if
82: if (j < n - 1) then
83: col = II + 1
84: PetscCallA(MatSetValue(A, row, col, val, INSERT_VALUES, ierr))
85: end if
86: col = II
87: val = 4.0
88: PetscCallA(MatSetValue(A, row, col, val, INSERT_VALUES, ierr))
89: val = 2.0
90: PetscCallA(MatSetValue(B, row, col, val, INSERT_VALUES, ierr))
91: end do
92: if (Istart == 0) then
93: row = 0
94: col = 0
95: val = 6.0
96: PetscCallA(MatSetValue(B, row, col, val, INSERT_VALUES, ierr))
97: row = 0
98: col = 1
99: val = -1.0
100: PetscCallA(MatSetValue(B, row, col, val, INSERT_VALUES, ierr))
101: row = 1
102: col = 0
103: val = -1.0
104: PetscCallA(MatSetValue(B, row, col, val, INSERT_VALUES, ierr))
105: row = 1
106: col = 1
107: val = 1.0
108: PetscCallA(MatSetValue(B, row, col, val, INSERT_VALUES, ierr))
109: end if
110: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
111: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
112: PetscCallA(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
113: PetscCallA(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Create eigensolver and set various options
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))
120: PetscCallA(EPSSetOperators(eps, A, B, ierr))
121: PetscCallA(EPSSetProblemType(eps, EPS_GHEP, ierr))
122: PetscCallA(EPSSetType(eps, EPSKRYLOVSCHUR, ierr))
124: ! ** Set interval and other settings for spectrum slicing
126: PetscCallA(EPSSetWhichEigenpairs(eps, EPS_ALL, ierr))
127: int0 = 1.1
128: int1 = 1.3
129: PetscCallA(EPSSetInterval(eps, int0, int1, ierr))
130: PetscCallA(EPSGetST(eps, st, ierr))
131: PetscCallA(STSetType(st, STSINVERT, ierr))
132: if (nprc > 0) then
133: npart = nprc
134: PetscCallA(EPSKrylovSchurSetPartitions(eps, npart, ierr))
135: end if
136: PetscCallA(EPSKrylovSchurGetKSP(eps, ksp, ierr))
137: PetscCallA(KSPGetPC(ksp, pc, ierr))
138: PetscCallA(KSPSetType(ksp, KSPPREONLY, ierr))
139: PetscCallA(PCSetType(pc, PCCHOLESKY, ierr))
141: ! ** Test interface functions of Krylov-Schur solver
143: PetscCallA(EPSKrylovSchurGetRestart(eps, keep, ierr))
144: if (rank == 0) then
145: write (*, '(a,f7.4)') ' Restart parameter before changing = ', keep
146: end if
147: keep = 0.4
148: PetscCallA(EPSKrylovSchurSetRestart(eps, keep, ierr))
149: PetscCallA(EPSKrylovSchurGetRestart(eps, keep, ierr))
150: if (rank == 0) then
151: write (*, '(a,f7.4)') ' ... changed to ', keep
152: end if
154: PetscCallA(EPSKrylovSchurGetLocking(eps, lock, ierr))
155: if (rank == 0) then
156: write (*, '(a,l4)') ' Locking flag before changing = ', lock
157: end if
158: PetscCallA(EPSKrylovSchurSetLocking(eps, PETSC_FALSE, ierr))
159: PetscCallA(EPSKrylovSchurGetLocking(eps, lock, ierr))
160: if (rank == 0) then
161: write (*, '(a,l4)') ' ... changed to ', lock
162: end if
164: PetscCallA(EPSKrylovSchurGetDimensions(eps, nev, ncv, mpd, ierr))
165: if (rank == 0) then
166: write (*, '(a,i2,a,i2,a,i2)') ' Sub-solve dimensions before changing: nev=', nev, ', ncv=', ncv, ', mpd=', mpd
167: end if
168: nev = 30
169: ncv = 60
170: mpd = 60
171: PetscCallA(EPSKrylovSchurSetDimensions(eps, nev, ncv, mpd, ierr))
172: PetscCallA(EPSKrylovSchurGetDimensions(eps, nev, ncv, mpd, ierr))
173: if (rank == 0) then
174: write (*, '(a,i2,a,i2,a,i2)') ' ... changed to: nev=', nev, ', ncv=', ncv, ', mpd=', mpd
175: end if
177: if (nprc > 0) then
178: PetscCallA(EPSKrylovSchurGetPartitions(eps, npart, ierr))
179: if (rank == 0) then
180: write (*, '(a,i2,a)') ' Using ', npart, ' partitions'
181: end if
182: PetscCheckA(npart <= MAXSUB, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, 'Too many subintervals')
183: subint(1) = int0
184: subint(npart + 1) = int1
185: do i = 2, npart
186: subint(i) = int0 + (i - 1)*(int1 - int0)/npart
187: end do
188: PetscCallA(EPSKrylovSchurSetSubintervals(eps, subint, ierr))
189: PetscCallA(EPSKrylovSchurGetSubintervals(eps, subint, ierr))
190: if (rank == 0) then
191: write (*, *) 'Using sub-interval separations ='
192: do i = 2, npart
193: write (*, '(f7.4)') subint(i)
194: end do
195: end if
196: end if
198: PetscCallA(EPSSetFromOptions(eps, ierr))
200: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: ! Compute all eigenvalues in interval and display info
202: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: PetscCallA(EPSSetUp(eps, ierr))
205: PetscCallA(EPSKrylovSchurGetInertias(eps, k, PETSC_NULL_REAL_ARRAY, PETSC_NULL_INTEGER_ARRAY, ierr))
206: PetscCheckA(k <= MAXSHI, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, 'Too many shifts')
207: PetscCallA(EPSKrylovSchurGetInertias(eps, k, shifts, inertias, ierr))
208: if (rank == 0) then
209: write (*, *) 'Inertias after EPSSetUp:'
210: do i = 1, k
211: write (*, '(a,f4.1,a,i3,a)') ' .. ', shifts(i), ' (', inertias(i), ')'
212: end do
213: end if
215: PetscCallA(EPSSolve(eps, ierr))
216: PetscCallA(EPSGetDimensions(eps, nev, ncv, mpd, ierr))
217: PetscCallA(EPSGetInterval(eps, int0, int1, ierr))
218: if (rank == 0) then
219: write (*, '(a,i2,a,f7.4,a,f7.4,a)') ' Found ', nev, ' eigenvalues in interval [', int0, ',', int1, ']'
220: end if
222: if (nprc > 0) then
223: PetscCallA(EPSKrylovSchurGetSubcommInfo(eps, k, nval, v, ierr))
224: if (rank == 0) then
225: write (*, '(a,i2,a,i2,a,i2,a)') ' Process ', rank, ' has worked in sub-interval ', k, ', containing ', nval, ' eigenvalues'
226: do i = 0, nval - 1
227: PetscCallA(EPSKrylovSchurGetSubcommPairs(eps, i, eval, v, ierr))
228: write (*, '(f7.4)') PetscRealPart(eval)
229: end do
230: end if
231: PetscCallA(VecDestroy(v, ierr))
233: PetscCallA(EPSKrylovSchurGetSubcommMats(eps, As, Bs, ierr))
234: PetscCallA(MatGetLocalSize(A, nloc, PETSC_NULL_INTEGER, ierr))
235: PetscCallA(MatGetLocalSize(As, nlocs, mlocs, ierr))
236: if (rank == 0) then
237: write (*, '(a,i2,a,i5,a,i5,a)') ' Process ', rank, ' owns ', nloc, ', rows of the global matrices, and ', nlocs, ' rows in the subcommunicator'
238: end if
240: ! ** Modify A on subcommunicators
241: PetscCallA(PetscObjectGetComm(As, comm, ierr))
242: PetscCallA(MatCreate(comm, Au, ierr))
243: PetscCallA(MatSetSizes(Au, nlocs, mlocs, m, m, ierr))
244: PetscCallA(MatSetFromOptions(Au, ierr))
245: PetscCallA(MatGetOwnershipRange(Au, Istart, Iend, ierr))
246: do II = Istart, Iend - 1
247: val = 0.5
248: PetscCallA(MatSetValue(Au, II, II, val, INSERT_VALUES, ierr))
249: end do
250: PetscCallA(MatAssemblyBegin(Au, MAT_FINAL_ASSEMBLY, ierr))
251: PetscCallA(MatAssemblyEnd(Au, MAT_FINAL_ASSEMBLY, ierr))
252: PetscCallA(EPSKrylovSchurUpdateSubcommMats(eps, one, mone, Au, zero, zero, PETSC_NULL_MAT, DIFFERENT_NONZERO_PATTERN, PETSC_TRUE, ierr))
253: PetscCallA(MatDestroy(Au, ierr))
254: end if
256: PetscCallA(EPSDestroy(eps, ierr))
257: PetscCallA(MatDestroy(A, ierr))
258: PetscCallA(MatDestroy(B, ierr))
260: PetscCallA(SlepcFinalize(ierr))
261: end program test17f
263: !/*TEST
264: !
265: ! test:
266: ! suffix: 1
267: ! nsize: 2
268: !
269: !TEST*/