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