parpack_solver.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2016 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__parpack_solver_h
#define dealii__parpack_solver_h

#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/vector.h>
#include <deal.II/base/index_set.h>

#include <cstring>


#ifdef DEAL_II_ARPACK_WITH_PARPACK

DEAL_II_NAMESPACE_OPEN

extern "C" {

  // http://www.mathkeisan.com/usersguide/man/pdnaupd.html
  void pdnaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which,
                int *nev, double *tol, double *resid, int *ncv,
                double *v, int *nloc, int *iparam, int *ipntr,
                double *workd, double *workl, int *lworkl,
                int *info);

  // http://www.mathkeisan.com/usersguide/man/pdsaupd.html
  void pdsaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which,
                int *nev, double *tol, double *resid, int *ncv,
                double *v, int *nloc, int *iparam, int *ipntr,
                double *workd, double *workl, int *lworkl,
                int *info);

  // http://www.mathkeisan.com/usersguide/man/pdneupd.html
  void pdneupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d,
                double *di, double *z, int *ldz, double *sigmar,
                double *sigmai, double *workev, char *bmat, int *n, char *which,
                int *nev, double *tol, double *resid, int *ncv,
                double *v, int *nloc, int *iparam, int *ipntr,
                double *workd, double *workl, int *lworkl, int *info);

  // http://www.mathkeisan.com/usersguide/man/pdseupd.html
  void pdseupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d,
                double *z, int *ldz, double *sigmar,
                char *bmat, int *n, char *which,
                int *nev, double *tol, double *resid, int *ncv,
                double *v, int *nloc, int *iparam, int *ipntr,
                double *workd, double *workl, int *lworkl, int *info);

  // other resources:
  //    http://acts.nersc.gov/superlu/example5/pnslac.c.html
  //    https://github.com/phpisciuneri/tijo/blob/master/dvr_parpack.cpp

}

template <typename VectorType>
class PArpackSolver : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  enum WhichEigenvalues
  {
    algebraically_largest,
    algebraically_smallest,
    largest_magnitude,
    smallest_magnitude,
    largest_real_part,
    smallest_real_part,
    largest_imaginary_part,
    smallest_imaginary_part,
    both_ends
  };

  template <typename MatrixType>
  class Shift : public ::Subscriptor
  {
  public:

    Shift (const MatrixType &A,
           const MatrixType &B,
           const double      sigma)
      :
      A(A),
      B(B),
      sigma(sigma)
    {}

    void vmult (VectorType &dst, const VectorType &src) const
    {
      B.vmult(dst,src);
      dst *= (-sigma);
      A.vmult_add(dst,src);
    }

    void Tvmult (VectorType &dst, const VectorType &src) const
    {
      B.Tvmult(dst,src);
      dst *= (-sigma);
      A.Tvmult_add(dst,src);
    }

  private:
    const MatrixType &A;
    const MatrixType &B;
    const double sigma;
  };

  struct AdditionalData
  {
    const unsigned int number_of_arnoldi_vectors;
    const WhichEigenvalues eigenvalue_of_interest;
    const bool symmetric;
    AdditionalData(
      const unsigned int number_of_arnoldi_vectors = 15,
      const WhichEigenvalues eigenvalue_of_interest = largest_magnitude,
      const bool symmetric = false);
  };

  SolverControl &control () const;

  PArpackSolver(SolverControl &control,
                const MPI_Comm &mpi_communicator,
                const AdditionalData &data = AdditionalData());

  void reinit(const ::IndexSet &locally_owned_dofs );

  void set_shift(const double s );

  template <typename MatrixType1,
            typename MatrixType2, typename INVERSE>
  void solve
  (const MatrixType1                  &A,
   const MatrixType2                  &B,
   const INVERSE                      &inverse,
   std::vector<std::complex<double> > &eigenvalues,
   std::vector<VectorType>            &eigenvectors,
   const unsigned int                  n_eigenvalues);

  std::size_t memory_consumption() const;

