Reference documentation for deal.II version 8.1.0
slepc_solver.h
1 // ---------------------------------------------------------------------
2 // @f$Id: slepc_solver.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2009 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #ifndef __deal2__slepc_solver_h
19 #define __deal2__slepc_solver_h
20 
21 #include <deal.II/base/config.h>
22 
23 #ifdef DEAL_II_WITH_SLEPC
24 
25 # include <deal.II/base/std_cxx1x/shared_ptr.h>
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/solver_control.h>
28 # include <deal.II/lac/petsc_matrix_base.h>
29 # include <deal.II/lac/slepc_spectral_transformation.h>
30 
31 # include <petscconf.h>
32 # include <petscksp.h>
33 # include <slepceps.h>
34 
36 
108 namespace SLEPcWrappers
109 {
110 
120  {
121  public:
127  const MPI_Comm &mpi_communicator);
128 
132  virtual ~SolverBase ();
133 
154  template <typename OutputVector>
155  void
157  std::vector<PetscScalar> &eigenvalues,
158  std::vector<OutputVector> &eigenvectors,
159  const unsigned int n_eigenpairs = 1);
160 
166  template <typename OutputVector>
167  void
169  const PETScWrappers::MatrixBase &B,
170  std::vector<PetscScalar> &eigenvalues,
171  std::vector<OutputVector> &eigenvectors,
172  const unsigned int n_eigenpairs = 1);
173 
179  template <typename OutputVector>
180  void
182  const PETScWrappers::MatrixBase &B,
183  std::vector<double> &real_eigenvalues,
184  std::vector<double> &imag_eigenvalues,
185  std::vector<OutputVector> &real_eigenvectors,
186  std::vector<OutputVector> &imag_eigenvectors,
187  const unsigned int n_eigenpairs = 1);
188 
192  void
194  (const PETScWrappers::VectorBase &this_initial_vector);
195 
199  void
201 
206  void
207  set_target_eigenvalue (const PetscScalar &this_target);
208 
215  void
216  set_which_eigenpairs (EPSWhich set_which);
217 
226  void
227  set_problem_type (EPSProblemType set_problem);
228 
234  void
236 
240  DeclException0 (ExcSLEPcWrappersUsageError);
241 
245  DeclException1 (ExcSLEPcError,
246  int,
247  << " An error with error number " << arg1
248  << " occurred while calling a SLEPc function");
249 
253  DeclException2 (ExcSLEPcEigenvectorConvergenceMismatchError,
254  int, int,
255  << " The number of converged eigenvectors is " << arg1
256  << " but " << arg2 << " were requested. ");
257 
261  SolverControl &control () const;
262 
263  protected:
264 
270 
274  const MPI_Comm mpi_communicator;
275 
281  virtual void set_solver_type (EPS &eps) const = 0;
282 
286  void
287  reset ();
288 
292  EPS *get_eps ();
293 
301  void
302  solve (const unsigned int n_eigenpairs,
303  unsigned int *n_converged);
304 
310  void
311  get_eigenpair (const unsigned int index,
312  PetscScalar &eigenvalues,
313  PETScWrappers::VectorBase &eigenvectors);
314 
320  void
321  get_eigenpair (const unsigned int index,
322  double &real_eigenvalues,
323  double &imag_eigenvalues,
324  PETScWrappers::VectorBase &real_eigenvectors,
325  PETScWrappers::VectorBase &imag_eigenvectors);
326 
331  void
333 
338  void
340  const PETScWrappers::MatrixBase &B);
341 
345  PetscScalar target_eigenvalue;
346 
350  EPSWhich set_which;
351 
355  EPSProblemType set_problem;
356 
357  private:
358 
365 
371 
376 
382 
389  static
390  int
391  convergence_test (EPS eps,
392  PetscScalar real_eigenvalue,
393  PetscScalar imag_eigenvalue,
394  PetscReal residual_norm,
395  PetscReal *estimated_error,
396  void *solver_control);
397 
406  struct SolverData
407  {
408 
412  ~SolverData ();
413 
417  EPS eps;
418 
422  EPSConvergedReason reason;
423  };
424 
428  std_cxx1x::shared_ptr<SolverData> solver_data;
429  };
430 
441  {
442  public:
443 
449  {};
450 
458  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
459  const AdditionalData &data = AdditionalData());
460 
461  protected:
462 
467 
472  virtual void set_solver_type (EPS &eps) const;
473  };
474 
483  class SolverArnoldi : public SolverBase
484  {
485  public:
491  {
496  AdditionalData (const bool delayed_reorthogonalization = false);
497 
502  };
503 
511  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
512  const AdditionalData &data = AdditionalData());
513 
514  protected:
515 
520 
525  virtual void set_solver_type (EPS &eps) const;
526  };
527 
536  class SolverLanczos : public SolverBase
537  {
538  public:
544  {};
545 
553  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
554  const AdditionalData &data = AdditionalData());
555 
556  protected:
557 
562 
567  virtual void set_solver_type (EPS &eps) const;
568  };
569 
579  class SolverPower : public SolverBase
580  {
581  public:
587  {};
588 
596  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
597  const AdditionalData &data = AdditionalData());
598 
599  protected:
600 
605 
610  virtual void set_solver_type (EPS &eps) const;
611  };
612 
622  {
623  public:
629  {};
630 
638  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
639  const AdditionalData &data = AdditionalData());
640 
641  protected:
642 
647 
652  virtual void set_solver_type (EPS &eps) const;
653  };
654 
664  {
665  public:
671  {};
672 
680  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
681  const AdditionalData &data = AdditionalData());
682 
683  protected:
684 
689 
694  virtual void set_solver_type (EPS &eps) const;
695  };
696 
697 
706  class SolverLAPACK : public SolverBase
707  {
708  public:
709 
715  {};
716 
724  const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
725  const AdditionalData &data = AdditionalData());
726 
727  protected:
728 
733 
738  virtual void set_solver_type (EPS &eps) const;
739  };
740 
741  // --------------------------- inline and template functions -----------
746  // todo: The logic of these functions can be simplified without breaking backward compatibility...
747 
748  template <typename OutputVector>
749  void
751  std::vector<PetscScalar> &eigenvalues,
752  std::vector<OutputVector> &eigenvectors,
753  const unsigned int n_eigenpairs)
754  {
755  // Panic if the number of eigenpairs wanted is out of bounds.
756  AssertThrow ((n_eigenpairs > 0) && (n_eigenpairs <= A.m ()),
757  ExcSLEPcWrappersUsageError());
758 
759  // Set the matrices of the problem
760  set_matrices (A);
761 
762  // and solve
763  unsigned int n_converged = 0;
764  solve (n_eigenpairs, &n_converged);
765 
766  if (n_converged > n_eigenpairs)
767  n_converged = n_eigenpairs;
768  AssertThrow (n_converged == n_eigenpairs,
769  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenpairs));
770 
771  AssertThrow (eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
772  eigenvectors.resize (n_converged, eigenvectors.front());
773  eigenvalues.resize (n_converged);
774 
775  for (unsigned int index=0; index<n_converged; ++index)
776  get_eigenpair (index, eigenvalues[index], eigenvectors[index]);
777  }
778 
779  template <typename OutputVector>
780  void
782  const PETScWrappers::MatrixBase &B,
783  std::vector<PetscScalar> &eigenvalues,
784  std::vector<OutputVector> &eigenvectors,
785  const unsigned int n_eigenpairs)
786  {
787  // Guard against incompatible matrix sizes:
788  AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m()));
789  AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n()));
790 
791  // Panic if the number of eigenpairs wanted is out of bounds.
792  AssertThrow ((n_eigenpairs>0) && (n_eigenpairs<=A.m ()),
793  ExcSLEPcWrappersUsageError());
794 
795  // Set the matrices of the problem
796  set_matrices (A, B);
797 
798  // and solve
799  unsigned int n_converged = 0;
800  solve (n_eigenpairs, &n_converged);
801 
802  if (n_converged>=n_eigenpairs)
803  n_converged = n_eigenpairs;
804 
805  AssertThrow (n_converged==n_eigenpairs,
806  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenpairs));
807  AssertThrow (eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
808 
809  eigenvectors.resize (n_converged, eigenvectors.front());
810  eigenvalues.resize (n_converged);
811 
812  for (unsigned int index=0; index<n_converged; ++index)
813  get_eigenpair (index, eigenvalues[index], eigenvectors[index]);
814  }
815 
816  template <typename OutputVector>
817  void
819  const PETScWrappers::MatrixBase &B,
820  std::vector<double> &real_eigenvalues,
821  std::vector<double> &imag_eigenvalues,
822  std::vector<OutputVector> &real_eigenvectors,
823  std::vector<OutputVector> &imag_eigenvectors,
824  const unsigned int n_eigenpairs)
825  {
826  // Guard against incompatible matrix sizes:
827  AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m()));
828  AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n()));
829 
830  // and incompatible eigenvalue/eigenvector sizes
831  AssertThrow (real_eigenvalues.size() == imag_eigenvalues.size(),
832  ExcDimensionMismatch(real_eigenvalues.size(), imag_eigenvalues.size()));
833  AssertThrow (real_eigenvectors.size() == imag_eigenvectors.size (),
834  ExcDimensionMismatch(real_eigenvectors.size(), imag_eigenvectors.size()));
835 
836  // Panic if the number of eigenpairs wanted is out of bounds.
837  AssertThrow ((n_eigenpairs>0) && (n_eigenpairs<=A.m ()),
838  ExcSLEPcWrappersUsageError());
839 
840  // Set the matrices of the problem
841  set_matrices (A, B);
842 
843  // and solve
844  unsigned int n_converged = 0;
845  solve (n_eigenpairs, &n_converged);
846 
847  if (n_converged>=n_eigenpairs)
848  n_converged = n_eigenpairs;
849 
850  AssertThrow (n_converged==n_eigenpairs,
851  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenpairs));
852  AssertThrow ((real_eigenvectors.size()!=0) && (imag_eigenvectors.size()!=0),
853  ExcSLEPcWrappersUsageError());
854 
855  real_eigenvectors.resize (n_converged, real_eigenvectors.front());
856  imag_eigenvectors.resize (n_converged, imag_eigenvectors.front());
857  real_eigenvalues.resize (n_converged);
858  imag_eigenvalues.resize (n_converged);
859 
860  for (unsigned int index=0; index<n_converged; ++index)
861  get_eigenpair (index,
862  real_eigenvalues[index], imag_eigenvalues[index],
863  real_eigenvectors[index], imag_eigenvectors[index]);
864  }
865 
866 }
867 
868 DEAL_II_NAMESPACE_CLOSE
869 
870 #endif // DEAL_II_WITH_SLEPC
871 
872 /*---------------------------- slepc_solver.h ---------------------------*/
873 
874 #endif
875 
876 /*---------------------------- slepc_solver.h ---------------------------*/
EPSProblemType set_problem
Definition: slepc_solver.h:355
const AdditionalData additional_data
Definition: slepc_solver.h:604
void set_initial_vector(const PETScWrappers::VectorBase &this_initial_vector)
static int convergence_test(EPS eps, PetscScalar real_eigenvalue, PetscScalar imag_eigenvalue, PetscReal residual_norm, PetscReal *estimated_error, void *solver_control)
size_type n() const
size_type m() const
const PETScWrappers::MatrixBase * opB
Definition: slepc_solver.h:370
DeclException0(ExcSLEPcWrappersUsageError)
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
const PETScWrappers::VectorBase * initial_vector
Definition: slepc_solver.h:375
void solve(const PETScWrappers::MatrixBase &A, std::vector< PetscScalar > &eigenvalues, std::vector< OutputVector > &eigenvectors, const unsigned int n_eigenpairs=1)
Definition: slepc_solver.h:750
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:274
SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
void set_problem_type(EPSProblemType set_problem)
DeclException2(ExcSLEPcEigenvectorConvergenceMismatchError, int, int,<< " The number of converged eigenvectors is "<< arg1<< " but "<< arg2<< " were requested. ")
virtual void set_solver_type(EPS &eps) const =0
void set_target_eigenvalue(const PetscScalar &this_target)
void set_matrices(const PETScWrappers::MatrixBase &A)
void set_which_eigenpairs(EPSWhich set_which)
const AdditionalData additional_data
Definition: slepc_solver.h:732
const AdditionalData additional_data
Definition: slepc_solver.h:561
SolverControl & control() const
const AdditionalData additional_data
Definition: slepc_solver.h:466
const PETScWrappers::MatrixBase * opA
Definition: slepc_solver.h:364
void get_eigenpair(const unsigned int index, PetscScalar &eigenvalues, PETScWrappers::VectorBase &eigenvectors)
const AdditionalData additional_data
Definition: slepc_solver.h:688
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const AdditionalData additional_data
Definition: slepc_solver.h:519
std_cxx1x::shared_ptr< SolverData > solver_data
Definition: slepc_solver.h:428
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
SolverControl & solver_control
Definition: slepc_solver.h:269
DeclException1(ExcSLEPcError, int,<< " An error with error number "<< arg1<< " occurred while calling a SLEPc function")
SLEPcWrappers::TransformationBase * transformation
Definition: slepc_solver.h:381
void get_solver_state(const SolverControl::State state)