slepc_solver.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2009 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
 
#ifndef dealii_slepc_solver_h
#  define dealii_slepc_solver_h
 
#  include <deal.II/base/config.h>
 
#  ifdef DEAL_II_WITH_SLEPC
 
#    include <deal.II/lac/petsc_matrix_base.h>
#    include <deal.II/lac/slepc_spectral_transformation.h>
#    include <deal.II/lac/solver_control.h>
 
#    include <petscconf.h>
#    include <petscksp.h>
 
#    include <slepceps.h>
 
#    include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
namespace SLEPcWrappers
{
  class SolverBase
  {
  public:
    SolverBase(SolverControl &cn, const MPI_Comm mpi_communicator);
 
    virtual ~SolverBase();
 
    template <typename OutputVector>
    void
    solve(const PETScWrappers::MatrixBase &A,
          std::vector<PetscScalar>        &eigenvalues,
          std::vector<OutputVector>       &eigenvectors,
          const unsigned int               n_eigenpairs = 1);
 
    template <typename OutputVector>
    void
    solve(const PETScWrappers::MatrixBase &A,
          const PETScWrappers::MatrixBase &B,
          std::vector<PetscScalar>        &eigenvalues,
          std::vector<OutputVector>       &eigenvectors,
          const unsigned int               n_eigenpairs = 1);
 
    template <typename OutputVector>
    void
    solve(const PETScWrappers::MatrixBase &A,
          const PETScWrappers::MatrixBase &B,
          std::vector<double>             &real_eigenvalues,
          std::vector<double>             &imag_eigenvalues,
          std::vector<OutputVector>       &real_eigenvectors,
          std::vector<OutputVector>       &imag_eigenvectors,
          const unsigned int               n_eigenpairs = 1);
 
    template <typename Vector>
    void
    set_initial_space(const std::vector<Vector> &initial_space);
 
    void
    set_transformation(SLEPcWrappers::TransformationBase &this_transformation);
 
    void
    set_target_eigenvalue(const PetscScalar &this_target);
 
    void
    set_which_eigenpairs(EPSWhich set_which);
 
    void
    set_problem_type(EPSProblemType set_problem);
 
    void
    get_solver_state(const SolverControl::State state);
 
    DeclException0(ExcSLEPcWrappersUsageError);
 
    DeclException1(ExcSLEPcError,
                   int,
                   << "    An error with error number " << arg1
                   << " occurred while calling a SLEPc function");
 
    DeclException2(ExcSLEPcEigenvectorConvergenceMismatchError,
                   int,
                   int,
                   << "    The number of converged eigenvectors is " << arg1
                   << " but " << arg2 << " were requested. ");
 
    SolverControl &
    control() const;
 
  protected:
    SolverControl &solver_control;
 
    const MPI_Comm mpi_communicator;
 
    void
    solve(const unsigned int n_eigenpairs, unsigned int *n_converged);
 
    void
    get_eigenpair(const unsigned int         index,
                  PetscScalar               &eigenvalues,
                  PETScWrappers::VectorBase &eigenvectors);
 
    void
    get_eigenpair(const unsigned int         index,
                  double                    &real_eigenvalues,
                  double                    &imag_eigenvalues,
                  PETScWrappers::VectorBase &real_eigenvectors,
                  PETScWrappers::VectorBase &imag_eigenvectors);
 
    void
    set_matrices(const PETScWrappers::MatrixBase &A);
 
    void
    set_matrices(const PETScWrappers::MatrixBase &A,
                 const PETScWrappers::MatrixBase &B);
 
  protected:
    EPS eps;
 
  private:
    EPSConvergedReason reason;
 
 
    static int
    convergence_test(EPS         eps,
                     PetscScalar real_eigenvalue,
                     PetscScalar imag_eigenvalue,
                     PetscReal   residual_norm,
                     PetscReal  *estimated_error,
                     void       *solver_control);
  };
 
 
 
