GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 52.8% 104 / 0 / 197
Functions: 60.0% 6 / 0 / 10
Branches: 23.5% 117 / 0 / 497
mform/m_mform_eigensolver.f90
Line Branch Exec Source
1 !> @brief Module that provides an interface to SLEPc for solving (generalized) eigenvalue problems
2 module m_mform_eigensolver
3 #include "petsc.fi"
4 use m_common, only: wp
5 use m_mform_basis, only: MFormSpace, size
6 use m_mform_constraint_global, only: MFormConstraint
7 use m_mform_matrix, only: AbstractMatrix
8 use m_mform_solver, only: get_constrained_petsc_matrix, SolverInitOptions, SolverInitTraits
9
10 implicit none
11
12 type, extends(SolverInitOptions) :: EigenInitOptions
13 !> Which part of the spectrum to compute (default: 'largest magnitude')
14 EPSWhich :: which = EPS_LARGEST_MAGNITUDE
15 !> The target value if 'which' is a target option (default: 0.0)
16 real(wp) :: target = 0.0_wp
17 !> The target interval if 'which' is 'all' (default: [0.0, 0.0])
18 real(wp) :: target_interval(2) = [0.0_wp, 0.0_wp]
19 !> The type of eigensolver to use (default: 'krylovschur')
20 EPSType :: eps_type = EPSKRYLOVSCHUR
21 !> Whether 'which' is of target type (default: false)
22 logical :: which_is_target = .false.
23 contains
24 procedure :: init_from_traits => init_eigen_options_from_traits
25 end type EigenInitOptions
26
27 !> @brief Type that wraps SLEPc functionality for solving (generalized) eigenvalue problems
28 type :: EigenSolver
29 !> The SLEPc eigensolver context
30 EPS :: eps
31
32 !> The PETSc matrix (possibly constrained)
33 Mat :: matA, matB
34
35 !> The m-form space of the residual
36 type(MFormSpace) :: space_row
37
38 !> The optional constraint
39 type(MFormConstraint), allocatable :: constraint
40
41 integer :: nr_pairs
42
43 type(EigenInitOptions) :: options
44 contains
45 procedure :: init => EigenSolver_init
46 procedure :: solve => EigenSolver_solve
47 procedure :: get_eigenpair => EigenSolver_get_eigenpair
48 procedure :: get_eigenvalues => EigenSolver_get_eigenvalues
49 procedure :: destroy => EigenSolver_destroy
50 end type EigenSolver
51
52 interface EigenInitOptions
53 module procedure init_EigenInitOptions
54 end interface EigenInitOptions
55 contains
56
57 !> @brief Initialize an EigenInitOptions object
58 !>
59 !> @param ksp_type _(optional)_ The KSP type to use (default: auto-select based on traits)
60 !> @param pc_type _(optional)_ The PC type to use (default: auto-select based on traits)
61 !> @param hypre_type _(optional)_ The HYPRE preconditioner type to use (default: auto-select based on traits)
62 !> @param rtol _(optional)_ The relative tolerance (default: 1.e-10)
63 !> @param atol _(optional)_ The absolute tolerance (default: 1.e-20)
64 !> @param max_iter _(optional)_ The maximum number of iterations (default: size of the space)
65 !> @param verbosity _(optional)_ The verbosity level (default: 'warn on failure')
66 !> - 'warn never'
67 !> - 'warn on failure'
68 !> - 'warn always'
69 !> - 'error on failure'
70 !> @param which _(optional)_ Which part of the spectrum to compute. Possible values are:
71 !> @param which _(optional)_ Which part of the spectrum to compute. Possible values are:
72 !> - 'largest magnitude' or 'lm' (default)
73 !> - 'smallest magnitude' or 'sm'
74 !> - 'largest real' or 'lr'
75 !> - 'smallest real' or 'sr'
76 !> - 'largest imaginary' or 'li'
77 !> - 'smallest imaginary' or 'si'
78 !> - 'target magnitude' or 'tm'
79 !> - 'target real' or 'tr'
80 !> - 'target imaginary' or 'ti'
81 !> - 'target interval' or 'all'
82 !> @param target _(optional)_ The target value if 'which' is a target option
83 !> @param target_interval _(optional)_ The target interval if 'which' is 'all'
84 !> @param eps_type _(optional)_ The type of eigensolver to use (default: 'krylovschur'). Possible values are:
85 !> - 'power'
86 !> - 'subspace'
87 !> - 'arnoldi'
88 !> - 'lanczos'
89 !> - 'krylovschur'
90 !> - 'gd'
91 !> - 'jd' or 'jacobi davidson'
92 !> - 'rqcg'
93 !> - 'lobpcg'
94 !> - 'ciss'
95 !> - 'lyapii'
96 !> - 'lapack'
97 !> - 'arpack'
98 !> - 'blopex'
99 !> - 'primme'
100 !> - 'feast'
101 !> - 'scalapack'
102 !> - 'elpa'
103 !> - 'elemental'
104 !> - 'evsl'
105 !> - 'chase'
106 22 function init_EigenInitOptions(ksp_type, pc_type, hypre_type, rtol, atol, max_iter, verbosity, &
107 which, target, target_interval, eps_type) result(options)
108 type(EigenInitOptions) :: options
109 character(*), optional, intent(in) :: ksp_type, pc_type, hypre_type
110 real(wp), optional, intent(in) :: rtol, atol
111 integer, optional, intent(in) :: max_iter
112 character(len=80), optional, intent(in) :: verbosity
113 character(*), optional, intent(in) :: which
114 real(wp), optional, intent(in) :: target
115 real(wp), optional, intent(in) :: target_interval(2)
116 character(*), optional, intent(in) :: eps_type
117
118 type(SolverInitOptions) :: base_options
119
120 base_options = SolverInitOptions(ksp_type=ksp_type, pc_type=pc_type, hypre_type=hypre_type, rtol=rtol, atol=atol, &
121
3/6
✓ Branch 7 → 8 taken 22 times.
✗ Branch 7 → 9 not taken.
✓ Branch 9 → 10 taken 22 times.
✗ Branch 9 → 11 not taken.
✓ Branch 11 → 12 taken 22 times.
✗ Branch 11 → 13 not taken.
 88 max_iter=max_iter, verbosity=verbosity)
