Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
slepc_solver.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2019 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 
24 # include <deal.II/lac/petsc_matrix_base.h>
25 # include <deal.II/lac/slepc_spectral_transformation.h>
26 # include <deal.II/lac/solver_control.h>
27 
28 # include <petscconf.h>
29 # include <petscksp.h>
30 
31 # include <slepceps.h>
32 
33 # include <memory>
34 
35 DEAL_II_NAMESPACE_OPEN
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 
216  DEAL_II_DEPRECATED
217  void
218  set_initial_vector(const PETScWrappers::VectorBase &this_initial_vector);
219 
226  template <typename Vector>
227  void
228  set_initial_space(const std::vector<Vector> &initial_space);
229 
233  void
235 
240  void
241  set_target_eigenvalue(const PetscScalar &this_target);
242 
249  void
250  set_which_eigenpairs(EPSWhich set_which);
251 
260  void
261  set_problem_type(EPSProblemType set_problem);
262 
268  void
270 
275 
280  int,
281  << " An error with error number " << arg1
282  << " occurred while calling a SLEPc function");
283 
288  int,
289  int,
290  << " The number of converged eigenvectors is " << arg1
291  << " but " << arg2 << " were requested. ");
292 
296  SolverControl &
297  control() const;
298 
299  protected:
305 
309  const MPI_Comm mpi_communicator;
310 
317  void
318  solve(const unsigned int n_eigenpairs, unsigned int *n_converged);
319 
325  void
326  get_eigenpair(const unsigned int index,
327  PetscScalar & eigenvalues,
328  PETScWrappers::VectorBase &eigenvectors);
329 
335  void
336  get_eigenpair(const unsigned int index,
337  double & real_eigenvalues,
338  double & imag_eigenvalues,
339  PETScWrappers::VectorBase &real_eigenvectors,
340  PETScWrappers::VectorBase &imag_eigenvectors);
341 
346  void
348 
353  void
355  const PETScWrappers::MatrixBase &B);
356 
357  protected:
361  EPS eps;
362 
363  private:
367  EPSConvergedReason reason;
368 
369 
376  static int
377  convergence_test(EPS eps,
378  PetscScalar real_eigenvalue,
379  PetscScalar imag_eigenvalue,
380  PetscReal residual_norm,
381  PetscReal * estimated_error,
382  void * solver_control);
383  };
384 
394  {
395  public:
401  {};
402 
409  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
410  const AdditionalData &data = AdditionalData());
411 
412  protected:
417  };
418 
427  class SolverArnoldi : public SolverBase
428  {
429  public:
435  {
440  AdditionalData(const bool delayed_reorthogonalization = false);
441 
446  };
447 
454  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
455  const AdditionalData &data = AdditionalData());
456 
457  protected:
462  };
463 
472  class SolverLanczos : public SolverBase
473  {
474  public:
480  {
484  EPSLanczosReorthogType reorthog;
485 
491  const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
492  };
493 
500  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
501  const AdditionalData &data = AdditionalData());
502 
503  protected:
508  };
509 
518  class SolverPower : public SolverBase
519  {
520  public:
526  {};
527 
534  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
535  const AdditionalData &data = AdditionalData());
536 
537  protected:
542  };
543 
553  {
554  public:
560  {
565 
569  AdditionalData(bool double_expansion = false);
570  };
571 
578  SolverControl & cn,
579  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
580  const AdditionalData &data = AdditionalData());
581 
582  protected:
587  };
588 
598  {
599  public:
605  {};
606 
613  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
614  const AdditionalData &data = AdditionalData());
615 
616  protected:
621  };
622 
623 
632  class SolverLAPACK : public SolverBase
633  {
634  public:
640  {};
641 
648  const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
649  const AdditionalData &data = AdditionalData());
650 
651  protected:
656  };
657 
658  // --------------------------- inline and template functions -----------
663  // todo: The logic of these functions can be simplified without breaking
664  // backward compatibility...
665 
666  template <typename OutputVector>
667  void
669  std::vector<PetscScalar> & eigenvalues,
670  std::vector<OutputVector> & eigenvectors,
671  const unsigned int n_eigenpairs)
672  {
673  // Panic if the number of eigenpairs wanted is out of bounds.
674  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
676 
677  // Set the matrices of the problem
678  set_matrices(A);
679 
680  // and solve
681  unsigned int n_converged = 0;
682  solve(n_eigenpairs, &n_converged);
683 
684  if (n_converged > n_eigenpairs)
685  n_converged = n_eigenpairs;
686  AssertThrow(n_converged == n_eigenpairs,
688  n_eigenpairs));
689 
690  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
691  eigenvectors.resize(n_converged, eigenvectors.front());
692  eigenvalues.resize(n_converged);
693 
694  for (unsigned int index = 0; index < n_converged; ++index)
695  get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
696  }
697 
698  template <typename OutputVector>
699  void
701  const PETScWrappers::MatrixBase &B,
702  std::vector<PetscScalar> & eigenvalues,
703  std::vector<OutputVector> & eigenvectors,
704  const unsigned int n_eigenpairs)
705  {
706  // Guard against incompatible matrix sizes:
707  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
708  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
709 
710  // Panic if the number of eigenpairs wanted is out of bounds.
711  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
713 
714  // Set the matrices of the problem
715  set_matrices(A, B);
716 
717  // and solve
718  unsigned int n_converged = 0;
719  solve(n_eigenpairs, &n_converged);
720 
721  if (n_converged >= n_eigenpairs)
722  n_converged = n_eigenpairs;
723 
724  AssertThrow(n_converged == n_eigenpairs,
726  n_eigenpairs));
727  AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
728 
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<double> & real_eigenvalues,
741  std::vector<double> & imag_eigenvalues,
742  std::vector<OutputVector> & real_eigenvectors,
743  std::vector<OutputVector> & imag_eigenvectors,
744  const unsigned int n_eigenpairs)
745  {
746  // Guard against incompatible matrix sizes:
747  AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
748  AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
749 
750  // and incompatible eigenvalue/eigenvector sizes
751  AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
752  ExcDimensionMismatch(real_eigenvalues.size(),
753  imag_eigenvalues.size()));
754  AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
755  ExcDimensionMismatch(real_eigenvectors.size(),
756  imag_eigenvectors.size()));
757 
758  // Panic if the number of eigenpairs wanted is out of bounds.
759  AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
761 
762  // Set the matrices of the problem
763  set_matrices(A, B);
764 
765  // and solve
766  unsigned int n_converged = 0;
767  solve(n_eigenpairs, &n_converged);
768 
769  if (n_converged >= n_eigenpairs)
770  n_converged = n_eigenpairs;
771 
772  AssertThrow(n_converged == n_eigenpairs,
774  n_eigenpairs));
775  AssertThrow((real_eigenvectors.size() != 0) &&
776  (imag_eigenvectors.size() != 0),
778 
779  real_eigenvectors.resize(n_converged, real_eigenvectors.front());
780  imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
781  real_eigenvalues.resize(n_converged);
782  imag_eigenvalues.resize(n_converged);
783 
784  for (unsigned int index = 0; index < n_converged; ++index)
785  get_eigenpair(index,
786  real_eigenvalues[index],
787  imag_eigenvalues[index],
788  real_eigenvectors[index],
789  imag_eigenvectors[index]);
790  }
791 
792  template <typename Vector>
793  void
794  SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
795  {
796  std::vector<Vec> vecs(this_initial_space.size());
797 
798  for (unsigned int i = 0; i < this_initial_space.size(); i++)
799  {
800  Assert(this_initial_space[i].l2_norm() > 0.0,
801  ExcMessage("Initial vectors should be nonzero."));
802  vecs[i] = this_initial_space[i];
803  }
804 
805  // if the eigensolver supports only a single initial vector, but several
806  // guesses are provided, then all except the first one will be discarded.
807  // One could still build a vector that is rich in the directions of all
808  // guesses, by taking a linear combination of them. (TODO: make function
809  // virtual?)
810 
811  const PetscErrorCode ierr =
812  EPSSetInitialSpace(eps, vecs.size(), vecs.data());
813  AssertThrow(ierr == 0, ExcSLEPcError(ierr));
814  }
815 
816 } // namespace SLEPcWrappers
817 
818 DEAL_II_NAMESPACE_CLOSE
819 
820 # endif // DEAL_II_WITH_SLEPC
821 
822 /*---------------------------- slepc_solver.h ---------------------------*/
823 
824 #endif
825 
826 /*---------------------------- 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:541
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::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
void set_which_eigenpairs(EPSWhich set_which)
void set_target_eigenvalue(const PetscScalar &this_target)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:309
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:1407
void set_initial_space(const std::vector< Vector > &initial_space)
Definition: slepc_solver.h:794
void set_initial_vector(const PETScWrappers::VectorBase &this_initial_vector)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const AdditionalData additional_data
Definition: slepc_solver.h:655
#define DeclException0(Exception0)
Definition: exceptions.h:473
EPSConvergedReason reason
Definition: slepc_solver.h:367
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:507
const AdditionalData additional_data
Definition: slepc_solver.h:416
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
Definition: slepc_solver.cc:90
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
Definition: slepc_solver.cc:35
const AdditionalData additional_data
Definition: slepc_solver.h:620
const AdditionalData additional_data
Definition: slepc_solver.h:461
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:668
SolverControl & solver_control
Definition: slepc_solver.h:304