Reference documentation for deal.II version 9.2.0
\(\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 
137 namespace SLEPcWrappers
138 {
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,
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 
387  {
388  public:
394  {};
395 
402  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
403  const AdditionalData &data = AdditionalData());
404 
405  protected:
410  };
411 
420  class SolverArnoldi : public SolverBase
421  {
422  public:
428  {
433  AdditionalData(const bool delayed_reorthogonalization = false);
434 
439  };
440 
447  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
448  const AdditionalData &data = AdditionalData());
449 
450  protected:
455  };
456 
465  class SolverLanczos : public SolverBase
466  {
467  public:
473  {
477  EPSLanczosReorthogType reorthog;
478 
484  const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
485  };
486 
493  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
494  const AdditionalData &data = AdditionalData());
495 
496  protected:
501  };
502 
511  class SolverPower : public SolverBase
512  {
513  public:
519  {};
520 
527  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
528  const AdditionalData &data = AdditionalData());
529 
530  protected:
535  };
536 
546  {
547  public:
553  {
558 
562  AdditionalData(bool double_expansion = false);
563  };
564 
571  SolverControl & cn,
572  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
573  const AdditionalData &data = AdditionalData());
574 
575  protected:
580  };
581 
591  {
592  public:
598  {};
599 
606  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
607  const AdditionalData &data = AdditionalData());
608 
609  protected:
614  };
615 
616 
625  class SolverLAPACK : public SolverBase
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  // --------------------------- inline and template functions -----------
656  // todo: The logic of these functions can be simplified without breaking
657  // backward compatibility...
658 
659  template <typename OutputVector>
660  void
662  std::vector<PetscScalar> & eigenvalues,
663  std::vector<OutputVector> & eigenvectors,
664  const unsigned int n_eigenpairs)
665  {
666  // Panic if the number of eigenpairs wanted is out of bounds.
667  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
669 
670  // Set the matrices of the problem
671  set_matrices(A);
672 
673  // and solve
674  unsigned int n_converged = 0;
675  solve(n_eigenpairs, &n_converged);
676 
677  if (n_converged > n_eigenpairs)
678  n_converged = n_eigenpairs;
679  AssertThrow(n_converged == n_eigenpairs,
681  n_eigenpairs));
682 
683  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
684  eigenvectors.resize(n_converged, eigenvectors.front());
685  eigenvalues.resize(n_converged);
686 
687  for (unsigned int index = 0; index < n_converged; ++index)
688  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
689  }
690 
691  template <typename OutputVector>
692  void
694  const PETScWrappers::MatrixBase &B,
695  std::vector<PetscScalar> & eigenvalues,
696  std::vector<OutputVector> & eigenvectors,
697  const unsigned int n_eigenpairs)
698  {
699  // Guard against incompatible matrix sizes:
700  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
701  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
702 
703  // Panic if the number of eigenpairs wanted is out of bounds.
704  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
706 
707  // Set the matrices of the problem
708  set_matrices(A, B);
709 
710  // and solve
711  unsigned int n_converged = 0;
712  solve(n_eigenpairs, &n_converged);
713 
714  if (n_converged >= n_eigenpairs)
715  n_converged = n_eigenpairs;
716 
717  AssertThrow(n_converged == n_eigenpairs,
719  n_eigenpairs));
720  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
721 
722  eigenvectors.resize(n_converged, eigenvectors.front());
723  eigenvalues.resize(n_converged);
724 
725  for (unsigned int index = 0; index < n_converged; ++index)
726  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
727  }
728 
729  template <typename OutputVector>
730  void
732  const PETScWrappers::MatrixBase &B,
733  std::vector<double> & real_eigenvalues,
734  std::vector<double> & imag_eigenvalues,
735  std::vector<OutputVector> & real_eigenvectors,
736  std::vector<OutputVector> & imag_eigenvectors,
737  const unsigned int n_eigenpairs)
738  {
739  // Guard against incompatible matrix sizes:
740  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
741  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
742 
743  // and incompatible eigenvalue/eigenvector sizes
744  AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
745  ExcDimensionMismatch(real_eigenvalues.size(),
746  imag_eigenvalues.size()));
747  AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
748  ExcDimensionMismatch(real_eigenvectors.size(),
749  imag_eigenvectors.size()));
750 
751  // Panic if the number of eigenpairs wanted is out of bounds.
752  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
754 
755  // Set the matrices of the problem
756  set_matrices(A, B);
757 
758  // and solve
759  unsigned int n_converged = 0;
760  solve(n_eigenpairs, &n_converged);
761 
762  if (n_converged >= n_eigenpairs)
763  n_converged = n_eigenpairs;
764 
765  AssertThrow(n_converged == n_eigenpairs,
767  n_eigenpairs));
768  AssertThrow((real_eigenvectors.size() != 0) &&
769  (imag_eigenvectors.size() != 0),
771 
772  real_eigenvectors.resize(n_converged, real_eigenvectors.front());
773  imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
774  real_eigenvalues.resize(n_converged);
775  imag_eigenvalues.resize(n_converged);
776 
777  for (unsigned int index = 0; index < n_converged; ++index)
778  get_eigenpair(index,
779  real_eigenvalues[index],
780  imag_eigenvalues[index],
781  real_eigenvectors[index],
782  imag_eigenvectors[index]);
783  }
784 
785  template <typename Vector>
786  void
787  SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
788  {
789  std::vector<Vec> vecs(this_initial_space.size());
790 
791  for (unsigned int i = 0; i < this_initial_space.size(); i++)
792  {
793  Assert(this_initial_space[i].l2_norm() > 0.0,
794  ExcMessage("Initial vectors should be nonzero."));
795  vecs[i] = this_initial_space[i];
796  }
797 
798  // if the eigensolver supports only a single initial vector, but several
799  // guesses are provided, then all except the first one will be discarded.
800  // One could still build a vector that is rich in the directions of all
801  // guesses, by taking a linear combination of them. (TODO: make function
802  // virtual?)
803 
804  const PetscErrorCode ierr =
805  EPSSetInitialSpace(eps, vecs.size(), vecs.data());
806  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
807  }
808 
809 } // namespace SLEPcWrappers
810 
812 
813 # endif // DEAL_II_WITH_SLEPC
814 
815 /*---------------------------- slepc_solver.h ---------------------------*/
816 
817 #endif
818 
819 /*---------------------------- slepc_solver.h ---------------------------*/
static ::ExceptionBase & ExcSLEPcWrappersUsageError()
void set_problem_type(EPSProblemType set_problem)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
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:534
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:1531
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:518
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:1419
void set_initial_space(const std::vector< Vector > &initial_space)
Definition: slepc_solver.h:787
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const AdditionalData additional_data
Definition: slepc_solver.h:648
#define DeclException0(Exception0)
Definition: exceptions.h:473
EPSConvergedReason reason
Definition: slepc_solver.h:360
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:500
const AdditionalData additional_data
Definition: slepc_solver.h:409
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:358
SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
Definition: slepc_solver.cc:35
const AdditionalData additional_data
Definition: slepc_solver.h:613
const AdditionalData additional_data
Definition: slepc_solver.h:454
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:661
Eigenvalue vector is filled.
SolverControl & solver_control
Definition: slepc_solver.h:297