122 22 options%ksp_type = base_options%ksp_type
123 22 options%pc_type = base_options%pc_type
124 22 options%hypre_type = base_options%hypre_type
125 22 options%ilu_level = base_options%ilu_level
126 22 options%rtol = base_options%rtol
127 22 options%atol = base_options%atol
128 22 options%max_iter = base_options%max_iter
129 22 options%verbosity = base_options%verbosity
130 #ifndef SLEPC_IS_LINKED
131 error stop "EigenInitOptions::init: SLEPc is not linked. Cannot use EigenInitOptions."
132 #else
133 options%which_is_target = .false.
134
1/2
✓ Branch 14 → 15 taken 22 times.
✗ Branch 14 → 30 not taken.
 22 if (present(which)) then
135 select case (which)
136 case ('largest magnitude', 'lm')
137 ✗ options%which = EPS_LARGEST_MAGNITUDE
138 case ('smallest magnitude', 'sm')
139 ✗ options%which = EPS_SMALLEST_MAGNITUDE
140 case ('largest real', 'lr')
141 ✗ options%which = EPS_LARGEST_REAL
142 case ('smallest real', 'sr')
143 16 options%which = EPS_SMALLEST_REAL
144 case ('largest imaginary', 'li')
145 ✗ options%which = EPS_LARGEST_IMAGINARY
146 case ('smallest imaginary', 'si')
147 ✗ options%which = EPS_SMALLEST_IMAGINARY
148 case ('target magnitude', 'tm')
149 ✗ options%which = EPS_TARGET_MAGNITUDE
150 ✗ options%which_is_target = .true.
151 case ('target real', 'tr')
152 6 options%which = EPS_TARGET_REAL
153 6 options%which_is_target = .true.
154 case ('target imaginary', 'ti')
155 ✗ options%which = EPS_TARGET_IMAGINARY
156 ✗ options%which_is_target = .true.
157 case ('target interval', 'all')
158 ✗ options%which = EPS_ALL
159 case default
160 ✗ write (*, '(A, A)') "Error: unknown option for 'which': ", which
161
2/11
✗ Branch 15 → 16 not taken.
✗ Branch 15 → 17 not taken.
✓ Branch 15 → 18 taken 16 times.
✗ Branch 15 → 19 not taken.
✗ Branch 15 → 20 not taken.
✗ Branch 15 → 21 not taken.
✓ Branch 15 → 22 taken 6 times.
✗ Branch 15 → 23 not taken.
✗ Branch 15 → 24 not taken.
✗ Branch 15 → 25 not taken.
✗ Branch 15 → 30 not taken.
 22 error stop
162 end select
163 else
164 options%which = EPS_LARGEST_MAGNITUDE
165 end if
166
167
1/2
✓ Branch 30 → 31 taken 22 times.
✗ Branch 30 → 56 not taken.
 22 if (present(eps_type)) then