protected:

  SolverControl &solver_control;

  const AdditionalData additional_data;

  // keep MPI communicator non-const as Arpack functions are not const either:

  MPI_Comm mpi_communicator;

  MPI_Fint mpi_communicator_fortran;

  // C++98 guarantees that the elements of a vector are stored contiguously

  int lworkl;

  std::vector<double> workl;

  std::vector<double> workd;

  int nloc;

  int ncv;


  int ldv;

  std::vector<double> v;

  std::vector<double> resid;

  int ldz;

  std::vector<double> z;

  int lworkev;

  std::vector<double> workev;

  std::vector<int> select;

  VectorType src,dst,tmp;

  std::vector< types::global_dof_index > local_indices;

  double shift_value;

private:

  DeclException2 (PArpackExcConvergedEigenvectors, int, int,
                  << arg1 << "eigenpairs were requested, but only"
                  << arg2 << " converged");

  DeclException2 (PArpackExcInvalidNumberofEigenvalues, int, int,
                  << "Number of wanted eigenvalues " << arg1
                  << " is larger that the size of the matrix " << arg2);

  DeclException2 (PArpackExcInvalidEigenvectorSize, int, int,
                  << "Number of wanted eigenvalues " << arg1
                  << " is larger that the size of eigenvectors " << arg2);

  DeclException2 (PArpackExcInvalidEigenvalueSize, int, int,
                  << "Number of wanted eigenvalues " << arg1
                  << " is larger that the size of eigenvalues " << arg2);

  DeclException2 (PArpackExcInvalidNumberofArnoldiVectors, int, int,
                  << "Number of Arnoldi vectors " << arg1
                  << " is larger that the size of the matrix " << arg2);

  DeclException2 (PArpackExcSmallNumberofArnoldiVectors, int, int,
                  << "Number of Arnoldi vectors " << arg1
                  << " is too small to obtain " << arg2
                  << " eigenvalues");

  DeclException1 (PArpackExcIdo, int, << "This ido " << arg1
                  << " is not supported. Check documentation of ARPACK");

  DeclException1 (PArpackExcMode, int, << "This mode " << arg1
                  << " is not supported. Check documentation of ARPACK");

  DeclException1 (PArpackExcInfoPdnaupd, int,
                  << "Error with Pdnaupd, info " << arg1
                  << ". Check documentation of ARPACK");

  DeclException1 (PArpackExcInfoPdneupd, int,
                  << "Error with Pdneupd, info " << arg1
                  << ". Check documentation of ARPACK");

  DeclException1 (PArpackExcInfoMaxIt, int,
                  << "Maximum number " << arg1
                  << " of iterations reached.");

  DeclException1 (PArpackExcNoShifts, int,
                  << "No shifts could be applied during implicit"
                  << " Arnoldi update, try increasing the number of"
                  << " Arnoldi vectors.");
};

template <typename VectorType>
std::size_t
PArpackSolver<VectorType>::memory_consumption() const
{
  return  MemoryConsumption::memory_consumption (double()) *
          (workl.size()  +
           workd.size()  +
           v.size()      +
           resid.size()  +
           z.size()      +
           workev.size()  )       +
          src.memory_consumption() +
          dst.memory_consumption() +
          tmp.memory_consumption() +
          MemoryConsumption::memory_consumption (types::global_dof_index()) * local_indices.size();
}

template <typename VectorType>
PArpackSolver<VectorType>::AdditionalData::
AdditionalData (const unsigned int     number_of_arnoldi_vectors,
                const WhichEigenvalues eigenvalue_of_interest,
                const bool             symmetric)
  :
  number_of_arnoldi_vectors(number_of_arnoldi_vectors),
  eigenvalue_of_interest(eigenvalue_of_interest),
  symmetric(symmetric)
{}

template <typename VectorType>
PArpackSolver<VectorType>::PArpackSolver (SolverControl        &control,
                                          const MPI_Comm       &mpi_communicator,
                                          const AdditionalData &data)
  :
  solver_control (control),
  additional_data (data),
  mpi_communicator( mpi_communicator ),
  mpi_communicator_fortran ( MPI_Comm_c2f( mpi_communicator ) ),
  shift_value(0.0)

{}

template <typename VectorType>
void PArpackSolver<VectorType>::set_shift(const double s )
{
  shift_value = s;
}

