Reference documentation for deal.II version Git 713825e468 2021-05-17 16:05:53 -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 {
147  {
148  public:
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 
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,
321  PETScWrappers::VectorBase &eigenvectors);
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  SolverControl & cn,
602  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
603  const AdditionalData &data = AdditionalData());
604 
605  protected:
610  };
611 
612 
613 
627  {
628  public:
634  {};
635 
642  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
643  const AdditionalData &data = AdditionalData());
644 
645  protected:
650  };
651 
652 
653 
666  class SolverLAPACK : public SolverBase
667  {
668  public:
674  {};
675 
682  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
683  const AdditionalData &data = AdditionalData());
684 
685  protected:
690  };
691 
692 
693 
694  // --------------------------- inline and template functions -----------
699  // todo: The logic of these functions can be simplified without breaking
700  // backward compatibility...
701 
702  template <typename OutputVector>
703  void
705  std::vector<PetscScalar> & eigenvalues,
706  std::vector<OutputVector> & eigenvectors,
707  const unsigned int n_eigenpairs)
708  {
709  // Panic if the number of eigenpairs wanted is out of bounds.
710  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
712 
713  // Set the matrices of the problem
714  set_matrices(A);
715 
716  // and solve
717  unsigned int n_converged = 0;
718  solve(n_eigenpairs, &n_converged);
719 
720  if (n_converged > n_eigenpairs)
721  n_converged = n_eigenpairs;
722  AssertThrow(n_converged == n_eigenpairs,
724  n_eigenpairs));
725 
726  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
727  eigenvectors.resize(n_converged, eigenvectors.front());
728  eigenvalues.resize(n_converged);
729 
730  for (unsigned int index = 0; index < n_converged; ++index)
731  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
732  }
733 
734  template <typename OutputVector>
735  void
737  const PETScWrappers::MatrixBase &B,
738  std::vector<PetscScalar> & eigenvalues,
739  std::vector<OutputVector> & eigenvectors,
740  const unsigned int n_eigenpairs)
741  {
742  // Guard against incompatible matrix sizes:
743  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
744  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
745 
746  // Panic if the number of eigenpairs wanted is out of bounds.
747  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
749 
750  // Set the matrices of the problem
751  set_matrices(A, B);
752 
753  // and solve
754  unsigned int n_converged = 0;
755  solve(n_eigenpairs, &n_converged);
756 
757  if (n_converged >= n_eigenpairs)
758  n_converged = n_eigenpairs;
759 
760  AssertThrow(n_converged == n_eigenpairs,
762  n_eigenpairs));
763  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
764 
765  eigenvectors.resize(n_converged, eigenvectors.front());
766  eigenvalues.resize(n_converged);
767 
768  for (unsigned int index = 0; index < n_converged; ++index)
769  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
770  }
771 
772  template <typename OutputVector>
773  void
775  const PETScWrappers::MatrixBase &B,
776  std::vector<double> & real_eigenvalues,
777  std::vector<double> & imag_eigenvalues,
778  std::vector<OutputVector> & real_eigenvectors,
779  std::vector<OutputVector> & imag_eigenvectors,
780  const unsigned int n_eigenpairs)
781  {
782  // Guard against incompatible matrix sizes:
783  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
784  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
785 
786  // and incompatible eigenvalue/eigenvector sizes
787  AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
788  ExcDimensionMismatch(real_eigenvalues.size(),
789  imag_eigenvalues.size()));
790  AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
791  ExcDimensionMismatch(real_eigenvectors.size(),
792  imag_eigenvectors.size()));
793 
794  // Panic if the number of eigenpairs wanted is out of bounds.
795  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
797 
798  // Set the matrices of the problem
799  set_matrices(A, B);
800 
801  // and solve
802  unsigned int n_converged = 0;
803  solve(n_eigenpairs, &n_converged);
804 
805  if (n_converged >= n_eigenpairs)
806  n_converged = n_eigenpairs;
807 
808  AssertThrow(n_converged == n_eigenpairs,
810  n_eigenpairs));
811  AssertThrow((real_eigenvectors.size() != 0) &&
812  (imag_eigenvectors.size() != 0),
814 
815  real_eigenvectors.resize(n_converged, real_eigenvectors.front());
816  imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
817  real_eigenvalues.resize(n_converged);
818  imag_eigenvalues.resize(n_converged);
819 
820  for (unsigned int index = 0; index < n_converged; ++index)
821  get_eigenpair(index,
822  real_eigenvalues[index],
823  imag_eigenvalues[index],
824  real_eigenvectors[index],
825  imag_eigenvectors[index]);
826  }
827 
828  template <typename Vector>
829  void
830  SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
831  {
832  std::vector<Vec> vecs(this_initial_space.size());
833 
834  for (unsigned int i = 0; i < this_initial_space.size(); i++)
835  {
836  Assert(this_initial_space[i].l2_norm() > 0.0,
837  ExcMessage("Initial vectors should be nonzero."));
838  vecs[i] = this_initial_space[i];
839  }
840 
841  // if the eigensolver supports only a single initial vector, but several
842  // guesses are provided, then all except the first one will be discarded.
843  // One could still build a vector that is rich in the directions of all
844  // guesses, by taking a linear combination of them. (TODO: make function
845  // virtual?)
846 
847  const PetscErrorCode ierr =
848  EPSSetInitialSpace(eps, vecs.size(), vecs.data());
849  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
850  }
851 
852 } // namespace SLEPcWrappers
853 
855 
856 # endif // DEAL_II_WITH_SLEPC
857 
858 /*---------------------------- slepc_solver.h ---------------------------*/
859 
860 #endif
861 
862 /*---------------------------- 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:558
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:1575
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:302
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:77
#define Assert(cond, exc)
Definition: exceptions.h:1465
void set_initial_space(const std::vector< Vector > &initial_space)
Definition: slepc_solver.h:830
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const AdditionalData additional_data
Definition: slepc_solver.h:689
#define DeclException0(Exception0)
Definition: exceptions.h:470
EPSConvergedReason reason
Definition: slepc_solver.h:360
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:518
const AdditionalData additional_data
Definition: slepc_solver.h:415
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
Definition: slepc_solver.cc:98
static const char A
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
Definition: slepc_solver.cc:35
const AdditionalData additional_data
Definition: slepc_solver.h:649
const AdditionalData additional_data
Definition: slepc_solver.h:466
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:704
Eigenvalue vector is filled.
SolverControl & solver_control
Definition: slepc_solver.h:297