168 ✗ select case (eps_type)
169 case ('power')
170 ✗ options%eps_type = EPSPOWER
171 case ('subspace')
172 ✗ options%eps_type = EPSSUBSPACE
173 case ('arnoldi')
174 8 options%eps_type = EPSARNOLDI
175 case ('lanczos')
176 ✗ options%eps_type = EPSLANCZOS
177 case ('krylovschur')
178 14 options%eps_type = EPSKRYLOVSCHUR
179 case ('gd')
180 ✗ options%eps_type = EPSGD
181 case ('jd', 'jacobi davidson')
182 ✗ options%eps_type = EPSJD
183 case ('rqcg')
184 ✗ options%eps_type = EPSRQCG
185 case ('lobpcg')
186 ✗ options%eps_type = EPSLOBPCG
187 case ('ciss')
188 ✗ options%eps_type = EPSCISS
189 case ('lyapii')
190 ✗ options%eps_type = EPSLYAPII
191 case ('lapack')
192 ✗ options%eps_type = EPSLAPACK
193 case ('arpack')
194 ✗ options%eps_type = EPSARPACK
195 case ('blopex')
196 ✗ options%eps_type = EPSBLOPEX
197 case ('primme')
198 ✗ options%eps_type = EPSPRIMME
199 case ('feast')
200 ✗ options%eps_type = EPSFEAST
201 case ('scalapack')
202 ✗ options%eps_type = EPSSCALAPACK
203 ! case ('elpa')
204 ! eps_type_ = EPSELPA
205 case ('elemental')
206 ✗ options%eps_type = EPSELEMENTAL
207 case ('evsl')
208 ✗ options%eps_type = EPSEVSL
209 ! case ('chase')
210 ! eps_type_ = EPSCHASE
211 case default
212 ✗ write (*, '(A, A)') "Error: unknown option for 'eps_type': ", eps_type
213
2/20
✗ Branch 31 → 32 not taken.
✗ Branch 31 → 33 not taken.
✓ Branch 31 → 34 taken 8 times.
✗ Branch 31 → 35 not taken.
✓ Branch 31 → 36 taken 14 times.
✗ Branch 31 → 37 not taken.
✗ Branch 31 → 38 not taken.
✗ Branch 31 → 39 not taken.
✗ Branch 31 → 40 not taken.
✗ Branch 31 → 41 not taken.
✗ Branch 31 → 42 not taken.
✗ Branch 31 → 43 not taken.
✗ Branch 31 → 44 not taken.
✗ Branch 31 → 45 not taken.
✗ Branch 31 → 46 not taken.
✗ Branch 31 → 47 not taken.
✗ Branch 31 → 48 not taken.
✗ Branch 31 → 49 not taken.
✗ Branch 31 → 50 not taken.
✗ Branch 31 → 51 not taken.
 22 error stop
214 end select
215 else
216 ✗ options%eps_type = EPSKRYLOVSCHUR
217 end if
218
219
2/2
✓ Branch 57 → 58 taken 6 times.
✓ Branch 57 → 63 taken 16 times.
 22 if (present(target)) then
220 6 options%target = target
221
1/2
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 63 taken 6 times.
 6 if (.not. options%which_is_target) then
222 ✗ write (*, '(A)') "EigenSolver::init Warning: 'target' specified but 'which' is not a target option"
223 end if
224 else
225 options%target = 0.0_wp
226 end if
227
228
1/2
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 71 taken 22 times.
 22 if (present(target_interval)) then
229 ✗ options%target_interval = target_interval
230 ✗ if (options%which /= EPS_ALL) then
231 ✗ write (*, '(A)') "EigenSolver::init Warning: 'target_interval' specified but 'which' is not 'all'"
232 end if
233 else
234
2/2
✓ Branch 71 → 72 taken 44 times.
✓ Branch 71 → 73 taken 22 times.
 66 options%target_interval = [0.0_wp, 0.0_wp]
235 end if
236 #endif
237
3/4
✓ Branch 3 → 4 taken 44 times.
✓ Branch 3 → 5 taken 22 times.
✓ Branch 5 → 6 taken 22 times.
✗ Branch 5 → 7 not taken.
 66 end function init_EigenInitOptions
238
239 !> @brief Select solver options for eigenvalue problems, favouring direct solvers (MUMPS)
240 !>
241 !> Eigenvalue solvers require highly accurate linear solves inside the spectral
242 !> transform. Iterative methods with standard preconditioners often lack the
243 !> accuracy needed for convergence of multiple eigenpairs. This override
244 !> defaults to a direct solver (PREONLY + Cholesky/LU with MUMPS) which works
245 !> both in serial and parallel.
246 !>
247 !> @param this The EigenInitOptions object
248 !> @param traits The solver traits derived from matrix properties
249 22 subroutine init_eigen_options_from_traits(this, traits)
250 class(EigenInitOptions), intent(inout) :: this
251 type(SolverInitTraits), intent(in) :: traits
252
253
1/2
✓ Branch 2 → 3 taken 22 times.
✗ Branch 2 → 4 not taken.
 22 if (len_trim(this%ksp_type) == 0) then
254 22 this%ksp_type = 'preonly'
255 end if
256
257
1/2
✓ Branch 4 → 5 taken 22 times.
✗ Branch 4 → 9 not taken.
 22 if (len_trim(this%pc_type) == 0) then
258
2/4
✓ Branch 5 → 6 taken 22 times.
✗ Branch 5 → 8 not taken.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 22 times.
 22 if (traits%is_symmetric .and. traits%is_spd) then
