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