slepc_solver.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2009 - 2018 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.
//
// ---------------------------------------------------------------------

#include <deal.II/lac/slepc_solver.h>

#ifdef DEAL_II_WITH_SLEPC

#  include <deal.II/lac/petsc_matrix_base.h>
#  include <deal.II/lac/petsc_vector_base.h>
#  include <deal.II/lac/slepc_spectral_transformation.h>

#  include <cmath>
#  include <vector>

#  include <petscversion.h>
#  include <slepcversion.h>

DEAL_II_NAMESPACE_OPEN

namespace SLEPcWrappers
{

  SolverBase::SolverBase (SolverControl  &cn,
                          const MPI_Comm &mpi_communicator)
    :
    solver_control (cn),
    mpi_communicator (mpi_communicator),
    reason (EPS_CONVERGED_ITERATING)
  {
    // create eigensolver context
    PetscErrorCode ierr = EPSCreate (mpi_communicator, &eps);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // hand over the absolute tolerance and the maximum number of
    // iteration steps to the SLEPc convergence criterion.
    ierr = EPSSetTolerances(eps, this->solver_control.tolerance(),
                            this->solver_control.max_steps());
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // default values:
    set_which_eigenpairs(EPS_LARGEST_MAGNITUDE);
    set_problem_type(EPS_GNHEP);

    // TODO:
    // By default, EPS initializes the starting vector or the initial subspace randomly.
  }

  SolverBase::~SolverBase ()
  {
    if (eps != nullptr)
      {
        // Destroy the solver object.
        const PetscErrorCode ierr = EPSDestroy (&eps);

        (void)ierr;
        AssertNothrow (ierr == 0, ExcSLEPcError(ierr));
      }
  }

  void
  SolverBase::set_matrices (const PETScWrappers::MatrixBase &A)
  {
    // standard eigenspectrum problem
    const PetscErrorCode ierr = EPSSetOperators (eps, A, nullptr);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  void
  SolverBase::set_matrices (const PETScWrappers::MatrixBase &A,
                            const PETScWrappers::MatrixBase &B)
  {
    // generalized eigenspectrum problem
    const PetscErrorCode ierr = EPSSetOperators (eps, A, B);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  void
  SolverBase::set_transformation (SLEPcWrappers::TransformationBase &transformation)
  {
    // set transformation type if any
    // STSetShift is called inside
    PetscErrorCode ierr = EPSSetST(eps,transformation.st);
    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));

#if DEAL_II_SLEPC_VERSION_GTE(3, 8, 0)
    // see https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2017-October/033649.html
    // From 3.8.0 SLEPc insists that when looking for smallest eigenvalues with shift-and-invert
    // users should (a) set target (b) use EPS_TARGET_MAGNITUDE
    // The former, however, needs to be applied to eps object and not spectral transformation.
    if (SLEPcWrappers::TransformationShiftInvert *sinv = dynamic_cast<SLEPcWrappers::TransformationShiftInvert *>(&transformation))
      {
        ierr = EPSSetTarget (eps, sinv->additional_data.shift_parameter);
        AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
      }
#endif

  }