259 ✗ this%pc_type = 'cholesky'
260 else
261 22 this%pc_type = 'lu'
262 end if
263 end if
264 22 end subroutine init_eigen_options_from_traits
265
266 !> @brief Initialize the eigensolver
267 !>
268 !> @param this The EigenSolver object
269 !> @param amat The matrix A in the eigenvalue problem: $A x = \lambda x$
270 !> @param bmat _(optional)_ The matrix B in the eigenvalue problem: $A x = \lambda B x$
271 !> @param[in] constraint1 _(optional)_ An optional constraint to be applied to the solution
272 !> @param[in] constraint2 _(optional)_ An optional constraint to be applied to the solution
273 !> @param[in] constraint3 _(optional)_ An optional constraint to be applied to the solution
274 !> @param[in] coord_transform _(optional)_ The coordinate transformation associated with the m-form spaces
275 !> @param options _(optional)_ The options for initializing the eigensolver (default: auto-select based on traits)
276 22 subroutine EigenSolver_init(this, amat, bmat, constraint1, constraint2, constraint3, coord_transform, options)
277 use m_mform_constraint_abstract, only: MFormConstraintLocal
278 use m_coord_transform_abstract, only: CoordTransformAbstract
279 implicit none
280
281 class(EigenSolver), intent(inout) :: this
282 class(AbstractMatrix), intent(in) :: amat
283 class(AbstractMatrix), optional, intent(in) :: bmat
284 class(MFormConstraintLocal), optional, intent(in) :: constraint1, constraint2, constraint3
285 class(CoordTransformAbstract), optional, intent(in) :: coord_transform
286 type(EigenInitOptions), optional, intent(in) :: options
287
288 integer :: ierr
289 EPSProblemType :: eps_problem_type
290 type(SolverInitTraits) :: traits
291 ST :: st
292 KSP :: ksp
293 PC :: pc
294
295 #ifndef SLEPC_IS_LINKED
296 error stop "EigenSolver::init: SLEPc is not linked. Cannot use EigenSolver."
297 #else
298 22 call this%destroy()
299
300
1/2
✓ Branch 3 → 4 taken 22 times.
✗ Branch 3 → 5 not taken.
 22 if (present(options)) this%options = options
301
302
5/12
✓ Branch 5 → 6 taken 14 times.
✓ Branch 5 → 7 taken 8 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 11 taken 14 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 8 times.
✗ Branch 8 → 9 not taken.
✗ Branch 8 → 11 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 33 taken 8 times.
✗ Branch 10 → 11 not taken.
✗ Branch 10 → 33 not taken.
 22 if (present(constraint1) .or. present(constraint2) .or. present(constraint3)) then
303
4/6
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 14 times.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 14 times.
✓ Branch 16 → 17 taken 70 times.
✓ Branch 16 → 18 taken 14 times.
 84 allocate (this%constraint)
304 14 call this%constraint%init(amat%space_col)
305
306
2/4
✓ Branch 19 → 20 taken 14 times.
✗ Branch 19 → 23 not taken.
✓ Branch 20 → 21 taken 14 times.
✗ Branch 20 → 23 not taken.
 14 if (present(constraint1)) call this%constraint%insert(constraint1, coord_transform=coord_transform)
307
2/4
✓ Branch 23 → 24 taken 14 times.
✗ Branch 23 → 27 not taken.
✓ Branch 24 → 25 taken 14 times.
✗ Branch 24 → 27 not taken.
 14 if (present(constraint2)) call this%constraint%insert(constraint2, coord_transform=coord_transform)
308
1/4
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 31 taken 14 times.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 31 not taken.
 14 if (present(constraint3)) call this%constraint%insert(constraint3, coord_transform=coord_transform)
309
310 14 call this%constraint%assemble()
311 end if
312
313
2/4
✓ Branch 33 → 34 taken 22 times.
✗ Branch 33 → 44 not taken.
✓ Branch 34 → 35 taken 22 times.
✗ Branch 34 → 44 not taken.
 22 if (present(bmat)) then
314
1/2
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 40 taken 22 times.
 22 if (size(amat%space_row) /= size(bmat%space_row)) then
315 ✗ write (*, '(A)') "Error: EigenSolver::init: the row space of A and B must have the same size"
316 ✗ error stop
317 end if
318
2/4
✓ Branch 40 → 41 taken 22 times.
✗ Branch 40 → 43 not taken.
✓ Branch 41 → 42 taken 22 times.
✗ Branch 41 → 43 not taken.
 22 if (amat%is_symmetric .and. bmat%is_symmetric) then
319 22 eps_problem_type = EPS_GHEP
320 else
321 ✗ eps_problem_type = EPS_GNHEP
322 end if
323 else
324 ✗ if (amat%is_symmetric) then
325 ✗ eps_problem_type = EPS_HEP
326 else
327 ✗ eps_problem_type = EPS_NHEP
328 end if
329 end if
330
331 #endif
332
333
5/8
✓ Branch 47 → 48 taken 22 times.
✗ Branch 47 → 51 not taken.
✓ Branch 48 → 49 taken 22 times.
✗ Branch 48 → 50 not taken.
✓ Branch 51 → 52 taken 22 times.
✗ Branch 51 → 54 not taken.
✓ Branch 52 → 53 taken 19 times.
✓ Branch 52 → 54 taken 3 times.
 22 this%space_row = amat%space_row
