Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16:50: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 - 2021 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 
408  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
409  const AdditionalData &data = AdditionalData());
410 
411  protected:
416  };
417 
418 
419 
432  class SolverArnoldi : public SolverBase
433  {
434  public:
440  {
445  AdditionalData(const bool delayed_reorthogonalization = false);
446 
451  };
452 
459  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
460  const AdditionalData &data = AdditionalData());
461 
462  protected:
467  };
468 
469 
470 
483  class SolverLanczos : public SolverBase
484  {
485  public:
491  {
495  EPSLanczosReorthogType reorthog;
496 
502  const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
503  };
504 
511  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
512  const AdditionalData &data = AdditionalData());
513 
514  protected:
519  };
520 
521 
522 
535  class SolverPower : public SolverBase
536  {
537  public:
543  {};
544 
551  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
552  const AdditionalData &data = AdditionalData());
553 
554  protected:
559  };
560 
561 
562 
576  {
577  public:
583  {
588 
592  AdditionalData(bool double_expansion = false);
593  };
594 
601  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
602  const AdditionalData &data = AdditionalData());
603 
604  protected:
609  };
610 
611 
612 
626  {
627  public:
633  {};
634 
641  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
642  const AdditionalData &data = AdditionalData());
643 
644  protected:
649  };
650 
651 
652 
665  class SolverLAPACK : public SolverBase
666  {
667  public:
673  {};
674 
681  const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
682  const AdditionalData &data = AdditionalData());
683 
684  protected:
689  };
690 
691 
692 
693  // --------------------------- inline and template functions -----------
698  // todo: The logic of these functions can be simplified without breaking
699  // backward compatibility...
700 
701  template <typename OutputVector>
702  void
704  std::vector<PetscScalar> & eigenvalues,
705  std::vector<OutputVector> & eigenvectors,
706  const unsigned int n_eigenpairs)
707  {
708  // Panic if the number of eigenpairs wanted is out of bounds.
709  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
711 
712  // Set the matrices of the problem
713  set_matrices(A);
714 
715  // and solve
716  unsigned int n_converged = 0;
717  solve(n_eigenpairs, &n_converged);
718 
719  if (n_converged > n_eigenpairs)
720  n_converged = n_eigenpairs;
721  AssertThrow(n_converged == n_eigenpairs,
723  n_eigenpairs));
724 
726  eigenvectors.resize(n_converged, eigenvectors.front());
727  eigenvalues.resize(n_converged);
728 
729  for (unsigned int index = 0; index < n_converged; ++index)
730  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
731  }
732 
733  template <typename OutputVector>
734  void
736  const PETScWrappers::MatrixBase &B,
737  std::vector<PetscScalar> & eigenvalues,
738  std::vector<OutputVector> & eigenvectors,
739  const unsigned int n_eigenpairs)
740  {
741  // Guard against incompatible matrix sizes:
742  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
743  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
744 
745  // Panic if the number of eigenpairs wanted is out of bounds.
746  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
748 
749  // Set the matrices of the problem
750  set_matrices(A, B);
751 
752  // and solve
753  unsigned int n_converged = 0;
754  solve(n_eigenpairs, &n_converged);
755 
756  if (n_converged >= n_eigenpairs)
757  n_converged = n_eigenpairs;
758 
759  AssertThrow(n_converged == n_eigenpairs,
761  n_eigenpairs));
763 
764  eigenvectors.resize(n_converged, eigenvectors.front());
765  eigenvalues.resize(n_converged);
766 
767  for (unsigned int index = 0; index < n_converged; ++index)
768  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
769  }
770 
771  template <typename OutputVector>
772  void
774  const PETScWrappers::MatrixBase &B,
775  std::vector<double> & real_eigenvalues,
776  std::vector<double> & imag_eigenvalues,
777  std::vector<OutputVector> & real_eigenvectors,
778  std::vector<OutputVector> & imag_eigenvectors,
779  const unsigned int n_eigenpairs)
780  {
781  // Guard against incompatible matrix sizes:
782  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
783  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
784 
785  // and incompatible eigenvalue/eigenvector sizes
786  AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
787  ExcDimensionMismatch(real_eigenvalues.size(),
788  imag_eigenvalues.size()));
789  AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
790  ExcDimensionMismatch(real_eigenvectors.size(),
791  imag_eigenvectors.size()));
792 
793  // Panic if the number of eigenpairs wanted is out of bounds.
794  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
796 
797  // Set the matrices of the problem
798  set_matrices(A, B);
799 
800  // and solve
801  unsigned int n_converged = 0;
802  solve(n_eigenpairs, &n_converged);
803 
804  if (n_converged >= n_eigenpairs)
805  n_converged = n_eigenpairs;
806 
807  AssertThrow(n_converged == n_eigenpairs,
809  n_eigenpairs));
810  AssertThrow((real_eigenvectors.size() != 0) &&
811  (imag_eigenvectors.size() != 0),
813 
814  real_eigenvectors.resize(n_converged, real_eigenvectors.front());
815  imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
816  real_eigenvalues.resize(n_converged);
817  imag_eigenvalues.resize(n_converged);
818 
819  for (unsigned int index = 0; index < n_converged; ++index)
820  get_eigenpair(index,
821  real_eigenvalues[index],
822  imag_eigenvalues[index],
823  real_eigenvectors[index],
824  imag_eigenvectors[index]);
825  }
826 
827  template <typename Vector>
828  void
829  SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
830  {
831  std::vector<Vec> vecs(this_initial_space.size());
832 
833  for (unsigned int i = 0; i < this_initial_space.size(); ++i)
834  {
835  Assert(this_initial_space[i].l2_norm() > 0.0,
836  ExcMessage("Initial vectors should be nonzero."));
837  vecs[i] = this_initial_space[i];
838  }
839 
840  // if the eigensolver supports only a single initial vector, but several
841  // guesses are provided, then all except the first one will be discarded.
842  // One could still build a vector that is rich in the directions of all
843  // guesses, by taking a linear combination of them. (TODO: make function
844  // virtual?)
845 
846  const PetscErrorCode ierr =
847  EPSSetInitialSpace(eps, vecs.size(), vecs.data());
848  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
849  }
850 
851 } // namespace SLEPcWrappers
852 
854 
855 # endif // DEAL_II_WITH_SLEPC
856 
857 /*---------------------------- slepc_solver.h ---------------------------*/
858 
859 #endif
860 
861 /*---------------------------- 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:466
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:829
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:703
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:648
const AdditionalData additional_data
Definition: slepc_solver.h:415
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:688
const AdditionalData additional_data
Definition: slepc_solver.h:518
SolverLanczos(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
Definition: slepc_solver.h:558
SolverPower(SolverControl &cn, const MPI_Comm mpi_communicator=PETSC_COMM_SELF, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcSLEPcError(int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:533
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
static ::ExceptionBase & ExcSLEPcWrappersUsageError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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)