  void
  SolverBase::set_initial_vector (const PETScWrappers::VectorBase &this_initial_vector)
  {
    Assert(this_initial_vector.l2_norm()>0.0,
           ExcMessage("Initial vector should be nonzero."));

    Vec vec = this_initial_vector;
    const PetscErrorCode ierr = EPSSetInitialSpace (eps, 1, &vec);

    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  void
  SolverBase::set_target_eigenvalue (const PetscScalar &this_target)
  {
    // set target eigenvalues to solve for
    // in all transformation except STSHIFT there is a direct connection between
    // the target and the shift, read more on p41 of SLEPc manual.
    const PetscErrorCode ierr = EPSSetTarget (eps, this_target );
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  void
  SolverBase::set_which_eigenpairs (const EPSWhich eps_which)
  {
    // set which portion of the eigenspectrum to solve for
    const PetscErrorCode ierr = EPSSetWhichEigenpairs (eps, eps_which);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  void
  SolverBase::set_problem_type (const EPSProblemType eps_problem)
  {
    const PetscErrorCode ierr = EPSSetProblemType (eps, eps_problem);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  void
  SolverBase::solve (const unsigned int  n_eigenpairs,
                     unsigned int       *n_converged)
  {
    // set number of eigenvectors to compute
    PetscErrorCode ierr = EPSSetDimensions (eps, n_eigenpairs,
                                            PETSC_DECIDE, PETSC_DECIDE);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // set the solve options to the eigenvalue problem solver context
    ierr = EPSSetFromOptions (eps);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // TODO breaks @ref step_36 "step-36"
    // force Krylov solver to use true residual instead of an estimate.
    //EPSSetTrueResidual(solver_data->eps, PETSC_TRUE);
    //AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // Set convergence test to be absolute
    ierr = EPSSetConvergenceTest (eps, EPS_CONV_ABS);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // TODO Set the convergence test function
    // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test,
    //              reinterpret_cast<void *>(&solver_control));
    // AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // solve the eigensystem
    ierr = EPSSolve (eps);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // get number of converged eigenstates
    ierr = EPSGetConverged (eps,
                            reinterpret_cast<PetscInt *>(n_converged));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    PetscInt n_iterations   = 0;
    PetscReal residual_norm = 0;

    // @todo Investigate elaborating on some of this to act on the
    // complete eigenspectrum
    {
      // get the number of solver iterations
      ierr = EPSGetIterationNumber (eps, &n_iterations);
      AssertThrow (ierr == 0, ExcSLEPcError(ierr));

      // get the maximum of residual norm among converged eigenvectors.
      for (unsigned int i = 0; i < *n_converged; i++)
        {
          double residual_norm_i = 0.0;
          // EPSComputeError (or, in older versions of SLEPc,
          // EPSComputeResidualNorm) uses an L2-norm and is not consistent
          // with the stopping criterion used during the solution process (see
          // the SLEPC manual, section 2.5). However, the norm that gives error
          // bounds (Saad, 1992, ch3) is (for Hermitian problems)
          // | \lambda - \widehat\lambda | <= ||r||_2
          //
          // Similarly, EPSComputeRelativeError may not be consistent with the
          // stopping criterion used in the solution process.
          //
          // EPSGetErrorEstimate is (according to the SLEPc manual) consistent
          // with the residual norm used during the solution process. However,
          // it is not guaranteed to be derived from the residual even when
          // EPSSetTrueResidual is set: see the discussion in the thread
          //
          // https://lists.mcs.anl.gov/pipermail/petsc-users/2014-November/023509.html
          //
          // for more information.
#if DEAL_II_SLEPC_VERSION_GTE(3, 6, 0)
          ierr = EPSComputeError(eps, i, EPS_ERROR_ABSOLUTE, &residual_norm_i);
#else
          ierr = EPSComputeResidualNorm(eps, i, &residual_norm_i);
#endif

          AssertThrow (ierr == 0, ExcSLEPcError(ierr));
          residual_norm = std::max (residual_norm, residual_norm_i);
        }

      // check the solver state
      const SolverControl::State state
        = solver_control.check (n_iterations, residual_norm);

      // get the solver state according to SLEPc
      get_solver_state (state);

      // as SLEPc uses different stopping criteria, we have to omit this step.
      // This can be checked only in conjunction with EPSGetErrorEstimate.
      // and in case of failure: throw exception
      // if (solver_control.last_check () != SolverControl::success)
      //   AssertThrow(false, SolverControl::NoConvergence (solver_control.last_step(),
      //                                                    solver_control.last_value()));
    }
  }

  void
  SolverBase::get_eigenpair (const unsigned int            index,
                             PetscScalar               &eigenvalues,
                             PETScWrappers::VectorBase &eigenvectors)
  {
    // get converged eigenpair
    const PetscErrorCode ierr = EPSGetEigenpair (eps, index,
                                                 &eigenvalues, nullptr,
                                                 eigenvectors, nullptr);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }


  void
  SolverBase::get_eigenpair (const unsigned int         index,
                             double                    &real_eigenvalues,
                             double                    &imag_eigenvalues,
                             PETScWrappers::VectorBase &real_eigenvectors,
                             PETScWrappers::VectorBase &imag_eigenvectors)
  {
#ifndef PETSC_USE_COMPLEX
    // get converged eigenpair
    const PetscErrorCode ierr = EPSGetEigenpair (eps, index,
                                                 &real_eigenvalues,
                                                 &imag_eigenvalues,
                                                 real_eigenvectors,
                                                 imag_eigenvectors);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
#else
    Assert ((false),
            ExcMessage ("Your PETSc/SLEPc installation was configured with scalar-type complex "
                        "but this function is not defined for complex types."));

    // Cast to void to silence compiler warnings
    (void) index;
    (void) real_eigenvalues;
    (void) imag_eigenvalues;
    (void) real_eigenvectors;
    (void) imag_eigenvectors;
#endif
  }

  void
  SolverBase::get_solver_state (const SolverControl::State state)
  {
    switch (state)
      {
      case ::SolverControl::iterate:
        reason = EPS_CONVERGED_ITERATING;
        break;

      case ::SolverControl::success:
        reason = static_cast<EPSConvergedReason>(1);
        break;

      case ::SolverControl::failure:
        if (solver_control.last_step() > solver_control.max_steps())
          reason = EPS_DIVERGED_ITS;
        else
          reason = EPS_DIVERGED_BREAKDOWN;
        break;

      default:
        Assert (false, ExcNotImplemented());
      }
  }

  /* ---------------------- SolverControls ----------------------- */
  SolverControl &
  SolverBase::control () const
  {
    return solver_control;
  }

  int
  SolverBase::convergence_test (EPS          /*eps             */,
                                PetscScalar  /*real_eigenvalue */,
                                PetscScalar  /*imag_eigenvalue */,
                                PetscReal    /*residual norm associated to the eigenpair   */,
                                PetscReal   */*(output) computed error estimate */,
                                void        */*solver_control_x*/)
  {
    // If the error estimate returned by the convergence test function is less
    // than the tolerance, then the eigenvalue is accepted as converged.
    // This function is undefined (future reference only).

    // return without failure.
    return 0;
  }

  /* ---------------------- SolverKrylovSchur ------------------------ */
  SolverKrylovSchur::SolverKrylovSchur (SolverControl        &cn,
                                        const MPI_Comm       &mpi_communicator,
                                        const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    const PetscErrorCode ierr = EPSSetType (eps,
                                            const_cast<char *>(EPSKRYLOVSCHUR));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  /* ---------------------- SolverArnoldi ------------------------ */
  SolverArnoldi::AdditionalData::
  AdditionalData (const bool delayed_reorthogonalization)
    :
    delayed_reorthogonalization (delayed_reorthogonalization)
  {}

  SolverArnoldi::SolverArnoldi (SolverControl        &cn,
                                const MPI_Comm       &mpi_communicator,
                                const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    PetscErrorCode ierr = EPSSetType (eps, const_cast<char *>(EPSARNOLDI));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    // if requested, set delayed reorthogonalization in the Arnoldi
    // iteration.
    if (additional_data.delayed_reorthogonalization)
      {
        ierr = EPSArnoldiSetDelayed (eps, PETSC_TRUE);
        AssertThrow (ierr == 0, ExcSLEPcError(ierr));
      }
  }


  /* ---------------------- Lanczos ------------------------ */
  SolverLanczos::AdditionalData::
  AdditionalData(const EPSLanczosReorthogType r)
    : reorthog(r)
  {}

  SolverLanczos::SolverLanczos (SolverControl        &cn,
                                const MPI_Comm       &mpi_communicator,
                                const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    PetscErrorCode ierr = EPSSetType (eps, const_cast<char *>(EPSLANCZOS));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    ierr = EPSLanczosSetReorthog(eps,additional_data.reorthog);
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  /* ----------------------- Power ------------------------- */
  SolverPower::SolverPower (SolverControl        &cn,
                            const MPI_Comm       &mpi_communicator,
                            const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    PetscErrorCode ierr = EPSSetType (eps, const_cast<char *>(EPSPOWER));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  /* ---------------- Generalized Davidson ----------------- */
  SolverGeneralizedDavidson::AdditionalData::
  AdditionalData(bool double_expansion)
    :  double_expansion(double_expansion)
  {}

  SolverGeneralizedDavidson::SolverGeneralizedDavidson (SolverControl        &cn,
                                                        const MPI_Comm       &mpi_communicator,
                                                        const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    PetscErrorCode ierr = EPSSetType (eps, const_cast<char *>(EPSGD));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));

    if (additional_data.double_expansion)
      {
        ierr = EPSGDSetDoubleExpansion (eps, PETSC_TRUE);
        AssertThrow (ierr == 0, ExcSLEPcError(ierr));
      }
  }

  /* ------------------ Jacobi Davidson -------------------- */
  SolverJacobiDavidson::SolverJacobiDavidson (SolverControl        &cn,
                                              const MPI_Comm       &mpi_communicator,
                                              const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    const PetscErrorCode ierr = EPSSetType (eps, const_cast<char *>(EPSJD));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
  }

  /* ---------------------- LAPACK ------------------------- */
  SolverLAPACK::SolverLAPACK (SolverControl        &cn,
                              const MPI_Comm       &mpi_communicator,
                              const AdditionalData &data)
    :
    SolverBase (cn, mpi_communicator),
    additional_data (data)
  {
    // 'Tis overwhelmingly likely that PETSc/SLEPc *always* has
    // BLAS/LAPACK, but let's be defensive.
#if PETSC_HAVE_BLASLAPACK
    const PetscErrorCode ierr = EPSSetType (eps, const_cast<char *>(EPSLAPACK));
    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
#else
    Assert ((false),
            ExcMessage ("Your PETSc/SLEPc installation was not configured with BLAS/LAPACK "
                        "but this is needed to use the LAPACK solver."));
#endif
  }
}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_SLEPC