334
1/2
✗ Branch 55 → 56 not taken.
✓ Branch 55 → 58 taken 22 times.
 22 SlepcCall(EPSCreate(this%space_row%tp_spaces(1)%domain%comm, this%eps, ierr))
335
336
1/2
✗ Branch 59 → 60 not taken.
✓ Branch 59 → 62 taken 22 times.
 22 SlepcCall(EPSSetType(this%eps, this%options%eps_type, ierr))
337
1/2
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 66 taken 22 times.
 22 SlepcCall(EPSSetProblemType(this%eps, eps_problem_type, ierr))
338
339
2/2
✓ Branch 66 → 67 taken 8 times.
✓ Branch 66 → 70 taken 14 times.
 22 if (.not. allocated(this%constraint)) then
340 ! NOTE a copy is required because the destructor of amat will destroy the matrix
341
1/2
✗ Branch 68 → 69 not taken.
✓ Branch 68 → 73 taken 8 times.
 8 PetscCall(MatConvert(amat%mat, "same", MAT_INITIAL_MATRIX, this%matA, ierr))
342 else
343
344 14 call get_constrained_petsc_matrix(this%matA, amat%mat, this%constraint, coeff=amat%constraint_coeff())
345 end if
346
2/4
✓ Branch 73 → 74 taken 22 times.
✗ Branch 73 → 75 not taken.
✗ Branch 74 → 75 not taken.
✓ Branch 74 → 79 taken 22 times.
 22 if (.not. present(bmat)) then
347 ✗ SlepcCall(EPSSetOperators(this%eps, this%matA, PETSC_NULL_MAT, ierr))
348 else
349
2/2
✓ Branch 79 → 80 taken 8 times.
✓ Branch 79 → 83 taken 14 times.
 22 if (.not. allocated(this%constraint)) then
350 ! NOTE a copy is required because the destructor of amat will destroy the matrix
351
1/2
✗ Branch 81 → 82 not taken.
✓ Branch 81 → 86 taken 8 times.
 8 PetscCall(MatConvert(bmat%mat, "same", MAT_INITIAL_MATRIX, this%matB, ierr))
352 else
353
354 14 call get_constrained_petsc_matrix(this%matB, bmat%mat, this%constraint, coeff=bmat%constraint_coeff())
355 end if
356
357
1/2
✗ Branch 87 → 88 not taken.
✓ Branch 87 → 90 taken 22 times.
 22 SlepcCall(EPSSetOperators(this%eps, this%matA, this%matB, ierr))
358 end if
359
360
1/2
✗ Branch 91 → 92 not taken.
✓ Branch 91 → 94 taken 22 times.
 22 SlepcCall(EPSSetWhichEigenpairs(this%eps, this%options%which, ierr))
361
2/2
✓ Branch 94 → 95 taken 6 times.
✓ Branch 94 → 99 taken 16 times.
 22 if (this%options%which_is_target) then
362
1/2
✗ Branch 96 → 97 not taken.
✓ Branch 96 → 99 taken 6 times.
 6 SlepcCall(EPSSetTarget(this%eps, this%options%target, ierr))
363 end if
364
365
1/2
✗ Branch 99 → 100 not taken.
✓ Branch 99 → 104 taken 22 times.
 22 if (this%options%which == EPS_ALL) then
366 ✗ SlepcCall(EPSSetInterval(this%eps, this%options%target_interval(1), this%options%target_interval(2), ierr))
367 end if
368
369
370 ! Build traits from matrix properties for automatic solver selection
371 22 traits%is_symmetric = amat%is_symmetric
372 ! The shifted operator (A - sigma*B) is generally indefinite, so not SPD
373 traits%is_spd = .false.
374 22 traits%is_parallel = amat%space_row%tp_spaces(1)%domain%nr_ranks > 1
375 traits%is_diffdiff = .false. ! Irrelevant
376 22 traits%m_form = amat%space_col%m
377 22 traits%allow_ams = .false.
378
379 22 call this%options%init_from_traits(traits)
380
381 ! Configure the Spectral Transform (ST) and its linear solver (KSP)
382
1/2
✗ Branch 106 → 107 not taken.
✓ Branch 106 → 109 taken 22 times.
 22 SlepcCall(EPSGetST(this%eps, st, ierr))
383
384
3/4
✓ Branch 109 → 110 taken 16 times.
✓ Branch 109 → 111 taken 6 times.
✗ Branch 110 → 111 not taken.
✓ Branch 110 → 115 taken 16 times.
 22 if (this%options%which_is_target .or. this%options%which == EPS_ALL) then
385 ! Target-based eigenvalue problems require shift-and-invert
386
1/2
✗ Branch 112 → 113 not taken.
✓ Branch 112 → 115 taken 6 times.
 6 SlepcCall(STSetType(st, STSINVERT, ierr))
387 end if
388
389
1/2
✗ Branch 116 → 117 not taken.
✓ Branch 116 → 119 taken 22 times.
 22 SlepcCall(STGetKSP(st, ksp, ierr))
