Actual source code: test3f.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: Simple example to test the PEP Fortran interface.
11: !
12: ! ----------------------------------------------------------------------
13: !
14: #include <slepc/finclude/slepcpep.h>
15: program test3f
16: use slepcpep
17: implicit none
19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: ! Declarations
21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
23: Mat :: A(3), B
24: PEP :: pep
25: ST :: st
26: KSP :: ksp
27: DS :: ds
28: PetscReal :: tol, tolabs, alpha, lambda
29: PetscScalar :: tget, val
30: PetscInt :: n, i, its, Istart, Iend, nev, ncv, mpd, nmat, np
31: PEPWhich :: which
32: PEPConvergedReason :: reason
33: PEPType :: tname
34: PEPExtract :: extr
35: PEPBasis :: basis
36: PEPScale :: scal
37: PEPRefine :: refine
38: PEPRefineScheme :: rscheme
39: PEPConv :: conv
40: PEPStop :: stp
41: PEPProblemType :: ptype
42: PetscMPIInt :: rank
43: PetscErrorCode :: ierr
44: PetscViewerAndFormat :: vf
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: ! Beginning of program
48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
51: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
52: n = 20
53: if (rank == 0) then
54: write (*, '(/a,i3,a)') 'Diagonal Quadratic Eigenproblem, n =', n, ' (Fortran)'
55: end if
57: PetscCallA(MatCreate(PETSC_COMM_WORLD, A(1), ierr))
58: PetscCallA(MatSetSizes(A(1), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
59: PetscCallA(MatSetFromOptions(A(1), ierr))
60: PetscCallA(MatGetOwnershipRange(A(1), Istart, Iend, ierr))
61: do i = Istart, Iend - 1
62: val = i + 1
63: PetscCallA(MatSetValue(A(1), i, i, val, INSERT_VALUES, ierr))
64: end do
65: PetscCallA(MatAssemblyBegin(A(1), MAT_FINAL_ASSEMBLY, ierr))
66: PetscCallA(MatAssemblyEnd(A(1), MAT_FINAL_ASSEMBLY, ierr))
68: PetscCallA(MatCreate(PETSC_COMM_WORLD, A(2), ierr))
69: PetscCallA(MatSetSizes(A(2), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
70: PetscCallA(MatSetFromOptions(A(2), ierr))
71: PetscCallA(MatGetOwnershipRange(A(2), Istart, Iend, ierr))
72: do i = Istart, Iend - 1
73: val = 1
74: PetscCallA(MatSetValue(A(2), i, i, val, INSERT_VALUES, ierr))
75: end do
76: PetscCallA(MatAssemblyBegin(A(2), MAT_FINAL_ASSEMBLY, ierr))
77: PetscCallA(MatAssemblyEnd(A(2), MAT_FINAL_ASSEMBLY, ierr))
79: PetscCallA(MatCreate(PETSC_COMM_WORLD, A(3), ierr))
80: PetscCallA(MatSetSizes(A(3), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
81: PetscCallA(MatSetFromOptions(A(3), ierr))
82: PetscCallA(MatGetOwnershipRange(A(3), Istart, Iend, ierr))
83: do i = Istart, Iend - 1
84: val = real(n, PETSC_REAL_KIND)/real(i + 1, PETSC_REAL_KIND)
85: PetscCallA(MatSetValue(A(3), i, i, val, INSERT_VALUES, ierr))
86: end do
87: PetscCallA(MatAssemblyBegin(A(3), MAT_FINAL_ASSEMBLY, ierr))
88: PetscCallA(MatAssemblyEnd(A(3), MAT_FINAL_ASSEMBLY, ierr))
90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Create eigensolver and test interface functions
92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: PetscCallA(PEPCreate(PETSC_COMM_WORLD, pep, ierr))
95: nmat = 3
96: PetscCallA(PEPSetOperators(pep, nmat, A, ierr))
97: PetscCallA(PEPGetNumMatrices(pep, nmat, ierr))
98: if (rank == 0) then
99: write (*, '(a,i2)') ' Polynomial of degree ', nmat - 1
100: end if
101: i = 0
102: PetscCallA(PEPGetOperators(pep, i, B, ierr))
103: PetscCallA(MatView(B, PETSC_NULL_VIEWER, ierr))
105: PetscCallA(PEPSetType(pep, PEPTOAR, ierr))
106: PetscCallA(PEPGetType(pep, tname, ierr))
107: if (rank == 0) then
108: write (*, '(a,a)') ' Type set to ', tname
109: end if
111: PetscCallA(PEPGetProblemType(pep, ptype, ierr))
112: if (rank == 0) then
113: write (*, '(a,i2)') ' Problem type before changing = ', ptype
114: end if
115: PetscCallA(PEPSetProblemType(pep, PEP_HERMITIAN, ierr))
116: PetscCallA(PEPGetProblemType(pep, ptype, ierr))
117: if (rank == 0) then
118: write (*, '(a,i2)') ' ... changed to ', ptype
119: end if
121: PetscCallA(PEPGetExtract(pep, extr, ierr))
122: if (rank == 0) then
123: write (*, '(a,i2)') ' Extraction before changing = ', extr
124: end if
125: PetscCallA(PEPSetExtract(pep, PEP_EXTRACT_STRUCTURED, ierr))
126: PetscCallA(PEPGetExtract(pep, extr, ierr))
127: if (rank == 0) then
128: write (*, '(a,i2)') ' ... changed to ', extr
129: end if
131: alpha = .1
132: its = 5
133: lambda = 1.
134: scal = PEP_SCALE_SCALAR
135: PetscCallA(PEPSetScale(pep, scal, alpha, PETSC_NULL_VEC, PETSC_NULL_VEC, its, lambda, ierr))
136: PetscCallA(PEPGetScale(pep, scal, alpha, PETSC_NULL_VEC, PETSC_NULL_VEC, its, lambda, ierr))
137: if (rank == 0) then
138: write (*, '(a,i2,a,f7.4,a,i2)') ' Scaling: ', scal, ', alpha=', alpha, ', its=', its
139: end if
141: basis = PEP_BASIS_CHEBYSHEV1
142: PetscCallA(PEPSetBasis(pep, basis, ierr))
143: PetscCallA(PEPGetBasis(pep, basis, ierr))
144: if (rank == 0) then
145: write (*, '(a,i2)') ' Polynomial basis: ', basis
146: end if
148: np = 1
149: tol = 1e-9
150: its = 2
151: refine = PEP_REFINE_SIMPLE
152: rscheme = PEP_REFINE_SCHEME_SCHUR
153: PetscCallA(PEPSetRefine(pep, refine, np, tol, its, rscheme, ierr))
154: PetscCallA(PEPGetRefine(pep, refine, np, tol, its, rscheme, ierr))
155: if (rank == 0) then
156: write (*, '(a,i2,a,f7.4,a,i2,a,i2)') ' Refinement: ', refine, ', tol=', tol, ', its=', its, ', schem=', rscheme
157: end if
159: tget = 4.8
160: PetscCallA(PEPSetTarget(pep, tget, ierr))
161: PetscCallA(PEPGetTarget(pep, tget, ierr))
162: PetscCallA(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE, ierr))
163: PetscCallA(PEPGetWhichEigenpairs(pep, which, ierr))
164: if (rank == 0) then
165: write (*, '(a,i2,a,f4.1)') ' Which = ', which, ', target = ', PetscRealPart(tget)
166: end if
168: nev = 4
169: PetscCallA(PEPSetDimensions(pep, nev, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
170: PetscCallA(PEPGetDimensions(pep, nev, ncv, mpd, ierr))
171: if (rank == 0) then
172: write (*, '(a,i2,a,i2,a,i2)') ' Dimensions: nev=', nev, ', ncv=', ncv, ', mpd=', mpd
173: end if
175: tol = 2.2e-4
176: its = 200
177: PetscCallA(PEPSetTolerances(pep, tol, its, ierr))
178: PetscCallA(PEPGetTolerances(pep, tol, its, ierr))
179: if (rank == 0) then
180: write (*, '(a,f8.5,a,i4)') ' Tolerance =', tol, ', max_its =', its
181: end if
183: PetscCallA(PEPSetConvergenceTest(pep, PEP_CONV_ABS, ierr))
184: PetscCallA(PEPGetConvergenceTest(pep, conv, ierr))
185: PetscCallA(PEPSetStoppingTest(pep, PEP_STOP_BASIC, ierr))
186: PetscCallA(PEPGetStoppingTest(pep, stp, ierr))
187: if (rank == 0) then
188: write (*, '(a,i2,a,i2)') ' Convergence test =', conv, ', stopping test =', stp
189: end if
191: PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, vf, ierr))
192: PetscCallA(PEPMonitorSet(pep, PEPMONITORFIRST, vf, PetscViewerAndFormatDestroy, ierr))
193: PetscCallA(PEPMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, PETSC_NULL, vf, ierr))
194: PetscCallA(PEPMonitorSet(pep, PEPMONITORCONVERGED, vf, PetscViewerAndFormatDestroy, ierr))
195: PetscCallA(PEPMonitorCancel(pep, ierr))
197: PetscCallA(PEPGetST(pep, st, ierr))
198: PetscCallA(STGetKSP(st, ksp, ierr))
199: tol = 1.e-8
200: tolabs = 1.e-35
201: PetscCallA(KSPSetTolerances(ksp, tol, tolabs, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
202: PetscCallA(STView(st, PETSC_NULL_VIEWER, ierr))
203: PetscCallA(PEPGetDS(pep, ds, ierr))
204: PetscCallA(DSView(ds, PETSC_NULL_VIEWER, ierr))
206: PetscCallA(PEPSetFromOptions(pep, ierr))
207: PetscCallA(PEPSolve(pep, ierr))
208: PetscCallA(PEPGetConvergedReason(pep, reason, ierr))
209: if (rank == 0) then
210: write (*, '(a,i2)') ' Finished - converged reason =', reason
211: end if
213: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: ! Display solution and clean up
215: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: PetscCallA(PEPErrorView(pep, PEP_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
217: PetscCallA(PEPDestroy(pep, ierr))
218: PetscCallA(MatDestroy(A(1), ierr))
219: PetscCallA(MatDestroy(A(2), ierr))
220: PetscCallA(MatDestroy(A(3), ierr))
222: PetscCallA(SlepcFinalize(ierr))
223: end program test3f
225: !/*TEST
226: !
227: ! test:
228: ! suffix: 1
229: ! args: -pep_tol 1e-5 -pep_ncv 22
230: ! filter: sed -e "s/[+-]0\.0*i//g"
231: !
232: !TEST*/