  class SolverKrylovSchur : public SolverBase
  {
  public:
    struct AdditionalData
    {};
 
    explicit SolverKrylovSchur(
      SolverControl        &cn,
      const MPI_Comm        mpi_communicator = PETSC_COMM_SELF,
      const AdditionalData &data             = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  class SolverArnoldi : public SolverBase
  {
  public:
    struct AdditionalData
    {
      explicit AdditionalData(const bool delayed_reorthogonalization = false);
 
      bool delayed_reorthogonalization;
    };
 
    explicit SolverArnoldi(SolverControl &cn,
                           const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
                           const AdditionalData &data      = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  class SolverLanczos : public SolverBase
  {
  public:
    struct AdditionalData
    {
      EPSLanczosReorthogType reorthog;
 
      explicit AdditionalData(
        const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
    };
 
    explicit SolverLanczos(SolverControl &cn,
                           const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
                           const AdditionalData &data      = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  class SolverPower : public SolverBase
  {
  public:
    struct AdditionalData
    {};
 
    explicit SolverPower(SolverControl &cn,
                         const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
                         const AdditionalData &data      = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  class SolverGeneralizedDavidson : public SolverBase
  {
  public:
    struct AdditionalData
    {
      bool double_expansion;
 
      explicit AdditionalData(bool double_expansion = false);
    };
 
    explicit SolverGeneralizedDavidson(
      SolverControl        &cn,
      const MPI_Comm        mpi_communicator = PETSC_COMM_SELF,
      const AdditionalData &data             = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  class SolverJacobiDavidson : public SolverBase
  {
  public:
    struct AdditionalData
    {};
 
    explicit SolverJacobiDavidson(
      SolverControl        &cn,
      const MPI_Comm        mpi_communicator = PETSC_COMM_SELF,
      const AdditionalData &data             = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  class SolverLAPACK : public SolverBase
  {
  public:
    struct AdditionalData
    {};
 
    explicit SolverLAPACK(SolverControl &cn,
                          const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
                          const AdditionalData &data      = AdditionalData());
 
  protected:
    const AdditionalData additional_data;
  };
 
 
 
  // --------------------------- inline and template functions -----------
  // todo: The logic of these functions can be simplified without breaking
  // backward compatibility...
 
  template <typename OutputVector>
  void
  SolverBase::solve(const PETScWrappers::MatrixBase &A,
                    std::vector<PetscScalar>        &eigenvalues,
                    std::vector<OutputVector>       &eigenvectors,
                    const unsigned int               n_eigenpairs)
  {
    // Panic if the number of eigenpairs wanted is out of bounds.
    AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
                ExcSLEPcWrappersUsageError());
 
    // Set the matrices of the problem
    set_matrices(A);
 
    // and solve
    unsigned int n_converged = 0;
    solve(n_eigenpairs, &n_converged);
 
    if (n_converged > n_eigenpairs)
      n_converged = n_eigenpairs;
    AssertThrow(n_converged == n_eigenpairs,
                ExcSLEPcEigenvectorConvergenceMismatchError(n_converged,
                                                            n_eigenpairs));
 
    AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
    eigenvectors.resize(n_converged, eigenvectors.front());
    eigenvalues.resize(n_converged);
 
    for (unsigned int index = 0; index < n_converged; ++index)
      get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
  }
 
  template <typename OutputVector>
  void
  SolverBase::solve(const PETScWrappers::MatrixBase &A,
                    const PETScWrappers::MatrixBase &B,
                    std::vector<PetscScalar>        &eigenvalues,
                    std::vector<OutputVector>       &eigenvectors,
                    const unsigned int               n_eigenpairs)
  {
    // Guard against incompatible matrix sizes:
    AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
    AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
 
    // Panic if the number of eigenpairs wanted is out of bounds.
    AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
                ExcSLEPcWrappersUsageError());
 
    // Set the matrices of the problem
    set_matrices(A, B);
 
    // and solve
    unsigned int n_converged = 0;
    solve(n_eigenpairs, &n_converged);
 
    if (n_converged >= n_eigenpairs)
      n_converged = n_eigenpairs;
 
    AssertThrow(n_converged == n_eigenpairs,
                ExcSLEPcEigenvectorConvergenceMismatchError(n_converged,
                                                            n_eigenpairs));
    AssertThrow(eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
 
    eigenvectors.resize(n_converged, eigenvectors.front());
    eigenvalues.resize(n_converged);
 
    for (unsigned int index = 0; index < n_converged; ++index)
      get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
  }
 
  template <typename OutputVector>
  void
  SolverBase::solve(const PETScWrappers::MatrixBase &A,
                    const PETScWrappers::MatrixBase &B,
                    std::vector<double>             &real_eigenvalues,
                    std::vector<double>             &imag_eigenvalues,
                    std::vector<OutputVector>       &real_eigenvectors,
                    std::vector<OutputVector>       &imag_eigenvectors,
                    const unsigned int               n_eigenpairs)
  {
    // Guard against incompatible matrix sizes:
    AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
    AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
 
    // and incompatible eigenvalue/eigenvector sizes
    AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
                ExcDimensionMismatch(real_eigenvalues.size(),
                                     imag_eigenvalues.size()));
    AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
                ExcDimensionMismatch(real_eigenvectors.size(),
                                     imag_eigenvectors.size()));
 
    // Panic if the number of eigenpairs wanted is out of bounds.
    AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
                ExcSLEPcWrappersUsageError());
 
    // Set the matrices of the problem
    set_matrices(A, B);
 
    // and solve
    unsigned int n_converged = 0;
    solve(n_eigenpairs, &n_converged);
 
    if (n_converged >= n_eigenpairs)
      n_converged = n_eigenpairs;
 
    AssertThrow(n_converged == n_eigenpairs,
                ExcSLEPcEigenvectorConvergenceMismatchError(n_converged,
                                                            n_eigenpairs));
    AssertThrow((real_eigenvectors.size() != 0) &&
                  (imag_eigenvectors.size() != 0),
                ExcSLEPcWrappersUsageError());
 
    real_eigenvectors.resize(n_converged, real_eigenvectors.front());
    imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
    real_eigenvalues.resize(n_converged);
    imag_eigenvalues.resize(n_converged);
 
    for (unsigned int index = 0; index < n_converged; ++index)
      get_eigenpair(index,
                    real_eigenvalues[index],
                    imag_eigenvalues[index],
                    real_eigenvectors[index],
                    imag_eigenvectors[index]);
  }
 
  template <typename VectorType>
  void
  SolverBase::set_initial_space(
    const std::vector<VectorType> &this_initial_space)
  {
    std::vector<Vec> vecs(this_initial_space.size());
 
    for (unsigned int i = 0; i < this_initial_space.size(); ++i)
      {
        Assert(this_initial_space[i].l2_norm() > 0.0,
               ExcMessage("Initial vectors should be nonzero."));
        vecs[i] = this_initial_space[i];
      }
 
    // if the eigensolver supports only a single initial vector, but several
    // guesses are provided, then all except the first one will be discarded.
    // One could still build a vector that is rich in the directions of all
    // guesses, by taking a linear combination of them. (TODO: make function
    // virtual?)
 
    const PetscErrorCode ierr =
      EPSSetInitialSpace(eps, vecs.size(), vecs.data());
    AssertThrow(ierr == 0, ExcSLEPcError(ierr));
  }
 
} // namespace SLEPcWrappers
 
DEAL_II_NAMESPACE_CLOSE
 
#  endif // DEAL_II_WITH_SLEPC
 
/*----------------------------   slepc_solver.h  ---------------------------*/
 
#endif
 
/*----------------------------   slepc_solver.h  ---------------------------*/