390
1/2
✗ Branch 120 → 121 not taken.
✓ Branch 120 → 122 taken 22 times.
 22 PetscCall(KSPSetType(ksp, this%options%ksp_type, ierr))
391
1/2
✗ Branch 123 → 124 not taken.
✓ Branch 123 → 125 taken 22 times.
 22 PetscCall(KSPGetPC(ksp, pc, ierr))
392
2/4
✓ Branch 125 → 126 taken 22 times.
✗ Branch 125 → 127 not taken.
✓ Branch 126 → 127 taken 22 times.
✗ Branch 126 → 134 not taken.
 22 if (this%options%pc_type == 'cholesky' .or. this%options%pc_type == 'lu') then
393
1/2
✗ Branch 128 → 129 not taken.
✓ Branch 128 → 130 taken 22 times.
 22 PetscCall(PCSetType(pc, this%options%pc_type, ierr))
394
1/2
✗ Branch 130 → 131 not taken.
✓ Branch 130 → 143 taken 22 times.
 22 if (traits%is_parallel) then
395 ✗ PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS, ierr))
396 end if
397 ✗ else if (.not. (this%options%pc_type == 'hypre' .and. this%options%hypre_type == 'ilu')) then
398 ✗ PetscCall(PCSetType(pc, this%options%pc_type, ierr))
399 ✗ if (this%options%pc_type == 'hypre') then
400 ✗ PetscCall(PCHYPRESetType(pc, this%options%hypre_type, ierr))
401 end if
402 end if
403
404 ! Allow command-line overrides after programmatic configuration
405
1/2
✗ Branch 144 → 145 not taken.
✓ Branch 144 → 147 taken 22 times.
 22 SlepcCall(EPSSetFromOptions(this%eps, ierr))
406 22 end subroutine EigenSolver_init
407
408 !> @brief Solve the eigenvalue problem
409 !>
410 !> @param this The EigenSolver object
411 !> @param rtol _(optional)_ The relative tolerance (default: 1.e-8)
412 !> @param max_iter _(optional)_ The maximum number of iterations (default: size of the space)
413 !> @param nr_pairs _(optional)_ The number of eigenpairs to compute (default: 6)
414 !> @param verbosity _(optional)_ The verbosity level (default: VERBOSITY_WARN_ON_FAILURE)
415 !> - VERBOSITY_WARN_NEVER
416 !> - VERBOSITY_WARN_ON_FAILURE
417 !> - VERBOSITY_WARN_ALWAYS
418 !> - VERBOSITY_ERROR_ON_FAILURE
419 22 subroutine EigenSolver_solve(this, rtol, max_iter, nr_pairs, verbosity)
420 use m_common, only: wp, VERBOSITY_WARN_NEVER, VERBOSITY_WARN_ON_FAILURE, VERBOSITY_WARN_ALWAYS, VERBOSITY_ERROR_ON_FAILURE
421 implicit none
422
423 class(EigenSolver), intent(inout) :: this
424 integer :: ierr
425 real(wp), optional, intent(in) :: rtol
426 integer, optional, intent(in) :: max_iter
427 integer, optional, intent(in) :: nr_pairs
428 integer, optional, intent(in) :: verbosity
429
430 real(wp) :: rtol_
431 integer :: max_iter_
432 integer :: nr_pairs_
433 integer :: verbosity_
434
435
1/2
✓ Branch 2 → 3 taken 22 times.
✗ Branch 2 → 4 not taken.
 22 if (present(rtol)) then
436 22 rtol_ = rtol
437 else
438 ✗ rtol_ = 1.e-8_wp
439 end if
440
441
1/2
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 22 times.
 22 if (present(max_iter)) then
442 ✗ max_iter_ = max_iter
443 else
444 22 max_iter_ = size(this%space_row)
445 end if
446
447
1/2
✓ Branch 8 → 9 taken 22 times.
✗ Branch 8 → 10 not taken.
 22 if (present(nr_pairs)) then
448 22 nr_pairs_ = nr_pairs
449 else
450 ✗ nr_pairs_ = min(6, size(this%space_row))
451 end if
452
453 22 if (this%space_row%tp_spaces(1)%domain%my_rank > 0) then
454 verbosity_ = VERBOSITY_WARN_NEVER
455 else
456 if (present(verbosity)) then
457 verbosity_ = verbosity
458 else
459 verbosity_ = VERBOSITY_WARN_ON_FAILURE
460 end if
461 end if
462
463
1/2
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 15 taken 22 times.
 22 SlepcCall(EPSSetTolerances(this%eps, rtol_, max_iter_, ierr))
464
1/2
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 19 taken 22 times.
 22 SlepcCall(EPSSetDimensions(this%eps, nr_pairs_, PETSC_DECIDE, PETSC_DECIDE, ierr))
465
466
1/2
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 23 taken 22 times.
 22 SlepcCall(EPSSolve(this%eps, ierr))
467
1/2
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 27 taken 22 times.
 22 SlepcCall(EPSGetConverged(this%eps, this%nr_pairs, ierr))
