ArpackSelfAdjointEigenSolver.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 David Harmon <dharmon@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
#ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
#define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
 
#include "../../../../Eigen/Dense"
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
namespace internal {
template <typename Scalar, typename RealScalar>
struct arpack_wrapper;
template <typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD>
struct OP;
}  // namespace internal
 
template <typename MatrixType, typename MatrixSolver = SimplicialLLT<MatrixType>, bool BisSPD = false>
class ArpackGeneralizedSelfAdjointEigenSolver {
 public:
  // typedef typename MatrixSolver::MatrixType MatrixType;

  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::Index Index;

  typedef typename NumTraits<Scalar>::Real RealScalar;

  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;

  ArpackGeneralizedSelfAdjointEigenSolver()
      : m_eivec(),
        m_eivalues(),
        m_isInitialized(false),
        m_eigenvectorsOk(false),
        m_nbrConverged(0),
        m_nbrIterations(0) {}

  ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType &A, const MatrixType &B, Index nbrEigenvalues,
                                          std::string eigs_sigma = "LM", int options = ComputeEigenvectors,
                                          RealScalar tol = 0.0)
      : m_eivec(),
        m_eivalues(),
        m_isInitialized(false),
        m_eigenvectorsOk(false),
        m_nbrConverged(0),
        m_nbrIterations(0) {
    compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
  }

 
  ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType &A, Index nbrEigenvalues, std::string eigs_sigma = "LM",
                                          int options = ComputeEigenvectors, RealScalar tol = 0.0)
      : m_eivec(),
        m_eivalues(),
        m_isInitialized(false),
        m_eigenvectorsOk(false),
        m_nbrConverged(0),
        m_nbrIterations(0) {
    compute(A, nbrEigenvalues, eigs_sigma, options, tol);
  }

  ArpackGeneralizedSelfAdjointEigenSolver &compute(const MatrixType &A, const MatrixType &B, Index nbrEigenvalues,
                                                   std::string eigs_sigma = "LM", int options = ComputeEigenvectors,
                                                   RealScalar tol = 0.0);

  ArpackGeneralizedSelfAdjointEigenSolver &compute(const MatrixType &A, Index nbrEigenvalues,
                                                   std::string eigs_sigma = "LM", int options = ComputeEigenvectors,
                                                   RealScalar tol = 0.0);

  const Matrix<Scalar, Dynamic, Dynamic> &eigenvectors() const {
    eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
    eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
    return m_eivec;
  }

  const Matrix<Scalar, Dynamic, 1> &eigenvalues() const {
    eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
    return m_eivalues;
  }

  Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const {
    eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
    eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
    return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
  }

  Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const {
    eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
    eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
    return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
  }

  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
    return m_info;
  }
 
  size_t getNbrConvergedEigenValues() const { return m_nbrConverged; }
 
  size_t getNbrIterations() const { return m_nbrIterations; }
 
 protected:
  Matrix<Scalar, Dynamic, Dynamic> m_eivec;
  Matrix<Scalar, Dynamic, 1> m_eivalues;
  ComputationInfo m_info;
  bool m_isInitialized;
  bool m_eigenvectorsOk;
 
  size_t m_nbrConverged;
  size_t m_nbrIterations;
};
 
template <typename MatrixType, typename MatrixSolver, bool BisSPD>
ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> &
ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>::compute(const MatrixType &A,
                                                                                   Index nbrEigenvalues,
                                                                                   std::string eigs_sigma, int options,
                                                                                   RealScalar tol) {
  MatrixType B(0, 0);
  compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
 
  return *this;
}
 
