Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
slepc_solver.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2018 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 
16 #include <deal.II/lac/slepc_solver.h>
17 
18 #ifdef DEAL_II_WITH_SLEPC
19 
20 # include <deal.II/lac/petsc_matrix_base.h>
21 # include <deal.II/lac/petsc_vector_base.h>
22 # include <deal.II/lac/slepc_spectral_transformation.h>
23 
24 # include <petscversion.h>
25 
26 # include <slepcversion.h>
27 
28 # include <cmath>
29 # include <vector>
30 
31 DEAL_II_NAMESPACE_OPEN
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 
61  {
62  if (eps != nullptr)
63  {
64  // Destroy the solver object.
65  const PetscErrorCode ierr = EPSDestroy(&eps);
66 
67  (void)ierr;
68  AssertNothrow(ierr == 0, ExcSLEPcError(ierr));
69  }
70  }
71 
72  void
74  {
75  // standard eigenspectrum problem
76  const PetscErrorCode ierr = EPSSetOperators(eps, A, nullptr);
77  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
78  }
79 
80  void
83  {
84  // generalized eigenspectrum problem
85  const PetscErrorCode ierr = EPSSetOperators(eps, A, B);
86  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
87  }
88 
89  void
91  SLEPcWrappers::TransformationBase &transformation)
92  {
93  // set transformation type if any
94  // STSetShift is called inside
95  PetscErrorCode ierr = EPSSetST(eps, transformation.st);
96  AssertThrow(ierr == 0, SolverBase::ExcSLEPcError(ierr));
97 
98 # if DEAL_II_SLEPC_VERSION_GTE(3, 8, 0)
99  // see
100  // https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2017-October/033649.html
101  // From 3.8.0 SLEPc insists that when looking for smallest eigenvalues with
102  // shift-and-invert users should (a) set target (b) use EPS_TARGET_MAGNITUDE
103  // The former, however, needs to be applied to eps object and not spectral
104  // transformation.
106  dynamic_cast<SLEPcWrappers::TransformationShiftInvert *>(
107  &transformation))
108  {
109  ierr = EPSSetTarget(eps, sinv->additional_data.shift_parameter);
110  AssertThrow(ierr == 0, SolverBase::ExcSLEPcError(ierr));
111  }
112 # endif
113  }
114 
115  void
117  const PETScWrappers::VectorBase &this_initial_vector)
118  {
119  Assert(this_initial_vector.l2_norm() > 0.0,
120  ExcMessage("Initial vector should be nonzero."));
121 
122  Vec vec = this_initial_vector;
123  const PetscErrorCode ierr = EPSSetInitialSpace(eps, 1, &vec);
124 
125  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
126  }
127 
128  void
129  SolverBase::set_target_eigenvalue(const PetscScalar &this_target)
130  {
131  // set target eigenvalues to solve for
132  // in all transformation except STSHIFT there is a direct connection between
133  // the target and the shift, read more on p41 of SLEPc manual.
134  const PetscErrorCode ierr = EPSSetTarget(eps, this_target);
135  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
136  }
137 
138  void
139  SolverBase::set_which_eigenpairs(const EPSWhich eps_which)
140  {
141  // set which portion of the eigenspectrum to solve for
142  const PetscErrorCode ierr = EPSSetWhichEigenpairs(eps, eps_which);
143  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
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  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
185  ierr = EPSGetConverged(eps, reinterpret_cast<PetscInt *>(n_converged));
186  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
187 
188  PetscInt n_iterations = 0;
189  PetscReal residual_norm = 0;
190 
191  // @todo Investigate elaborating on some of this to act on the
192  // complete eigenspectrum
193  {
194  // get the number of solver iterations
195  ierr = EPSGetIterationNumber(eps, &n_iterations);
196  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
197 
198  // get the maximum of residual norm among converged eigenvectors.
199  for (unsigned int i = 0; i < *n_converged; i++)
200  {
201  double residual_norm_i = 0.0;
202  // EPSComputeError (or, in older versions of SLEPc,
203  // EPSComputeResidualNorm) uses an L2-norm and is not consistent
204  // with the stopping criterion used during the solution process (see
205  // the SLEPC manual, section 2.5). However, the norm that gives error
206  // bounds (Saad, 1992, ch3) is (for Hermitian problems)
207  // | \lambda - \widehat\lambda | <= ||r||_2
208  //
209  // Similarly, EPSComputeRelativeError may not be consistent with the
210  // stopping criterion used in the solution process.
211  //
212  // EPSGetErrorEstimate is (according to the SLEPc manual) consistent
213  // with the residual norm used during the solution process. However,
214  // it is not guaranteed to be derived from the residual even when
215  // EPSSetTrueResidual is set: see the discussion in the thread
216  //
217  // https://lists.mcs.anl.gov/pipermail/petsc-users/2014-November/023509.html
218  //
219  // for more information.
220 # if DEAL_II_SLEPC_VERSION_GTE(3, 6, 0)
221  ierr = EPSComputeError(eps, i, EPS_ERROR_ABSOLUTE, &residual_norm_i);
222 # else
223  ierr = EPSComputeResidualNorm(eps, i, &residual_norm_i);
224 # endif
225 
226  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
227  residual_norm = std::max(residual_norm, residual_norm_i);
228  }
229 
230  // check the solver state
231  const SolverControl::State state =
232  solver_control.check(n_iterations, residual_norm);
233 
234  // get the solver state according to SLEPc
235  get_solver_state(state);
236 
237  // as SLEPc uses different stopping criteria, we have to omit this step.
238  // This can be checked only in conjunction with EPSGetErrorEstimate.
239  // and in case of failure: throw exception
240  // if (solver_control.last_check () != SolverControl::success)
241  // AssertThrow(false, SolverControl::NoConvergence
242  // (solver_control.last_step(),
243  // solver_control.last_value()));
244  }
245  }
246 
247  void
248  SolverBase::get_eigenpair(const unsigned int index,
249  PetscScalar & eigenvalues,
251  {
252  // get converged eigenpair
253  const PetscErrorCode ierr =
254  EPSGetEigenpair(eps, index, &eigenvalues, nullptr, eigenvectors, nullptr);
255  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
256  }
257 
258 
259  void
260  SolverBase::get_eigenpair(const unsigned int index,
261  double & real_eigenvalues,
262  double & imag_eigenvalues,
263  PETScWrappers::VectorBase &real_eigenvectors,
264  PETScWrappers::VectorBase &imag_eigenvectors)
265  {
266 # ifndef PETSC_USE_COMPLEX
267  // get converged eigenpair
268  const PetscErrorCode ierr = EPSGetEigenpair(eps,
269  index,
270  &real_eigenvalues,
271  &imag_eigenvalues,
272  real_eigenvectors,
273  imag_eigenvectors);
274  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
275 # else
276  Assert(
277  (false),
278  ExcMessage(
279  "Your PETSc/SLEPc installation was configured with scalar-type complex "
280  "but this function is not defined for complex types."));
281 
282  // Cast to void to silence compiler warnings
283  (void)index;
284  (void)real_eigenvalues;
285  (void)imag_eigenvalues;
286  (void)real_eigenvectors;
287  (void)imag_eigenvectors;
288 # endif
289  }
290 
291  void
293  {
294  switch (state)
295  {
296  case ::SolverControl::iterate:
297  reason = EPS_CONVERGED_ITERATING;
298  break;
299 
300  case ::SolverControl::success:
301  reason = static_cast<EPSConvergedReason>(1);
302  break;
303 
304  case ::SolverControl::failure:
306  reason = EPS_DIVERGED_ITS;
307  else
308  reason = EPS_DIVERGED_BREAKDOWN;
309  break;
310 
311  default:
312  Assert(false, ExcNotImplemented());
313  }
314  }
315 
316  /* ---------------------- SolverControls ----------------------- */
317  SolverControl &
319  {
320  return solver_control;
321  }
322 
323  int
325  EPS /*eps */,
326  PetscScalar /*real_eigenvalue */,
327  PetscScalar /*imag_eigenvalue */,
328  PetscReal /*residual norm associated to the eigenpair */,
329  PetscReal * /*(output) computed error estimate */,
330  void * /*solver_control_x*/)
331  {
332  // If the error estimate returned by the convergence test function is less
333  // than the tolerance, then the eigenvalue is accepted as converged.
334  // This function is undefined (future reference only).
335 
336  // return without failure.
337  return 0;
338  }
339 
340  /* ---------------------- SolverKrylovSchur ------------------------ */
342  const MPI_Comm & mpi_communicator,
343  const AdditionalData &data)
344  : SolverBase(cn, mpi_communicator)
345  , additional_data(data)
346  {
347  const PetscErrorCode ierr =
348  EPSSetType(eps, const_cast<char *>(EPSKRYLOVSCHUR));
349  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
350  }
351 
352  /* ---------------------- SolverArnoldi ------------------------ */
354  const bool delayed_reorthogonalization)
355  : delayed_reorthogonalization(delayed_reorthogonalization)
356  {}
357 
359  const MPI_Comm & mpi_communicator,
360  const AdditionalData &data)
362  , additional_data(data)
363  {
364  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSARNOLDI));
365  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
366 
367  // if requested, set delayed reorthogonalization in the Arnoldi
368  // iteration.
370  {
371  ierr = EPSArnoldiSetDelayed(eps, PETSC_TRUE);
372  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
373  }
374  }
375 
376 
377  /* ---------------------- Lanczos ------------------------ */
378  SolverLanczos::AdditionalData::AdditionalData(const EPSLanczosReorthogType r)
379  : reorthog(r)
380  {}
381 
383  const MPI_Comm & mpi_communicator,
384  const AdditionalData &data)
386  , additional_data(data)
387  {
388  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSLANCZOS));
389  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
390 
391  ierr = EPSLanczosSetReorthog(eps, additional_data.reorthog);
392  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
393  }
394 
395  /* ----------------------- Power ------------------------- */
397  const MPI_Comm & mpi_communicator,
398  const AdditionalData &data)
399  : SolverBase(cn, mpi_communicator)
400  , additional_data(data)
401  {
402  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSPOWER));
403  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
404  }
405 
406  /* ---------------- Generalized Davidson ----------------- */
408  bool double_expansion)
409  : double_expansion(double_expansion)
410  {}
411 
413  SolverControl & cn,
414  const MPI_Comm & mpi_communicator,
415  const AdditionalData &data)
417  , additional_data(data)
418  {
419  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSGD));
420  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
421 
423  {
424  ierr = EPSGDSetDoubleExpansion(eps, PETSC_TRUE);
425  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
426  }
427  }
428 
429  /* ------------------ Jacobi Davidson -------------------- */
431  const MPI_Comm &mpi_communicator,
432  const AdditionalData &data)
433  : SolverBase(cn, mpi_communicator)
434  , additional_data(data)
435  {
436  const PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSJD));
437  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
438  }
439 
440  /* ---------------------- LAPACK ------------------------- */
442  const MPI_Comm & mpi_communicator,
443  const AdditionalData &data)
444  : SolverBase(cn, mpi_communicator)
445  , additional_data(data)
446  {
447  // 'Tis overwhelmingly likely that PETSc/SLEPc *always* has
448  // BLAS/LAPACK, but let's be defensive.
449 # if PETSC_HAVE_BLASLAPACK
450  const PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSLAPACK));
451  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
452 # else
453  Assert(
454  (false),
455  ExcMessage(
456  "Your PETSc/SLEPc installation was not configured with BLAS/LAPACK "
457  "but this is needed to use the LAPACK solver."));
458 # endif
459  }
460 } // namespace SLEPcWrappers
461 
462 DEAL_II_NAMESPACE_CLOSE
463 
464 #endif // DEAL_II_WITH_SLEPC
void set_problem_type(EPSProblemType set_problem)
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1471
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())
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
void set_which_eigenpairs(EPSWhich set_which)
void set_target_eigenvalue(const PetscScalar &this_target)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:309
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:73
#define Assert(cond, exc)
Definition: exceptions.h:1407
void set_initial_vector(const PETScWrappers::VectorBase &this_initial_vector)
SolverGeneralizedDavidson(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
EPSConvergedReason reason
Definition: slepc_solver.h:367
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:507
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:90
double tolerance() const
unsigned int last_step() const
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
static ::ExceptionBase & ExcNotImplemented()
AdditionalData(const EPSLanczosReorthogType r=EPS_LANCZOS_REORTHOG_FULL)
const AdditionalData additional_data
Definition: slepc_solver.h:461
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:668
SolverControl & solver_control
Definition: slepc_solver.h:304
SolverArnoldi(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())