468 ! TODO check for convergence
469 22 end subroutine EigenSolver_solve
470
471 !> @brief Get an eigenpair
472 !>
473 !> @param this The EigenSolver object
474 !> @param index The index of the eigenpair to get (1-based)
475 !> @param val_r The real part of the eigenvalue
476 !> @param fun_r The real part of the eigenfunction
477 !> @param val_i _(optional)_ The imaginary part of the eigenvalue
478 !> @param fun_i _(optional)_ The imaginary part of the eigenfunction
479 !> @param err_est _(optional)_ The error estimate for the eigenpair
480 ✗ subroutine EigenSolver_get_eigenpair(this, index, val_r, fun_r, val_i, fun_i, err_est)
481 use m_mform_constraint_global, only: CONSTRAINT_EXTRACTION
482 use m_mform_basis, only: MFormFun
483 implicit none
484
485 class(EigenSolver), intent(in) :: this
486 integer, intent(in) :: index
487 real(wp), intent(out) :: val_r
488 real(wp), optional, intent(out) :: val_i
489 type(MFormFun), intent(out) :: fun_r
490 type(MFormFun), optional, intent(out) :: fun_i
491 real(wp), optional, intent(out) :: err_est
492
493 Vec :: vec_r, vec_i
494 Vec :: vec_r_proj, vec_i_proj
495 integer :: ierr
496
497 ✗ PetscCall(MatCreateVecs(this%matA, vec_r_proj, PETSC_NULL_VEC, ierr))
498 ✗ if (present(fun_i)) then
499 ✗ PetscCall(VecDuplicate(vec_r_proj, vec_i_proj, ierr))
500 end if
501
502 if (present(fun_i)) then
503 ✗ SlepcCall(EPSGetEigenpair(this%eps, index - 1, val_r, val_i, vec_r_proj, vec_i_proj, ierr))
504 else
505 ✗ SlepcCall(EPSGetEigenpair(this%eps, index - 1, val_r, PETSC_NULL_SCALAR, vec_r_proj, PETSC_NULL_VEC, ierr))
506 end if
507
508 ✗ if (allocated(this%constraint)) then
509 ✗ if (this%constraint%mode == CONSTRAINT_EXTRACTION) then
510 ✗ PetscCall(VecCreate(this%space_row%tp_spaces(1)%domain%comm, vec_r, ierr))
511 ✗ PetscCall(VecSetSizes(vec_r, size(this%space_row, my_rank=.true.), size(this%space_row), ierr))
512 ✗ PetscCall(VecSetFromOptions(vec_r, ierr))
513 ✗ PetscCall(MatMult(this%constraint%mat, vec_r_proj, vec_r, ierr))
514
515 ✗ if (present(fun_i)) then
516 ✗ PetscCall(VecDuplicate(vec_r, vec_i, ierr))
517 ✗ PetscCall(MatMult(this%constraint%mat, vec_i_proj, vec_i, ierr))
518 end if
519 else
520 ✗ PetscCall(VecCopy(vec_r_proj, vec_r, ierr))
521 ✗ if (present(fun_i)) then
522 ✗ PetscCall(VecCopy(vec_i_proj, vec_i, ierr))
523 end if
524 end if
525 else
526 ✗ vec_r = vec_r_proj
527 ✗ if (present(fun_i)) then
528 ✗ vec_i = vec_i_proj
529 end if
530 end if
531
532 ✗ call fun_r%init(this%space_row, vec=vec_r)
533 ✗ PetscCall(VecDestroy(vec_r, ierr))
534 ✗ if (allocated(this%constraint)) then
535 ✗ PetscCall(VecDestroy(vec_r_proj, ierr))
536 end if
537
538 ✗ if (present(fun_i)) then
539 ✗ call fun_i%init(this%space_row, vec=vec_i)
540 ✗ PetscCall(VecDestroy(vec_i, ierr))
541 ✗ if (allocated(this%constraint)) then
542 ✗ PetscCall(VecDestroy(vec_i_proj, ierr))
543 end if
544 end if
545
546 ✗ if (present(err_est)) then
547 ✗ SlepcCall(EPSGetErrorEstimate(this%eps, index - 1, err_est, ierr))
548 end if
549
550 ✗ end subroutine EigenSolver_get_eigenpair
551
552 !> @brief Get all computed eigenvalues
553 !>
554 !> @param this The EigenSolver object
555 !> @param vals_r The real parts of the eigenvalues
556 !> @param vals_i _(optional)_ The imaginary parts of the eigenvalues
557 22 subroutine EigenSolver_get_eigenvalues(this, vals_r, vals_i)
558 implicit none
559
560 class(EigenSolver), intent(in) :: this
561 real(wp), allocatable, intent(out) :: vals_r(:)
562 real(wp), optional, allocatable, intent(out) :: vals_i(:)
563
564 integer :: i, ierr
565
566
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 14 taken 22 times.
 22 if (this%nr_pairs <= 0) then