template <typename MatrixType, typename MatrixSolver, bool BisSPD>
ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> &
ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>::compute(const MatrixType &A,
                                                                                   const MatrixType &B,
                                                                                   Index nbrEigenvalues,
                                                                                   std::string eigs_sigma, int options,
                                                                                   RealScalar tol) {
  eigen_assert(A.cols() == A.rows());
  eigen_assert(B.cols() == B.rows());
  eigen_assert(B.rows() == 0 || A.cols() == B.rows());
  eigen_assert((options & ~(EigVecMask | GenEigMask)) == 0 && (options & EigVecMask) != EigVecMask &&
               "invalid option parameter");
 
  bool isBempty = (B.rows() == 0) || (B.cols() == 0);
 
  // For clarity, all parameters match their ARPACK name
  //
  // Always 0 on the first call
  //
  int ido = 0;
 
  int n = (int)A.cols();
 
  // User options: "LA", "SA", "SM", "LM", "BE"
  //
  char whch[3] = "LM";
 
  // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 }
  //
  RealScalar sigma = 0.0;
 
  if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) {
    eigs_sigma[0] = toupper(eigs_sigma[0]);
    eigs_sigma[1] = toupper(eigs_sigma[1]);
 
    // In the following special case we're going to invert the problem, since solving
    // for larger magnitude is much much faster
    // i.e., if 'SM' is specified, we're going to really use 'LM', the default
    //
    if (eigs_sigma.substr(0, 2) != "SM") {
      whch[0] = eigs_sigma[0];
      whch[1] = eigs_sigma[1];
    }
  } else {
    eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!");
 
    // If it's not scalar values, then the user may be explicitly
    // specifying the sigma value to cluster the evs around
    //
    sigma = atof(eigs_sigma.c_str());
 
    // If atof fails, it returns 0.0, which is a fine default
    //
  }
 
  // "I" means normal eigenvalue problem, "G" means generalized
  //
  char bmat[2] = "I";
  if (eigs_sigma.substr(0, 2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
    bmat[0] = 'G';
 
  // Now we determine the mode to use
  //
  int mode = (bmat[0] == 'G') + 1;
  if (eigs_sigma.substr(0, 2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))) {
    // We're going to use shift-and-invert mode, and basically find
    // the largest eigenvalues of the inverse operator
    //
    mode = 3;
  }
 
  // The user-specified number of eigenvalues/vectors to compute
  //
  int nev = (int)nbrEigenvalues;
 
  // Allocate space for ARPACK to store the residual
  //
  Scalar *resid = new Scalar[n];
 
  // Number of Lanczos vectors, must satisfy nev < ncv <= n
  // Note that this indicates that nev != n, and we cannot compute
  // all eigenvalues of a mtrix
  //
  int ncv = std::min(std::max(2 * nev, 20), n);
 
  // The working n x ncv matrix, also store the final eigenvectors (if computed)
  //
  Scalar *v = new Scalar[n * ncv];
  int ldv = n;
 
  // Working space
  //
  Scalar *workd = new Scalar[3 * n];
  int lworkl = ncv * ncv + 8 * ncv;  // Must be at least this length
  Scalar *workl = new Scalar[lworkl];
 
  int *iparam = new int[11];
  iparam[0] = 1;  // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it
  iparam[2] = std::max(300, (int)std::ceil(2 * n / std::max(ncv, 1)));
  iparam[6] = mode;  // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert
 
  // Used during reverse communicate to notify where arrays start
  //
  int *ipntr = new int[11];
 
  // Error codes are returned in here, initial value of 0 indicates a random initial
  // residual vector is used, any other values means resid contains the initial residual
  // vector, possibly from a previous run
  //
  int info = 0;
 
  Scalar scale = 1.0;
  // if (!isBempty)
  //{
  // Scalar scale = B.norm() / std::sqrt(n);
  // scale = std::pow(2, std::floor(std::log(scale+1)));
  // for (size_t i=0; i<(size_t)B.outerSize(); i++)
  //     for (typename MatrixType::InnerIterator it(B, i); it; ++it)
  //         it.valueRef() /= scale;
  // }
 
  MatrixSolver OP;
  if (mode == 1 || mode == 2) {
    if (!isBempty) OP.compute(B);
  } else if (mode == 3) {
    if (sigma == 0.0) {
      OP.compute(A);
    } else {
      // Note: We will never enter here because sigma must be 0.0
      //
      if (isBempty) {
        MatrixType AminusSigmaB(A);
        for (Index i = 0; i < A.rows(); ++i) AminusSigmaB.coeffRef(i, i) -= sigma;
 
        OP.compute(AminusSigmaB);
      } else {
        MatrixType AminusSigmaB = A - sigma * B;
        OP.compute(AminusSigmaB);
      }
    }
  }
 
  if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() != Success) {
    m_info = OP.info() delete[] v;
    delete[] iparam;
    delete[] ipntr;
    delete[] workd;
    delete[] workl;
    delete[] resid;
    m_isInitialized = false;
    return *this;
  }
 
  do {
    internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid, &ncv, v, &ldv, iparam,
                                                        ipntr, workd, workl, &lworkl, &info);
 
    if (ido == -1 || ido == 1) {
      Scalar *in = workd + ipntr[0] - 1;
      Scalar *out = workd + ipntr[1] - 1;
 
      if (ido == 1 && mode != 2) {
        Scalar *out2 = workd + ipntr[2] - 1;
        if (isBempty || mode == 1)
          Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
        else
          Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
 
        in = workd + ipntr[2] - 1;
      }
 
      if (mode == 1) {
        if (isBempty) {
          // OP = A
          //
          Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n);
        } else {
          // OP = L^{-1}AL^{-T}
          //
          internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
        }
      } else if (mode == 2) {
        if (ido == 1) Matrix<Scalar, Dynamic, 1>::Map(in, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n);
 
        // OP = B^{-1} A
        //
        Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
      } else if (mode == 3) {
        // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I)
        // The B * in is already computed and stored at in if ido == 1
        //
        if (ido == 1 || isBempty)
          Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
        else
          Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
      }
    } else if (ido == 2) {
      Scalar *in = workd + ipntr[0] - 1;
      Scalar *out = workd + ipntr[1] - 1;
 
      if (isBempty || mode == 1)
        Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
      else
        Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
    }
  } while (ido != 99);
 
  if (info == 1)
    m_info = NoConvergence;
  else if (info == 3)
    m_info = NumericalIssue;
  else if (info < 0)
    m_info = InvalidInput;
  else if (info != 0)
    eigen_assert(false && "Unknown ARPACK return value!");
  else {
    // Do we compute eigenvectors or not?
    //
    int rvec = (options & ComputeEigenvectors) == ComputeEigenvectors;
 
    // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK))
    //
    char howmny[2] = "A";
 
    // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK)
    //
    int *select = new int[ncv];
 
    // Final eigenvalues
    //
    m_eivalues.resize(nev, 1);
 
    internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv, &sigma, bmat,
                                                        &n, whch, &nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr,
                                                        workd, workl, &lworkl, &info);
 
    if (info == -14)
      m_info = NoConvergence;
    else if (info != 0)
      m_info = InvalidInput;
    else {
      if (rvec) {
        m_eivec.resize(A.rows(), nev);
        for (int i = 0; i < nev; i++)
          for (int j = 0; j < n; j++) m_eivec(j, i) = v[i * n + j] / scale;
 
        if (mode == 1 && !isBempty && BisSPD)
          internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
 
        m_eigenvectorsOk = true;
      }
 
      m_nbrIterations = iparam[2];
      m_nbrConverged = iparam[4];
 
      m_info = Success;
    }
 
    delete[] select;
  }
 
  delete[] v;
  delete[] iparam;
  delete[] ipntr;
  delete[] workd;
  delete[] workl;
  delete[] resid;
 
  m_isInitialized = (m_info == Success);
 
  return *this;
}
 
