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