Reference documentation for deal.II version 9.3.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\}}\)
slepc_solver.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2020 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 
33 namespace 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,
47  this->solver_control.tolerance(),
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
99  SLEPcWrappers::TransformationBase &transformation)
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.
114  dynamic_cast<SLEPcWrappers::TransformationShiftInvert *>(
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),
291  ExcMessage(
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 ----------------------- */
334  SolverControl &
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 
384  const MPI_Comm & mpi_communicator,
385  const AdditionalData &data)
386  : SolverBase(cn, mpi_communicator)
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 
411  const MPI_Comm & mpi_communicator,
412  const AdditionalData &data)
413  : SolverBase(cn, mpi_communicator)
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,
448  const MPI_Comm & mpi_communicator,
449  const AdditionalData &data)
450  : SolverBase(cn, mpi_communicator)
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),
493  ExcMessage(
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
void set_problem_type(EPSProblemType set_problem)
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1528
static int convergence_test(EPS eps, PetscScalar real_eigenvalue, PetscScalar imag_eigenvalue, PetscReal residual_norm, PetscReal *estimated_error, void *solver_control)
virtual State check(const unsigned int step, const double check_value)
SolverJacobiDavidson(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
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)
SolverPower(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
void set_which_eigenpairs(EPSWhich set_which)
void set_target_eigenvalue(const PetscScalar &this_target)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:302
void get_eigenpair(const unsigned int index, PetscScalar &eigenvalues, PETScWrappers::VectorBase &eigenvectors)
SolverKrylovSchur(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
AdditionalData(const bool delayed_reorthogonalization=false)
static ::ExceptionBase & ExcMessage(std::string arg1)
void get_solver_state(const SolverControl::State state)
void set_matrices(const PETScWrappers::MatrixBase &A)
Definition: slepc_solver.cc:77
#define Assert(cond, exc)
Definition: exceptions.h:1465
SolverGeneralizedDavidson(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
EPSConvergedReason reason
Definition: slepc_solver.h:360
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:518
SolverLanczos(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
Definition: slepc_solver.cc:98
double tolerance() const
unsigned int last_step() const
static const char A
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
unsigned int max_steps() const
SolverLAPACK(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
Definition: slepc_solver.cc:35
const AdditionalData additional_data
Definition: slepc_solver.h:649
static ::ExceptionBase & ExcNotImplemented()
AdditionalData(const EPSLanczosReorthogType r=EPS_LANCZOS_REORTHOG_FULL)
const AdditionalData additional_data
Definition: slepc_solver.h:466
static ::ExceptionBase & ExcSLEPcError(int arg1)
void solve(const PETScWrappers::MatrixBase &A, std::vector< PetscScalar > &eigenvalues, std::vector< OutputVector > &eigenvectors, const unsigned int n_eigenpairs=1)
Definition: slepc_solver.h:704
Eigenvalue vector is filled.
SolverControl & solver_control
Definition: slepc_solver.h:297
SolverArnoldi(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())