Reference documentation for deal.II version Git 0096380e24 2020-08-06 21:03:09 -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.h
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 
16 
17 #ifndef dealii_slepc_solver_h
18 # define dealii_slepc_solver_h
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_SLEPC
23 
27 
28 # include <petscconf.h>
29 # include <petscksp.h>
30 
31 # include <slepceps.h>
32 
33 # include <memory>
34 
36 
131 namespace SLEPcWrappers
132 {
141  {
142  public:
147  SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator);
148 
152  virtual ~SolverBase();
153 
172  template <typename OutputVector>
173  void
175  std::vector<PetscScalar> & eigenvalues,
176  std::vector<OutputVector> & eigenvectors,
177  const unsigned int n_eigenpairs = 1);
178 
184  template <typename OutputVector>
185  void
187  const PETScWrappers::MatrixBase &B,
188  std::vector<PetscScalar> & eigenvalues,
189  std::vector<OutputVector> & eigenvectors,
190  const unsigned int n_eigenpairs = 1);
191 
197  template <typename OutputVector>
198  void
200  const PETScWrappers::MatrixBase &B,
201  std::vector<double> & real_eigenvalues,
202  std::vector<double> & imag_eigenvalues,
203  std::vector<OutputVector> & real_eigenvectors,
204  std::vector<OutputVector> & imag_eigenvectors,
205  const unsigned int n_eigenpairs = 1);
206 
213  template <typename Vector>
214  void
215  set_initial_space(const std::vector<Vector> &initial_space);
216 
220  void
222 
227  void
228  set_target_eigenvalue(const PetscScalar &this_target);
229 
236  void
237  set_which_eigenpairs(EPSWhich set_which);
238 
247  void
248  set_problem_type(EPSProblemType set_problem);
249 
255  void
257 
262 
267  int,
268  << " An error with error number " << arg1
269  << " occurred while calling a SLEPc function");
270 
275  int,
276  int,
277  << " The number of converged eigenvectors is " << arg1
278  << " but " << arg2 << " were requested. ");
279 
283  SolverControl &
284  control() const;
285 
286  protected:
292 
296  const MPI_Comm mpi_communicator;
297 
304  void
305  solve(const unsigned int n_eigenpairs, unsigned int *n_converged);
306 
312  void
313  get_eigenpair(const unsigned int index,
314  PetscScalar & eigenvalues,
315  PETScWrappers::VectorBase &eigenvectors);
316 
322  void
323  get_eigenpair(const unsigned int index,
324  double & real_eigenvalues,
325  double & imag_eigenvalues,
326  PETScWrappers::VectorBase &real_eigenvectors,
327  PETScWrappers::VectorBase &imag_eigenvectors);
328 
333  void
335 
340  void
342  const PETScWrappers::MatrixBase &B);
343 
344  protected:
348  EPS eps;
349 
350  private:
354  EPSConvergedReason reason;
355 
356 
363  static int
364  convergence_test(EPS eps,
365  PetscScalar real_eigenvalue,
366  PetscScalar imag_eigenvalue,
367  PetscReal residual_norm,
368  PetscReal * estimated_error,
369  void * solver_control);
370  };
371 
379  {
380  public:
386  {};
387 
394  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
395  const AdditionalData &data = AdditionalData());
396 
397  protected:
402  };
403 
410  class SolverArnoldi : public SolverBase
411  {
412  public:
418  {
423  AdditionalData(const bool delayed_reorthogonalization = false);
424 
429  };
430 
437  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
438  const AdditionalData &data = AdditionalData());
439 
440  protected:
445  };
446 
453  class SolverLanczos : public SolverBase
454  {
455  public:
461  {
465  EPSLanczosReorthogType reorthog;
466 
472  const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
473  };
474 
481  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
482  const AdditionalData &data = AdditionalData());
483 
484  protected:
489  };
490 
497  class SolverPower : public SolverBase
498  {
499  public:
505  {};
506 
513  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
514  const AdditionalData &data = AdditionalData());
515 
516  protected:
521  };
522 
530  {
531  public:
537  {
542 
546  AdditionalData(bool double_expansion = false);
547  };
548 
555  SolverControl & cn,
556  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
557  const AdditionalData &data = AdditionalData());
558 
559  protected:
564  };
565 
573  {
574  public:
580  {};
581 
588  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
589  const AdditionalData &data = AdditionalData());
590 
591  protected:
596  };
597 
598 
605  class SolverLAPACK : public SolverBase
606  {
607  public:
613  {};
614 
621  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
622  const AdditionalData &data = AdditionalData());
623 
624  protected:
629  };
630 
631  // --------------------------- inline and template functions -----------
636  // todo: The logic of these functions can be simplified without breaking
637  // backward compatibility...
638 
639  template <typename OutputVector>
640  void
642  std::vector<PetscScalar> & eigenvalues,
643  std::vector<OutputVector> & eigenvectors,
644  const unsigned int n_eigenpairs)
645  {
646  // Panic if the number of eigenpairs wanted is out of bounds.
647  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
649 
650  // Set the matrices of the problem
651  set_matrices(A);
652 
653  // and solve
654  unsigned int n_converged = 0;
655  solve(n_eigenpairs, &n_converged);
656 
657  if (n_converged > n_eigenpairs)
658  n_converged = n_eigenpairs;
659  AssertThrow(n_converged == n_eigenpairs,
661  n_eigenpairs));
662 
663  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
664  eigenvectors.resize(n_converged, eigenvectors.front());
665  eigenvalues.resize(n_converged);
666 
667  for (unsigned int index = 0; index < n_converged; ++index)
668  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
669  }
670 
671  template <typename OutputVector>
672  void
674  const PETScWrappers::MatrixBase &B,
675  std::vector<PetscScalar> & eigenvalues,
676  std::vector<OutputVector> & eigenvectors,
677  const unsigned int n_eigenpairs)
678  {
679  // Guard against incompatible matrix sizes:
680  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
681  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
682 
683  // Panic if the number of eigenpairs wanted is out of bounds.
684  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
686 
687  // Set the matrices of the problem
688  set_matrices(A, B);
689 
690  // and solve
691  unsigned int n_converged = 0;
692  solve(n_eigenpairs, &n_converged);
693 
694  if (n_converged >= n_eigenpairs)
695  n_converged = n_eigenpairs;
696 
697  AssertThrow(n_converged == n_eigenpairs,
699  n_eigenpairs));
700  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
701 
702  eigenvectors.resize(n_converged, eigenvectors.front());
703  eigenvalues.resize(n_converged);
704 
705  for (unsigned int index = 0; index < n_converged; ++index)
706  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
707  }
708 
709  template <typename OutputVector>
710  void
712  const PETScWrappers::MatrixBase &B,
713  std::vector<double> & real_eigenvalues,
714  std::vector<double> & imag_eigenvalues,
715  std::vector<OutputVector> & real_eigenvectors,
716  std::vector<OutputVector> & imag_eigenvectors,
717  const unsigned int n_eigenpairs)
718  {
719  // Guard against incompatible matrix sizes:
720  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
721  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
722 
723  // and incompatible eigenvalue/eigenvector sizes
724  AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
725  ExcDimensionMismatch(real_eigenvalues.size(),
726  imag_eigenvalues.size()));
727  AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
728  ExcDimensionMismatch(real_eigenvectors.size(),
729  imag_eigenvectors.size()));
730 
731  // Panic if the number of eigenpairs wanted is out of bounds.
732  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
734 
735  // Set the matrices of the problem
736  set_matrices(A, B);
737 
738  // and solve
739  unsigned int n_converged = 0;
740  solve(n_eigenpairs, &n_converged);
741 
742  if (n_converged >= n_eigenpairs)
743  n_converged = n_eigenpairs;
744 
745  AssertThrow(n_converged == n_eigenpairs,
747  n_eigenpairs));
748  AssertThrow((real_eigenvectors.size() != 0) &&
749  (imag_eigenvectors.size() != 0),
751 
752  real_eigenvectors.resize(n_converged, real_eigenvectors.front());
753  imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
754  real_eigenvalues.resize(n_converged);
755  imag_eigenvalues.resize(n_converged);
756 
757  for (unsigned int index = 0; index < n_converged; ++index)
758  get_eigenpair(index,
759  real_eigenvalues[index],
760  imag_eigenvalues[index],
761  real_eigenvectors[index],
762  imag_eigenvectors[index]);
763  }
764 
765  template <typename Vector>
766  void
767  SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
768  {
769  std::vector<Vec> vecs(this_initial_space.size());
770 
771  for (unsigned int i = 0; i < this_initial_space.size(); i++)
772  {
773  Assert(this_initial_space[i].l2_norm() > 0.0,
774  ExcMessage("Initial vectors should be nonzero."));
775  vecs[i] = this_initial_space[i];
776  }
777 
778  // if the eigensolver supports only a single initial vector, but several
779  // guesses are provided, then all except the first one will be discarded.
780  // One could still build a vector that is rich in the directions of all
781  // guesses, by taking a linear combination of them. (TODO: make function
782  // virtual?)
783 
784  const PetscErrorCode ierr =
785  EPSSetInitialSpace(eps, vecs.size(), vecs.data());
786  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
787  }
788 
789 } // namespace SLEPcWrappers
790 
792 
793 # endif // DEAL_II_WITH_SLEPC
794 
795 /*---------------------------- slepc_solver.h ---------------------------*/
796 
797 #endif
798 
799 /*---------------------------- slepc_solver.h ---------------------------*/
static ::ExceptionBase & ExcSLEPcWrappersUsageError()
void set_problem_type(EPSProblemType set_problem)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
static int convergence_test(EPS eps, PetscScalar real_eigenvalue, PetscScalar imag_eigenvalue, PetscReal residual_norm, PetscReal *estimated_error, void *solver_control)
const AdditionalData additional_data
Definition: slepc_solver.h:520
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)
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
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)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
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
void set_initial_space(const std::vector< Vector > &initial_space)
Definition: slepc_solver.h:767
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const AdditionalData additional_data
Definition: slepc_solver.h:628
#define DeclException0(Exception0)
Definition: exceptions.h:470
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
const AdditionalData additional_data
Definition: slepc_solver.h:401
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
Definition: slepc_solver.cc:90
static const char A
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
Definition: slepc_solver.cc:35
const AdditionalData additional_data
Definition: slepc_solver.h:595
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