Reference documentation for deal.II version Git e1b67c38dc 2020-08-13 12:19:15 -0400
\(\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 
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
116  SolverBase::set_target_eigenvalue(const PetscScalar &this_target)
117  {
118  // set target eigenvalues to solve for
119  // in all transformation except STSHIFT there is a direct connection between
120  // the target and the shift, read more on p41 of SLEPc manual.
121  const PetscErrorCode ierr = EPSSetTarget(eps, this_target);
122  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
123  }
124 
125  void
126  SolverBase::set_which_eigenpairs(const EPSWhich eps_which)
127  {
128  // set which portion of the eigenspectrum to solve for
129  const PetscErrorCode ierr = EPSSetWhichEigenpairs(eps, eps_which);
130  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
131  }
132 
133  void
134  SolverBase::set_problem_type(const EPSProblemType eps_problem)
135  {
136  const PetscErrorCode ierr = EPSSetProblemType(eps, eps_problem);
137  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
138  }
139 
140  void
141  SolverBase::solve(const unsigned int n_eigenpairs, unsigned int *n_converged)
142  {
143  // set number of eigenvectors to compute
144  PetscErrorCode ierr =
145  EPSSetDimensions(eps, n_eigenpairs, PETSC_DECIDE, PETSC_DECIDE);
146  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
147 
148  // set the solve options to the eigenvalue problem solver context
149  ierr = EPSSetFromOptions(eps);
150  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
151 
152  // TODO breaks @ref step_36 "step-36"
153  // force Krylov solver to use true residual instead of an estimate.
154  // EPSSetTrueResidual(solver_data->eps, PETSC_TRUE);
155  // AssertThrow (ierr == 0, ExcSLEPcError(ierr));
156 
157  // Set convergence test to be absolute
158  ierr = EPSSetConvergenceTest(eps, EPS_CONV_ABS);
159  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
160 
161  // TODO Set the convergence test function
162  // ierr = EPSSetConvergenceTestFunction (solver_data->eps,
163  // &convergence_test,
164  // reinterpret_cast<void *>(&solver_control));
165  // AssertThrow (ierr == 0, ExcSLEPcError(ierr));
166 
167  // solve the eigensystem
168  ierr = EPSSolve(eps);
169  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
170 
171  // get number of converged eigenstates
172  ierr = EPSGetConverged(eps, reinterpret_cast<PetscInt *>(n_converged));
173  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
174 
175  PetscInt n_iterations = 0;
176  PetscReal residual_norm = 0;
177 
178  // @todo Investigate elaborating on some of this to act on the
179  // complete eigenspectrum
180  {
181  // get the number of solver iterations
182  ierr = EPSGetIterationNumber(eps, &n_iterations);
183  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
184 
185  // get the maximum of residual norm among converged eigenvectors.
186  for (unsigned int i = 0; i < *n_converged; i++)
187  {
188  double residual_norm_i = 0.0;
189  // EPSComputeError (or, in older versions of SLEPc,
190  // EPSComputeResidualNorm) uses an L2-norm and is not consistent
191  // with the stopping criterion used during the solution process (see
192  // the SLEPC manual, section 2.5). However, the norm that gives error
193  // bounds (Saad, 1992, ch3) is (for Hermitian problems)
194  // | \lambda - \widehat\lambda | <= ||r||_2
195  //
196  // Similarly, EPSComputeRelativeError may not be consistent with the
197  // stopping criterion used in the solution process.
198  //
199  // EPSGetErrorEstimate is (according to the SLEPc manual) consistent
200  // with the residual norm used during the solution process. However,
201  // it is not guaranteed to be derived from the residual even when
202  // EPSSetTrueResidual is set: see the discussion in the thread
203  //
204  // https://lists.mcs.anl.gov/pipermail/petsc-users/2014-November/023509.html
205  //
206  // for more information.
207 # if DEAL_II_SLEPC_VERSION_GTE(3, 6, 0)
208  ierr = EPSComputeError(eps, i, EPS_ERROR_ABSOLUTE, &residual_norm_i);
209 # else
210  ierr = EPSComputeResidualNorm(eps, i, &residual_norm_i);
211 # endif
212 
213  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
214  residual_norm = std::max(residual_norm, residual_norm_i);
215  }
216 
217  // check the solver state
218  const SolverControl::State state =
219  solver_control.check(n_iterations, residual_norm);
220 
221  // get the solver state according to SLEPc
222  get_solver_state(state);
223 
224  // as SLEPc uses different stopping criteria, we have to omit this step.
225  // This can be checked only in conjunction with EPSGetErrorEstimate.
226  // and in case of failure: throw exception
227  // if (solver_control.last_check () != SolverControl::success)
228  // AssertThrow(false, SolverControl::NoConvergence
229  // (solver_control.last_step(),
230  // solver_control.last_value()));
231  }
232  }
233 
234  void
235  SolverBase::get_eigenpair(const unsigned int index,
236  PetscScalar & eigenvalues,
238  {
239  // get converged eigenpair
240  const PetscErrorCode ierr =
241  EPSGetEigenpair(eps, index, &eigenvalues, nullptr, eigenvectors, nullptr);
242  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
243  }
244 
245 
246  void
247  SolverBase::get_eigenpair(const unsigned int index,
248  double & real_eigenvalues,
249  double & imag_eigenvalues,
250  PETScWrappers::VectorBase &real_eigenvectors,
251  PETScWrappers::VectorBase &imag_eigenvectors)
252  {
253 # ifndef PETSC_USE_COMPLEX
254  // get converged eigenpair
255  const PetscErrorCode ierr = EPSGetEigenpair(eps,
256  index,
257  &real_eigenvalues,
258  &imag_eigenvalues,
259  real_eigenvectors,
260  imag_eigenvectors);
261  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
262 # else
263  Assert(
264  (false),
265  ExcMessage(
266  "Your PETSc/SLEPc installation was configured with scalar-type complex "
267  "but this function is not defined for complex types."));
268 
269  // Cast to void to silence compiler warnings
270  (void)index;
271  (void)real_eigenvalues;
272  (void)imag_eigenvalues;
273  (void)real_eigenvectors;
274  (void)imag_eigenvectors;
275 # endif
276  }
277 
278  void
280  {
281  switch (state)
282  {
283  case ::SolverControl::iterate:
284  reason = EPS_CONVERGED_ITERATING;
285  break;
286 
287  case ::SolverControl::success:
288  reason = static_cast<EPSConvergedReason>(1);
289  break;
290 
291  case ::SolverControl::failure:
293  reason = EPS_DIVERGED_ITS;
294  else
295  reason = EPS_DIVERGED_BREAKDOWN;
296  break;
297 
298  default:
299  Assert(false, ExcNotImplemented());
300  }
301  }
302 
303  /* ---------------------- SolverControls ----------------------- */
304  SolverControl &
306  {
307  return solver_control;
308  }
309 
310  int
312  EPS /*eps */,
313  PetscScalar /*real_eigenvalue */,
314  PetscScalar /*imag_eigenvalue */,
315  PetscReal /*residual norm associated to the eigenpair */,
316  PetscReal * /*(output) computed error estimate */,
317  void * /*solver_control_x*/)
318  {
319  // If the error estimate returned by the convergence test function is less
320  // than the tolerance, then the eigenvalue is accepted as converged.
321  // This function is undefined (future reference only).
322 
323  // return without failure.
324  return 0;
325  }
326 
327  /* ---------------------- SolverKrylovSchur ------------------------ */
329  const MPI_Comm & mpi_communicator,
330  const AdditionalData &data)
331  : SolverBase(cn, mpi_communicator)
332  , additional_data(data)
333  {
334  const PetscErrorCode ierr =
335  EPSSetType(eps, const_cast<char *>(EPSKRYLOVSCHUR));
336  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
337  }
338 
339  /* ---------------------- SolverArnoldi ------------------------ */
341  const bool delayed_reorthogonalization)
342  : delayed_reorthogonalization(delayed_reorthogonalization)
343  {}
344 
346  const MPI_Comm & mpi_communicator,
347  const AdditionalData &data)
348  : SolverBase(cn, mpi_communicator)
349  , additional_data(data)
350  {
351  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSARNOLDI));
352  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
353 
354  // if requested, set delayed reorthogonalization in the Arnoldi
355  // iteration.
357  {
358  ierr = EPSArnoldiSetDelayed(eps, PETSC_TRUE);
359  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
360  }
361  }
362 
363 
364  /* ---------------------- Lanczos ------------------------ */
365  SolverLanczos::AdditionalData::AdditionalData(const EPSLanczosReorthogType r)
366  : reorthog(r)
367  {}
368 
370  const MPI_Comm & mpi_communicator,
371  const AdditionalData &data)
372  : SolverBase(cn, mpi_communicator)
373  , additional_data(data)
374  {
375  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSLANCZOS));
376  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
377 
378  ierr = EPSLanczosSetReorthog(eps, additional_data.reorthog);
379  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
380  }
381 
382  /* ----------------------- Power ------------------------- */
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 *>(EPSPOWER));
390  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
391  }
392 
393  /* ---------------- Generalized Davidson ----------------- */
395  bool double_expansion)
396  : double_expansion(double_expansion)
397  {}
398 
400  SolverControl & cn,
401  const MPI_Comm & mpi_communicator,
402  const AdditionalData &data)
403  : SolverBase(cn, mpi_communicator)
404  , additional_data(data)
405  {
406  PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSGD));
407  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
408 
410  {
411  ierr = EPSGDSetDoubleExpansion(eps, PETSC_TRUE);
412  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
413  }
414  }
415 
416  /* ------------------ Jacobi Davidson -------------------- */
418  const MPI_Comm &mpi_communicator,
419  const AdditionalData &data)
420  : SolverBase(cn, mpi_communicator)
421  , additional_data(data)
422  {
423  const PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSJD));
424  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
425  }
426 
427  /* ---------------------- LAPACK ------------------------- */
429  const MPI_Comm & mpi_communicator,
430  const AdditionalData &data)
431  : SolverBase(cn, mpi_communicator)
432  , additional_data(data)
433  {
434  // 'Tis overwhelmingly likely that PETSc/SLEPc *always* has
435  // BLAS/LAPACK, but let's be defensive.
436 # if PETSC_HAVE_BLASLAPACK
437  const PetscErrorCode ierr = EPSSetType(eps, const_cast<char *>(EPSLAPACK));
438  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
439 # else
440  Assert(
441  (false),
442  ExcMessage(
443  "Your PETSc/SLEPc installation was not configured with BLAS/LAPACK "
444  "but this is needed to use the LAPACK solver."));
445 # endif
446  }
447 } // namespace SLEPcWrappers
448 
450 
451 #endif // DEAL_II_WITH_SLEPC
void set_problem_type(EPSProblemType set_problem)
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1474
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:1521
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:296
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:1411
SolverGeneralizedDavidson(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
EPSConvergedReason reason
Definition: slepc_solver.h:354
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:488
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
static const char A
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:595
static ::ExceptionBase & ExcNotImplemented()
AdditionalData(const EPSLanczosReorthogType r=EPS_LANCZOS_REORTHOG_FULL)
const AdditionalData additional_data
Definition: slepc_solver.h:444
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:641
Eigenvalue vector is filled.
SolverControl & solver_control
Definition: slepc_solver.h:291
T max(const T &t, const MPI_Comm &mpi_communicator)
SolverArnoldi(SolverControl &cn, const MPI_Comm &mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())