Reference documentation for deal.II version GIT c14369f203 2023-10-01 07:40:02+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\}}\)
slepc_solver.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2023 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 {
147  {
148  public:
153  SolverBase(SolverControl &cn, const MPI_Comm mpi_communicator);
154 
158  virtual ~SolverBase();
159 
178  template <typename OutputVector>
179  void
181  std::vector<PetscScalar> &eigenvalues,
182  std::vector<OutputVector> &eigenvectors,
183  const unsigned int n_eigenpairs = 1);
184 
190  template <typename OutputVector>
191  void
193  const PETScWrappers::MatrixBase &B,
194  std::vector<PetscScalar> &eigenvalues,
195  std::vector<OutputVector> &eigenvectors,
196  const unsigned int n_eigenpairs = 1);
197 
203  template <typename OutputVector>
204  void
206  const PETScWrappers::MatrixBase &B,
207  std::vector<double> &real_eigenvalues,
208  std::vector<double> &imag_eigenvalues,
209  std::vector<OutputVector> &real_eigenvectors,
210  std::vector<OutputVector> &imag_eigenvectors,
211  const unsigned int n_eigenpairs = 1);
212 
219  template <typename Vector>
220  void
221  set_initial_space(const std::vector<Vector> &initial_space);
222 
226  void
228 
233  void
234  set_target_eigenvalue(const PetscScalar &this_target);
235 
242  void
243  set_which_eigenpairs(EPSWhich set_which);
244 
253  void
254  set_problem_type(EPSProblemType set_problem);
255 
261  void
263 
268 
273  int,
274  << " An error with error number " << arg1
275  << " occurred while calling a SLEPc function");
276 
281  int,
282  int,
283  << " The number of converged eigenvectors is " << arg1
284  << " but " << arg2 << " were requested. ");
285 
289  SolverControl &
290  control() const;
291 
292  protected:
298 
302  const MPI_Comm mpi_communicator;
303 
310  void
311  solve(const unsigned int n_eigenpairs, unsigned int *n_converged);
312 
318  void
319  get_eigenpair(const unsigned int index,
320  PetscScalar &eigenvalues,
322 
328  void
329  get_eigenpair(const unsigned int index,
330  double &real_eigenvalues,
331  double &imag_eigenvalues,
332  PETScWrappers::VectorBase &real_eigenvectors,
333  PETScWrappers::VectorBase &imag_eigenvectors);
334 
339  void
341 
346  void
348  const PETScWrappers::MatrixBase &B);
349 
350  protected:
354  EPS eps;
355 
356  private:
360  EPSConvergedReason reason;
361 
362 
369  static int
370  convergence_test(EPS eps,
371  PetscScalar real_eigenvalue,
372  PetscScalar imag_eigenvalue,
373  PetscReal residual_norm,
374  PetscReal *estimated_error,
375  void *solver_control);
376  };
377 
378 
379 
393  {
394  public:
400  {};
401 
407  explicit SolverKrylovSchur(
408  SolverControl &cn,
409  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
410  const AdditionalData &data = AdditionalData());
411 
412  protected:
417  };
418 
419 
420 
433  class SolverArnoldi : public SolverBase
434  {
435  public:
441  {
446  explicit AdditionalData(const bool delayed_reorthogonalization = false);
447 
452  };
453 
459  explicit SolverArnoldi(SolverControl &cn,
460  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
461  const AdditionalData &data = AdditionalData());
462 
463  protected:
468  };
469 
470 
471 
484  class SolverLanczos : public SolverBase
485  {
486  public:
492  {
496  EPSLanczosReorthogType reorthog;
497 
502  explicit AdditionalData(
503  const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
504  };
505 
511  explicit SolverLanczos(SolverControl &cn,
512  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
513  const AdditionalData &data = AdditionalData());
514 
515  protected:
520  };
521 
522 
523 
536  class SolverPower : public SolverBase
537  {
538  public:
544  {};
545 
551  explicit SolverPower(SolverControl &cn,
552  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
553  const AdditionalData &data = AdditionalData());
554 
555  protected:
560  };
561 
562 
563 
577  {
578  public:
584  {
589 
593  explicit AdditionalData(bool double_expansion = false);
594  };
595 
601  explicit SolverGeneralizedDavidson(
602  SolverControl &cn,
603  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
604  const AdditionalData &data = AdditionalData());
605 
606  protected:
611  };
612 
613 
614 
628  {
629  public:
635  {};
636 
642  explicit SolverJacobiDavidson(
643  SolverControl &cn,
644  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
645  const AdditionalData &data = AdditionalData());
646 
647  protected:
652  };
653 
654 
655 
668  class SolverLAPACK : public SolverBase
669  {
670  public:
676  {};
677 
683  explicit SolverLAPACK(SolverControl &cn,
684  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
685  const AdditionalData &data = AdditionalData());
686 
687  protected:
692  };
693 
694 
695 
696  // --------------------------- inline and template functions -----------
701  // todo: The logic of these functions can be simplified without breaking
702  // backward compatibility...
703 
704  template <typename OutputVector>
705  void
707  std::vector<PetscScalar> &eigenvalues,
708  std::vector<OutputVector> &eigenvectors,
709  const unsigned int n_eigenpairs)
710  {
711  // Panic if the number of eigenpairs wanted is out of bounds.
712  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
714 
715  // Set the matrices of the problem
716  set_matrices(A);
717 
718  // and solve
719  unsigned int n_converged = 0;
720  solve(n_eigenpairs, &n_converged);
721 
722  if (n_converged > n_eigenpairs)
723  n_converged = n_eigenpairs;
724  AssertThrow(n_converged == n_eigenpairs,
726  n_eigenpairs));
727 
729  eigenvectors.resize(n_converged, eigenvectors.front());
730  eigenvalues.resize(n_converged);
731 
732  for (unsigned int index = 0; index < n_converged; ++index)
733  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
734  }
735 
736  template <typename OutputVector>
737  void
739  const PETScWrappers::MatrixBase &B,
740  std::vector<PetscScalar> &eigenvalues,
741  std::vector<OutputVector> &eigenvectors,
742  const unsigned int n_eigenpairs)
743  {
744  // Guard against incompatible matrix sizes:
745  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
746  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
747 
748  // Panic if the number of eigenpairs wanted is out of bounds.
749  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
751 
752  // Set the matrices of the problem
753  set_matrices(A, B);
754 
755  // and solve
756  unsigned int n_converged = 0;
757  solve(n_eigenpairs, &n_converged);
758 
759  if (n_converged >= n_eigenpairs)
760  n_converged = n_eigenpairs;
761 
762  AssertThrow(n_converged == n_eigenpairs,
764  n_eigenpairs));
766 
767  eigenvectors.resize(n_converged, eigenvectors.front());
768  eigenvalues.resize(n_converged);
769 
770  for (unsigned int index = 0; index < n_converged; ++index)
771  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
772  }
773 
774  template <typename OutputVector>
775  void
777  const PETScWrappers::MatrixBase &B,
778  std::vector<double> &real_eigenvalues,
779  std::vector<double> &imag_eigenvalues,
780  std::vector<OutputVector> &real_eigenvectors,
781  std::vector<OutputVector> &imag_eigenvectors,
782  const unsigned int n_eigenpairs)
783  {
784  // Guard against incompatible matrix sizes:
785  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
786  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
787 
788  // and incompatible eigenvalue/eigenvector sizes
789  AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
790  ExcDimensionMismatch(real_eigenvalues.size(),
791  imag_eigenvalues.size()));
792  AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
793  ExcDimensionMismatch(real_eigenvectors.size(),
794  imag_eigenvectors.size()));
795 
796  // Panic if the number of eigenpairs wanted is out of bounds.
797  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
799 
800  // Set the matrices of the problem
801  set_matrices(A, B);
802 
803  // and solve
804  unsigned int n_converged = 0;
805  solve(n_eigenpairs, &n_converged);
806 
807  if (n_converged >= n_eigenpairs)
808  n_converged = n_eigenpairs;
809 
810  AssertThrow(n_converged == n_eigenpairs,
812  n_eigenpairs));
813  AssertThrow((real_eigenvectors.size() != 0) &&
814  (imag_eigenvectors.size() != 0),
816 
817  real_eigenvectors.resize(n_converged, real_eigenvectors.front());
818  imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
819  real_eigenvalues.resize(n_converged);
820  imag_eigenvalues.resize(n_converged);
821 
822  for (unsigned int index = 0; index < n_converged; ++index)
823  get_eigenpair(index,
824  real_eigenvalues[index],
825  imag_eigenvalues[index],
826  real_eigenvectors[index],
827  imag_eigenvectors[index]);
828  }
829 
830  template <typename Vector>
831  void
832  SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
833  {
834  std::vector<Vec> vecs(this_initial_space.size());
835 
836  for (unsigned int i = 0; i < this_initial_space.size(); ++i)
837  {
838  Assert(this_initial_space[i].l2_norm() > 0.0,
839  ExcMessage("Initial vectors should be nonzero."));
840  vecs[i] = this_initial_space[i];
841  }
842 
843  // if the eigensolver supports only a single initial vector, but several
844  // guesses are provided, then all except the first one will be discarded.
845  // One could still build a vector that is rich in the directions of all
846  // guesses, by taking a linear combination of them. (TODO: make function
847  // virtual?)
848 
849  const PetscErrorCode ierr =
850  EPSSetInitialSpace(eps, vecs.size(), vecs.data());
851  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
852  }
853 
854 } // namespace SLEPcWrappers
855 
857 
858 # endif // DEAL_II_WITH_SLEPC
859 
860 /*---------------------------- slepc_solver.h ---------------------------*/
861 
862 #endif
863 
864 /*---------------------------- slepc_solver.h ---------------------------*/
SolverArnoldi(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
Definition: slepc_solver.h:467
SolverControl & control() const
EPSConvergedReason reason
Definition: slepc_solver.h:360
void set_target_eigenvalue(const PetscScalar &this_target)
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:302
void set_matrices(const PETScWrappers::MatrixBase &A)
Definition: slepc_solver.cc:77
SolverBase(SolverControl &cn, const MPI_Comm mpi_communicator)
Definition: slepc_solver.cc:35
void get_solver_state(const SolverControl::State state)
void set_initial_space(const std::vector< Vector > &initial_space)
Definition: slepc_solver.h:832
void set_problem_type(EPSProblemType set_problem)
SolverControl & solver_control
Definition: slepc_solver.h:297
void solve(const PETScWrappers::MatrixBase &A, std::vector< PetscScalar > &eigenvalues, std::vector< OutputVector > &eigenvectors, const unsigned int n_eigenpairs=1)
Definition: slepc_solver.h:706
void get_eigenpair(const unsigned int index, PetscScalar &eigenvalues, PETScWrappers::VectorBase &eigenvectors)
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
Definition: slepc_solver.cc:98
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())
const AdditionalData additional_data
Definition: slepc_solver.h:651
const AdditionalData additional_data
Definition: slepc_solver.h:416
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
Definition: slepc_solver.h:691
const AdditionalData additional_data
Definition: slepc_solver.h:519
SolverLanczos(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
Definition: slepc_solver.h:559
SolverPower(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define DeclException0(Exception0)
Definition: exceptions.h:467
static ::ExceptionBase & ExcSLEPcError(int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:535
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcSLEPcWrappersUsageError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
@ eigenvalues
Eigenvalue vector is filled.
static const char A
AdditionalData(const bool delayed_reorthogonalization=false)
AdditionalData(const EPSLanczosReorthogType r=EPS_LANCZOS_REORTHOG_FULL)
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)