// Single precision
//
extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which, int *nev, float *tol, float *resid, int *ncv,
                        float *v, int *ldv, int *iparam, int *ipntr, float *workd, float *workl, int *lworkl,
                        int *info);
 
extern "C" void sseupd_(int *rvec, char *All, int *select, float *d, float *z, int *ldz, float *sigma, char *bmat,
                        int *n, char *which, int *nev, float *tol, float *resid, int *ncv, float *v, int *ldv,
                        int *iparam, int *ipntr, float *workd, float *workl, int *lworkl, int *ierr);
 
// Double precision
//
extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv,
                        double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl,
                        int *info);
 
extern "C" void dseupd_(int *rvec, char *All, int *select, double *d, double *z, int *ldz, double *sigma, char *bmat,
                        int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv,
                        int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *ierr);
 
namespace internal {
 
template <typename Scalar, typename RealScalar>
struct arpack_wrapper {
  static inline void saupd(int *ido, char *bmat, int *n, char *which, int *nev, RealScalar *tol, Scalar *resid,
                           int *ncv, Scalar *v, int *ldv, int *iparam, int *ipntr, Scalar *workd, Scalar *workl,
                           int *lworkl, int *info) {
    EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
  }
 
  static inline void seupd(int *rvec, char *All, int *select, Scalar *d, Scalar *z, int *ldz, RealScalar *sigma,
                           char *bmat, int *n, char *which, int *nev, RealScalar *tol, Scalar *resid, int *ncv,
                           Scalar *v, int *ldv, int *iparam, int *ipntr, Scalar *workd, Scalar *workl, int *lworkl,
                           int *ierr) {
    EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
  }
};
 
template <>
struct arpack_wrapper<float, float> {
  static inline void saupd(int *ido, char *bmat, int *n, char *which, int *nev, float *tol, float *resid, int *ncv,
                           float *v, int *ldv, int *iparam, int *ipntr, float *workd, float *workl, int *lworkl,
                           int *info) {
    ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
  }
 