template <typename VectorType>
void PArpackSolver<VectorType>::reinit(const ::IndexSet &locally_owned_dofs)
{
  // store local indices to write to vectors
  locally_owned_dofs.fill_index_vector(local_indices);

  // scalars
  nloc = locally_owned_dofs.n_elements ();
  ncv  = additional_data.number_of_arnoldi_vectors;

  Assert (local_indices.size() == nloc, ExcInternalError() );

  // vectors
  ldv = nloc;
  v.resize (ldv*ncv, 0.0);

  // TODO: add optional input for resid
  resid.resize(nloc, 1.0);

  // work arrays for ARPACK
  workd.resize(3*nloc,0.0);

  lworkl = additional_data.symmetric ?
           ncv*ncv + 8*ncv
           :
           3*ncv*ncv+6*ncv;
  workl.resize (lworkl, 0.);

  ldz = nloc;
  z.resize (ldz*ncv, 0.); // TODO we actually need only ldz*nev

  // WORKEV  Double precision  work array of dimension 3*NCV.
  lworkev = additional_data.symmetric ?
            0 /*not used in symmetric case*/
            :
            3*ncv;
  workev.resize (lworkev, 0.);

  select.resize (ncv, 0);

  // deal.II vectors:
  src.reinit (locally_owned_dofs,mpi_communicator);
  dst.reinit (locally_owned_dofs,mpi_communicator);
  tmp.reinit (locally_owned_dofs,mpi_communicator);

}

