slepc_solver.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2009 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#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/solver_control.h>
#  include <deal.II/lac/petsc_matrix_base.h>
#  include <deal.II/lac/slepc_spectral_transformation.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);

    DEAL_II_DEPRECATED
    void
    set_initial_vector
    (const PETScWrappers::VectorBase &this_initial_vector);

    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
    {};

    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
    {
      AdditionalData (const bool delayed_reorthogonalization = false);

      bool delayed_reorthogonalization;
    };

    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;

      AdditionalData(const EPSLanczosReorthogType r  = EPS_LANCZOS_REORTHOG_FULL);
    };

    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
    {};

    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;

      AdditionalData(bool double_expansion = false);
    };

    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
    {};

    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
    {};

    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 Vector>
  void
  SolverBase::set_initial_space(const std::vector<Vector> &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));
  }

}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_SLEPC

/*----------------------------   slepc_solver.h  ---------------------------*/

#endif

/*----------------------------   slepc_solver.h  ---------------------------*/