567 ✗ allocate (vals_r(0))
568 ✗ if (present(vals_i)) allocate (vals_i(0))
569 ✗ return
570 end if
571
572
6/12
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 22 times.
✓ Branch 16 → 15 taken 22 times.
✗ Branch 16 → 17 not taken.
✓ Branch 17 → 18 taken 22 times.
✗ Branch 17 → 19 not taken.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 22 times.
✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 22 times.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 22 times.
 66 allocate (vals_r(this%nr_pairs))
573
6/12
✓ Branch 25 → 26 taken 22 times.
✗ Branch 25 → 42 not taken.
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 22 times.
✓ Branch 28 → 27 taken 22 times.
✗ Branch 28 → 29 not taken.
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 22 times.
✗ Branch 31 → 32 not taken.
✓ Branch 31 → 33 taken 22 times.
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 22 times.
 44 if (present(vals_i)) allocate (vals_i(this%nr_pairs))
574
575 if (present(vals_i)) then
576
2/2
✓ Branch 36 → 37 taken 295 times.
✓ Branch 36 → 49 taken 22 times.
 317 do i = 1, this%nr_pairs
577
1/2
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 41 taken 295 times.
 317 SlepcCall(EPSGetEigenvalue(this%eps, i - 1, vals_r(i), vals_i(i), ierr))
578 end do
579 else
580 ✗ do i = 1, this%nr_pairs
581 ✗ SlepcCall(EPSGetEigenvalue(this%eps, i - 1, vals_r(i), PETSC_NULL_SCALAR, ierr))
582 end do
583 end if
584
585 end subroutine EigenSolver_get_eigenvalues
586
587 !> @brief Destroy the eigensolver and free associated resources
588 !>
589 !> @param this The EigenSolver object
590 25 subroutine EigenSolver_destroy(this)
591 implicit none
592
593 class(EigenSolver), intent(inout) :: this
594 integer :: ierr
595
596
1/2
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 6 taken 25 times.
 25 SlepcCall(EPSDestroy(this%eps, ierr))
597
1/2
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 25 times.
 25 PetscCall(MatDestroy(this%matA, ierr))
598
1/2
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 25 times.
 25 PetscCall(MatDestroy(this%matB, ierr))
599
2/2
✓ Branch 12 → 13 taken 14 times.
✓ Branch 12 → 87 taken 11 times.
 25 if (allocated(this%constraint)) then
600 14 call this%constraint%destroy()
601
12/68
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 14 times.
✓ Branch 16 → 17 taken 14 times.
✗ Branch 16 → 18 not taken.
✓ Branch 19 → 20 taken 70 times.
✓ Branch 19 → 63 taken 14 times.
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 23 taken 70 times.
✗ Branch 21 → 22 not taken.
✗ Branch 21 → 23 not taken.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 62 taken 70 times.
✗ Branch 24 → 25 not taken.
✗ Branch 24 → 26 not taken.
✗ Branch 26 → 27 not taken.
✗ Branch 26 → 39 not taken.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 38 not taken.
✗ Branch 29 → 30 not taken.
✗ Branch 29 → 31 not taken.
✗ Branch 31 → 32 not taken.
✗ Branch 31 → 33 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 39 → 40 not taken.
✗ Branch 39 → 50 not taken.
✗ Branch 41 → 42 not taken.
✗ Branch 41 → 49 not taken.
✗ Branch 42 → 43 not taken.
✗ Branch 42 → 44 not taken.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 46 → 47 not taken.
✗ Branch 46 → 48 not taken.
✗ Branch 50 → 51 not taken.
✗ Branch 50 → 61 not taken.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 60 not taken.
✗ Branch 53 → 54 not taken.
✗ Branch 53 → 55 not taken.
✗ Branch 55 → 56 not taken.
✗ Branch 55 → 57 not taken.
✗ Branch 57 → 58 not taken.
✗ Branch 57 → 59 not taken.
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 65 taken 14 times.
✗ Branch 65 → 66 not taken.
✓ Branch 65 → 67 taken 14 times.
✗ Branch 67 → 68 not taken.
✓ Branch 67 → 69 taken 14 times.
✗ Branch 69 → 70 not taken.
✓ Branch 69 → 71 taken 14 times.
✗ Branch 71 → 72 not taken.
✓ Branch 71 → 73 taken 14 times.
✗ Branch 73 → 74 not taken.
✓ Branch 73 → 86 taken 14 times.
✗ Branch 75 → 76 not taken.
✗ Branch 75 → 85 not taken.
✗ Branch 76 → 77 not taken.
✗ Branch 76 → 78 not taken.
✗ Branch 78 → 79 not taken.
✗ Branch 78 → 80 not taken.
✗ Branch 80 → 81 not taken.
✗ Branch 80 → 82 not taken.
✗ Branch 82 → 83 not taken.
✗ Branch 82 → 84 not taken.
 98 deallocate (this%constraint)
602 end if
603
604 end subroutine EigenSolver_destroy
605
606 ✗ end module m_mform_eigensolver
607