template <typename VectorType>
template <typename MatrixType1,typename MatrixType2, typename INVERSE>
void PArpackSolver<VectorType>::solve
(const MatrixType1                  &/*system_matrix*/,
 const MatrixType2                  &mass_matrix,
 const INVERSE                      &inverse,
 std::vector<std::complex<double> > &eigenvalues,
 std::vector<VectorType>            &eigenvectors,
 const unsigned int                  n_eigenvalues)
{

  Assert (n_eigenvalues <= eigenvectors.size(),
          PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size()));

  Assert (n_eigenvalues <= eigenvalues.size(),
          PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));


  Assert (n_eigenvalues < mass_matrix.m(),
          PArpackExcInvalidNumberofEigenvalues(n_eigenvalues, mass_matrix.m()));

  Assert (additional_data.number_of_arnoldi_vectors < mass_matrix.m(),
          PArpackExcInvalidNumberofArnoldiVectors(
            additional_data.number_of_arnoldi_vectors, mass_matrix.m()));

  Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1,
          PArpackExcSmallNumberofArnoldiVectors(
            additional_data.number_of_arnoldi_vectors, n_eigenvalues));
  // ARPACK mode for dnaupd, here only
  //  Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
  //c           ===> OP = (inv[K - sigma*M])*M  and  B = M.
  //c           ===> Shift-and-Invert mode
  int mode = 3;

  // reverse communication parameter
  // must be zero on the first call to pdnaupd
  int ido = 0;

  // 'G' generalized eigenvalue problem
  // 'I' standard eigenvalue problem
  char bmat[2] = "G";

  // Specify the eigenvalues of interest, possible parameters:
  // "LA" algebraically largest
  // "SA" algebraically smallest
  // "LM" largest magnitude
  // "SM" smallest magnitude
  // "LR" largest real part
  // "SR" smallest real part
  // "LI" largest imaginary part
  // "SI" smallest imaginary part
  // "BE" both ends of spectrum simultaneous

  char which[3];
  switch (additional_data.eigenvalue_of_interest)
    {
    case algebraically_largest:
      std::strcpy (which, "LA");
      break;
    case algebraically_smallest:
      std::strcpy (which, "SA");
      break;
    case largest_magnitude:
      std::strcpy (which, "LM");
      break;
    case smallest_magnitude:
      std::strcpy (which, "SM");
      break;
    case largest_real_part:
      std::strcpy (which, "LR");
      break;
    case smallest_real_part:
      std::strcpy (which, "SR");
      break;
    case largest_imaginary_part:
      std::strcpy (which, "LI");
      break;
    case smallest_imaginary_part:
      std::strcpy (which, "SI");
      break;
    case both_ends:
      std::strcpy (which, "BE");
      break;
    }

  // tolerance for ARPACK
  double tol = control().tolerance();

  //information to the routines
  std::vector<int> iparam (11, 0);

  iparam[0] = 1;
  // shift strategy: exact shifts with respect to the current Hessenberg matrix H.

  // maximum number of iterations
  iparam[2] = control().max_steps();

  // Parpack currently works only for NB = 1
  iparam[3] = 1;

  // Sets the mode of dsaupd:
  // 1 is exact shifting,
  // 2 is user-supplied shifts,
  // 3 is shift-invert mode,
  // 4 is buckling mode,
  // 5 is Cayley mode.

  iparam[6] = mode;
  std::vector<int> ipntr (14, 0);

  //information out of the iteration
  //  If INFO .EQ. 0, a random initial residual vector is used.
  //  If INFO .NE. 0, RESID contains the initial residual vector,
  //  possibly from a previous run.
  // Typical choices in this situation might be to use the final value
  // of the starting vector from the previous eigenvalue calculation
  int info = 1;

  // Number of eigenvalues of OP to be computed. 0 < NEV < N.
  int nev = n_eigenvalues;
  int n_inside_arpack = nloc;

  while (ido != 99)
    {
      // call of ARPACK pdnaupd routine
      if (additional_data.symmetric)
        pdsaupd_(&mpi_communicator_fortran,&ido, bmat, &n_inside_arpack, which, &nev, &tol,
                 &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
                 &workd[0], &workl[0], &lworkl, &info);
      else
        pdnaupd_(&mpi_communicator_fortran,&ido, bmat, &n_inside_arpack, which, &nev, &tol,
                 &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
                 &workd[0], &workl[0], &lworkl, &info);

      if (ido == 99)
        break;

      switch (mode)
        {
//        OP = (inv[K - sigma*M])*M
        case 3:
        {
          switch (ido)
            {
//            compute  Y = OP * X  where
//            IPNTR(1) is the pointer into WORKD for X,
//            IPNTR(2) is the pointer into WORKD for Y.
            case -1:
            {
              const int shift_x = ipntr[0]-1;
              const int shift_y = ipntr[1]-1;
              Assert (shift_x>=0, ::ExcInternalError() );
              Assert (shift_x+nloc <= workd.size(), ::ExcInternalError() );
              Assert (shift_y>=0, ::ExcInternalError() );
              Assert (shift_y+nloc <= workd.size(), ::ExcInternalError() );

              src = 0.0;
              src.add (nloc,
                       &local_indices[0],
                       &workd[0]+shift_x );
              src.compress (VectorOperation::add);

              // multiplication with mass matrix M
              mass_matrix.vmult(tmp, src);
              // solving linear system
              inverse.vmult(dst,tmp);

              // store the result
              dst.extract_subvector_to (local_indices.begin(),
                                        local_indices.end(),
                                        &workd[0]+shift_y  );
            }
            break;

//            compute  Y = OP * X where
//            IPNTR(1) is the pointer into WORKD for X,
//            IPNTR(2) is the pointer into WORKD for Y.
//            In mode 3,4 and 5, the vector B * X is already
//            available in WORKD(ipntr(3)).  It does not
//            need to be recomputed in forming OP * X.
            case  1:
            {
              const int shift_x   = ipntr[0]-1;
              const int shift_y   = ipntr[1]-1;
              const int shift_b_x = ipntr[2]-1;

              Assert (shift_x>=0, ::ExcInternalError() );
              Assert (shift_x+nloc <= workd.size(), ::ExcInternalError() );
              Assert (shift_y>=0, ::ExcInternalError() );
              Assert (shift_y+nloc <= workd.size(), ::ExcInternalError() );
              Assert (shift_b_x>=0, ::ExcInternalError() );
              Assert (shift_b_x+nloc <= workd.size(), ::ExcInternalError() );
              Assert (shift_y>=0, ::ExcInternalError() );
              Assert (shift_y+nloc <= workd.size(), ::ExcInternalError() );

              src = 0.0; // B*X
              src.add (nloc,
                       &local_indices[0],
                       &workd[0]+shift_b_x );

              tmp = 0.0; // X
              tmp.add (nloc,
                       &local_indices[0],
                       &workd[0]+shift_x);

              src.compress (VectorOperation::add);
              tmp.compress (VectorOperation::add);

              // solving linear system
              inverse.vmult(dst,src);

              // store the result
              dst.extract_subvector_to (local_indices.begin(),
                                        local_indices.end(),
                                        &workd[0]+shift_y  );

            }
            break;

//            compute  Y = B * X  where
//            IPNTR(1) is the pointer into WORKD for X,
//            IPNTR(2) is the pointer into WORKD for Y.
            case  2:
            {

              const int shift_x = ipntr[0]-1;
              const int shift_y = ipntr[1]-1;
              Assert (shift_x>=0, ::ExcInternalError() );
              Assert (shift_x+nloc <= workd.size(), ::ExcInternalError() );
              Assert (shift_y>=0, ::ExcInternalError() );
              Assert (shift_y+nloc <= workd.size(), ::ExcInternalError() );

              src = 0.0;
              src.add (nloc,
                       &local_indices[0],
                       &workd[0]+shift_x );
              src.compress (VectorOperation::add);

              // Multiplication with mass matrix M
              mass_matrix.vmult(dst, src);

              // store the result
              dst.extract_subvector_to (local_indices.begin(),
                                        local_indices.end(),
                                        &workd[0]+shift_y);

            }
            break;

            default:
              Assert (false, PArpackExcIdo(ido));
              break;
            }
        }
        break;
        default:
          Assert (false, PArpackExcMode(mode));
          break;
        }
    }

  if (info<0)
    {
      Assert (false, PArpackExcInfoPdnaupd(info));
    }
  else
    {
      // 1 - compute eigenvectors,
      // 0 - only eigenvalues
      int rvec = 1;

      // which eigenvectors
      char howmany[4] = "All";

      double sigmar = shift_value; // real part of the shift
      double sigmai = 0.0; // imaginary part of the shift

      std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
      std::vector<double> eigenvalues_im (n_eigenvalues, 0.);

      // call of ARPACK pdneupd routine
      if (additional_data.symmetric)
        pdseupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0],
                 &z[0], &ldz, &sigmar,
                 bmat, &n_inside_arpack, which, &nev, &tol,
                 &resid[0], &ncv, &v[0], &ldv,
                 &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
      else
        pdneupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0],
                 &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai,
                 &workev[0], bmat, &n_inside_arpack, which, &nev, &tol,
                 &resid[0], &ncv, &v[0], &ldv,
                 &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);

      if (info == 1)
        {
          Assert (false, PArpackExcInfoMaxIt(control().max_steps()));
        }
      else if (info == 3)
        {
          Assert (false, PArpackExcNoShifts(1));
        }
      else if (info!=0)
        {
          Assert (false, PArpackExcInfoPdneupd(info));
        }

      for (size_type i=0; i<n_eigenvalues; ++i)
        {
          eigenvectors[i] = 0.0;
          Assert (i*nloc + nloc <= v.size(), ::ExcInternalError() );

          eigenvectors[i].add (nloc,
                               &local_indices[0],
                               &v[i*nloc] );
          eigenvectors[i].compress (VectorOperation::add);
        }

      for (size_type i=0; i<n_eigenvalues; ++i)
        eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
                                               eigenvalues_im[i]);
    }

  Assert (iparam[4] == n_eigenvalues,
          PArpackExcConvergedEigenvectors(iparam[4], n_eigenvalues));

  // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
  // with respect to a semi-inner product defined by M.

  // resid likely contains residual with respect to M-norm.
  {

    tmp = 0.0;
    tmp.add (nloc,
             &local_indices[0],
             &resid[0]);
    solver_control.check  ( iparam[2], tmp.l2_norm() );
  }


}

template <typename VectorType>
SolverControl &PArpackSolver<VectorType>::control () const
{
  return solver_control;
}

DEAL_II_NAMESPACE_CLOSE


#endif
#endif