  static inline void seupd(int *rvec, char *All, int *select, float *d, float *z, int *ldz, float *sigma, char *bmat,
                           int *n, char *which, int *nev, float *tol, float *resid, int *ncv, float *v, int *ldv,
                           int *iparam, int *ipntr, float *workd, float *workl, int *lworkl, int *ierr) {
    sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd,
            workl, lworkl, ierr);
  }
};
 
template <>
struct arpack_wrapper<double, double> {
  static inline void saupd(int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv,
                           double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl,
                           int *info) {
    dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
  }
 
  static inline void seupd(int *rvec, char *All, int *select, double *d, double *z, int *ldz, double *sigma, char *bmat,
                           int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv,
                           int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *ierr) {
    dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd,
            workl, lworkl, ierr);
  }
};
 
template <typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD>
struct OP {
  static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out);
  static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs);
};
 
template <typename MatrixSolver, typename MatrixType, typename Scalar>
struct OP<MatrixSolver, MatrixType, Scalar, true> {
  static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) {
    // OP = L^{-1} A L^{-T}  (B = LL^T)
    //
    // First solve L^T out = in
    //
    Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixU().solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
    Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationPinv() * Matrix<Scalar, Dynamic, 1>::Map(out, n);
 
    // Then compute out = A out
    //
    Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(out, n);
 
    // Then solve L out = out
    //
    Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationP() * Matrix<Scalar, Dynamic, 1>::Map(out, n);
    Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixL().solve(Matrix<Scalar, Dynamic, 1>::Map(out, n));
  }
 
  static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) {
    // Solve L^T out = in
    //
    Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) =
        OP.matrixU().solve(Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k));
    Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) =
        OP.permutationPinv() * Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k);
  }
};
 
template <typename MatrixSolver, typename MatrixType, typename Scalar>
struct OP<MatrixSolver, MatrixType, Scalar, false> {
  static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) {
    eigen_assert(false && "Should never be in here...");
  }
 
  static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) {
    eigen_assert(false && "Should never be in here...");
  }
};
 
}  // end namespace internal
 
}  // end namespace Eigen
 
#endif  // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H