deal.II version GIT relicensing-3517-g3d7778a52c 2025-06-18 11:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
slepc_solver.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
17#ifdef DEAL_II_WITH_SLEPC
18
22
23# include <petscversion.h>
24
25# include <slepcversion.h>
26
27# include <cmath>
28# include <vector>
29
31
32namespace SLEPcWrappers
33{
34 SolverBase::SolverBase(SolverControl &cn, const MPI_Comm mpi_communicator)
35 : solver_control(cn)
36 , mpi_communicator(mpi_communicator)
37 , reason(EPS_CONVERGED_ITERATING)
38 {
39 // create eigensolver context
40 PetscErrorCode ierr = EPSCreate(mpi_communicator, &eps);
41 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
42
43 // hand over the absolute tolerance and the maximum number of
44 // iteration steps to the SLEPc convergence criterion.
45 ierr = EPSSetTolerances(eps,
47 this->solver_control.max_steps());
48 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
49
50 // default values:
51 set_which_eigenpairs(EPS_LARGEST_MAGNITUDE);
52 set_problem_type(EPS_GNHEP);
53
54 // TODO:
55 // By default, EPS initializes the starting vector or the initial subspace
56 // randomly.
57 }
58
59
60
62 {
63 if (eps != nullptr)
64 {
65 // Destroy the solver object.
66 const PetscErrorCode ierr = EPSDestroy(&eps);
67 AssertNothrow(ierr == 0, ExcSLEPcError(ierr));
68 }
69 }
70
71
72
73 void
75 {
76 // standard eigenspectrum problem
77 const PetscErrorCode ierr = EPSSetOperators(eps, A, nullptr);
78 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
79 }
80
81
82
83 void
86 {
87 // generalized eigenspectrum problem
88 const PetscErrorCode ierr = EPSSetOperators(eps, A, B);
89 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
90 }
91
92
93
94 void
97 {
98 // set transformation type if any
99 // STSetShift is called inside
100 PetscErrorCode ierr = EPSSetST(eps, transformation.st);
101 AssertThrow(ierr == 0, SolverBase::ExcSLEPcError(ierr));
102
103# if DEAL_II_SLEPC_VERSION_GTE(3, 8, 0)
104 // see
105 // https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2017-October/033649.html
106 // From 3.8.0 onward, SLEPc insists that when looking for smallest
107 // eigenvalues with shift-and-invert, users should (a) set target,
108 // (b) use EPS_TARGET_MAGNITUDE. The former, however, needs to be
109 // applied to the 'eps' object and not the spectral transformation.
112 &transformation))
113 {
114 ierr = EPSSetTarget(eps, sinv->additional_data.shift_parameter);
115 AssertThrow(ierr == 0, SolverBase::ExcSLEPcError(ierr));
116 }
117# endif
118 }
119
120
121
122 void
123 SolverBase::set_target_eigenvalue(const PetscScalar &this_target)
124 {
125 // set target eigenvalues to solve for
126 // in all transformation except STSHIFT there is a direct connection between
127 // the target and the shift, read more on p41 of SLEPc manual.
128 const PetscErrorCode ierr = EPSSetTarget(eps, this_target);
129 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
130 }
131
132
133
134 void
135 SolverBase::set_which_eigenpairs(const EPSWhich eps_which)
136 {
137 // set which portion of the eigenspectrum to solve for
138 const PetscErrorCode ierr = EPSSetWhichEigenpairs(eps, eps_which);
139 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
140 }
141
142
143
144 void
145 SolverBase::set_problem_type(const EPSProblemType eps_problem)
146 {
147 const PetscErrorCode ierr = EPSSetProblemType(eps, eps_problem);
148 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
149 }
150
151
152
153 void
154 SolverBase::solve(const unsigned int n_eigenpairs, unsigned int *n_converged)
155 {
156 // set number of eigenvectors to compute
157 PetscErrorCode ierr =
158 EPSSetDimensions(eps, n_eigenpairs, PETSC_DECIDE, PETSC_DECIDE);
159 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
160
161 // set the solve options to the eigenvalue problem solver context
162 ierr = EPSSetFromOptions(eps);
163 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
164
165 // TODO breaks @ref step_36 "step-36"
166 // force Krylov solver to use true residual instead of an estimate.
167 // EPSSetTrueResidual(solver_data->eps, PETSC_TRUE);
168 // AssertThrow (ierr == 0, ExcSLEPcError(ierr));
169
170 // Set convergence test to be absolute
171 ierr = EPSSetConvergenceTest(eps, EPS_CONV_ABS);
172 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
173
174 // TODO Set the convergence test function
175 // ierr = EPSSetConvergenceTestFunction (solver_data->eps,
176 // &convergence_test,
177 // reinterpret_cast<void *>(&solver_control));
178 // AssertThrow (ierr == 0, ExcSLEPcError(ierr));
179
180 // solve the eigensystem
181 ierr = EPSSolve(eps);
182 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
183
184 // Get number of converged eigenstates. We need to go around with a
185 // temporary variable once because the function wants to have a
186 // PetscInt as second argument whereas the `n_converged` argument
187 // to this function is just an unsigned int.
188 {
189 PetscInt petsc_n_converged = *n_converged;
190 ierr = EPSGetConverged(eps, &petsc_n_converged);
191 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
192 *n_converged = petsc_n_converged;
193 }
194
195 PetscInt n_iterations = 0;
196 PetscReal residual_norm = 0;
197
198 // @todo Investigate elaborating on some of this to act on the
199 // complete eigenspectrum
200 {
201 // get the number of solver iterations
202 ierr = EPSGetIterationNumber(eps, &n_iterations);
203 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
204
205 // get the maximum of residual norm among converged eigenvectors.
206 for (unsigned int i = 0; i < *n_converged; ++i)
207 {
208 double residual_norm_i = 0.0;
209 // EPSComputeError (or, in older versions of SLEPc,
210 // EPSComputeResidualNorm) uses an L2-norm and is not consistent
211 // with the stopping criterion used during the solution process (see
212 // the SLEPC manual, section 2.5). However, the norm that gives error
213 // bounds (Saad, 1992, ch3) is (for Hermitian problems)
214 // | \lambda - \widehat\lambda | <= ||r||_2
215 //
216 // Similarly, EPSComputeRelativeError may not be consistent with the
217 // stopping criterion used in the solution process.
218 //
219 // EPSGetErrorEstimate is (according to the SLEPc manual) consistent
220 // with the residual norm used during the solution process. However,
221 // it is not guaranteed to be derived from the residual even when
222 // EPSSetTrueResidual is set: see the discussion in the thread
223 //
224 // https://lists.mcs.anl.gov/pipermail/petsc-users/2014-November/023509.html
225 //
226 // for more information.
227# if DEAL_II_SLEPC_VERSION_GTE(3, 6, 0)
228 ierr = EPSComputeError(eps, i, EPS_ERROR_ABSOLUTE, &residual_norm_i);
229# else
230 ierr = EPSComputeResidualNorm(eps, i, &residual_norm_i);
231# endif
232
233 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
234 residual_norm = std::max(residual_norm, residual_norm_i);
235 }
236
237 // check the solver state
238 const SolverControl::State state =
239 solver_control.check(n_iterations, residual_norm);
240
241 // get the solver state according to SLEPc
242 get_solver_state(state);
243
244 // as SLEPc uses different stopping criteria, we have to omit this step.
245 // This can be checked only in conjunction with EPSGetErrorEstimate.
246 // and in case of failure: throw exception
247 // if (solver_control.last_check () != SolverControl::success)
248 // AssertThrow(false, SolverControl::NoConvergence
249 // (solver_control.last_step(),
250 // solver_control.last_value()));
251 }
252 }
253
254
255
256 void
257 SolverBase::get_eigenpair(const unsigned int index,
258 PetscScalar &eigenvalues,
260 {
261 // get converged eigenpair
262 const PetscErrorCode ierr =
263 EPSGetEigenpair(eps, index, &eigenvalues, nullptr, eigenvectors, nullptr);
264 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
265 }
266
267
268
269 void
270 SolverBase::get_eigenpair(const unsigned int index,
271 double &real_eigenvalues,
272 double &imag_eigenvalues,
273 PETScWrappers::VectorBase &real_eigenvectors,
274 PETScWrappers::VectorBase &imag_eigenvectors)
275 {
276# ifndef PETSC_USE_COMPLEX
277 // get converged eigenpair
278 const PetscErrorCode ierr = EPSGetEigenpair(eps,
279 index,
280 &real_eigenvalues,
281 &imag_eigenvalues,
282 real_eigenvectors,
283 imag_eigenvectors);
284 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
285# else
286 Assert(
287 (false),
289 "Your PETSc/SLEPc installation was configured with scalar-type complex "
290 "but this function is not defined for complex types."));
291
292 // Cast to void to silence compiler warnings
293 (void)index;
294 (void)real_eigenvalues;
295 (void)imag_eigenvalues;
296 (void)real_eigenvectors;
297 (void)imag_eigenvectors;
298# endif
299 }
300
301
302
303 void
305 {
306 switch (state)
307 {
308 case ::SolverControl::iterate:
309 reason = EPS_CONVERGED_ITERATING;
310 break;
311
312 case ::SolverControl::success:
313 reason = static_cast<EPSConvergedReason>(1);
314 break;
315
316 case ::SolverControl::failure:
318 reason = EPS_DIVERGED_ITS;
319 else
320 reason = EPS_DIVERGED_BREAKDOWN;
321 break;
322
323 default:
325 }
326 }
327
328
329
330 /* ---------------------- SolverControls ----------------------- */
333 {
334 return solver_control;
335 }
336
337
338
339 int
341 EPS /*eps */,
342 PetscScalar /*real_eigenvalue */,
343 PetscScalar /*imag_eigenvalue */,
344 PetscReal /*residual norm associated to the eigenpair */,
345 PetscReal * /*(output) computed error estimate */,
346 void * /*solver_control_x*/)
347 {
348 // If the error estimate returned by the convergence test function is less
349 // than the tolerance, then the eigenvalue is accepted as converged.
350 // This function is undefined (future reference only).
351
352 // return without failure.
353 return 0;
354 }
355
356
357
358 /* ---------------------- SolverKrylovSchur ------------------------ */
360 const MPI_Comm mpi_communicator,
361 const AdditionalData &data)
362 : SolverBase(cn, mpi_communicator)
363 , additional_data(data)
364 {
365 const PetscErrorCode ierr =
366 EPSSetType(eps, const_cast<char *>(EPSKRYLOVSCHUR));
367 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
368 }
369
370
371
372 /* ---------------------- SolverArnoldi ------------------------ */
374 const bool delayed_reorthogonalization)
375 : delayed_reorthogonalization(delayed_reorthogonalization)
376 {}
377
378
379
382 const AdditionalData &data)
385 {
386 PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSARNOLDI));
387 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
388
389 // if requested, set delayed reorthogonalization in the Arnoldi
390 // iteration.
392 {
393 ierr = EPSArnoldiSetDelayed(eps, PETSC_TRUE);
394 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
395 }
396 }
397
398
399
400 /* ---------------------- Lanczos ------------------------ */
401 SolverLanczos::AdditionalData::AdditionalData(const EPSLanczosReorthogType r)
402 : reorthog(r)
403 {}
404
405
406
409 const AdditionalData &data)
412 {
413 PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSLANCZOS));
414 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
415
416 ierr = EPSLanczosSetReorthog(eps, additional_data.reorthog);
417 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
418 }
419
420
421
422 /* ----------------------- Power ------------------------- */
424 const MPI_Comm mpi_communicator,
425 const AdditionalData &data)
426 : SolverBase(cn, mpi_communicator)
427 , additional_data(data)
428 {
429 PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSPOWER));
430 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
431 }
432
433
434
435 /* ---------------- Generalized Davidson ----------------- */
437 bool double_expansion)
438 : double_expansion(double_expansion)
439 {}
440
441
442
444 SolverControl &cn,
446 const AdditionalData &data)
449 {
450 PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSGD));
451 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
452
454 {
455 ierr = EPSGDSetDoubleExpansion(eps, PETSC_TRUE);
456 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
457 }
458 }
459
460
461
462 /* ------------------ Jacobi Davidson -------------------- */
464 const MPI_Comm mpi_communicator,
465 const AdditionalData &data)
466 : SolverBase(cn, mpi_communicator)
467 , additional_data(data)
468 {
469 const PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSJD));
470 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
471 }
472
473
474
475 /* ---------------------- LAPACK ------------------------- */
477 const MPI_Comm mpi_communicator,
478 const AdditionalData &data)
479 : SolverBase(cn, mpi_communicator)
480 , additional_data(data)
481 {
482 // 'Tis overwhelmingly likely that PETSc/SLEPc *always* has
483 // BLAS/LAPACK, but let's be defensive.
484# if PETSC_HAVE_BLASLAPACK
485 const PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSLAPACK));
486 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
487# else
488 Assert(
489 (false),
491 "Your PETSc/SLEPc installation was not configured with BLAS/LAPACK "
492 "but this is needed to use the LAPACK solver."));
493# endif
494 }
495} // namespace SLEPcWrappers
496
498
499#endif // DEAL_II_WITH_SLEPC
SolverArnoldi(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverControl & control() const
EPSConvergedReason reason
void set_target_eigenvalue(const PetscScalar &this_target)
const MPI_Comm mpi_communicator
void set_matrices(const PETScWrappers::MatrixBase &A)
SolverBase(SolverControl &cn, const MPI_Comm mpi_communicator)
void get_solver_state(const SolverControl::State state)
void set_problem_type(EPSProblemType set_problem)
SolverControl & solver_control
void solve(const PETScWrappers::MatrixBase &A, std::vector< PetscScalar > &eigenvalues, std::vector< OutputVector > &eigenvectors, const unsigned int n_eigenpairs=1)
void get_eigenpair(const unsigned int index, PetscScalar &eigenvalues, PETScWrappers::VectorBase &eigenvectors)
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
static int convergence_test(EPS eps, PetscScalar real_eigenvalue, PetscScalar imag_eigenvalue, PetscReal residual_norm, PetscReal *estimated_error, void *solver_control)
void set_which_eigenpairs(EPSWhich set_which)
SolverGeneralizedDavidson(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
SolverJacobiDavidson(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
SolverKrylovSchur(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
SolverLAPACK(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverLanczos(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
SolverPower(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
unsigned int last_step() const
virtual State check(const unsigned int step, const double check_value)
unsigned int max_steps() const
double tolerance() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcSLEPcError(int arg1)
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:750
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
AdditionalData(const bool delayed_reorthogonalization=false)
AdditionalData(const EPSLanczosReorthogType r=EPS_LANCZOS_REORTHOG_FULL)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)