10#ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
11#define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
13#include "../../../../Eigen/Dense"
18 template<
typename Scalar,
typename RealScalar>
struct arpack_wrapper;
19 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar,
bool BisSPD>
struct OP;
24template<
typename MatrixType,
typename MatrixSolver=SimplicialLLT<MatrixType>,
bool BisSPD=false>
25class ArpackGeneralizedSelfAdjointEigenSolver
31 typedef typename MatrixType::Scalar Scalar;
32 typedef typename MatrixType::Index
Index;
40 typedef typename NumTraits<Scalar>::Real RealScalar;
47 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
55 ArpackGeneralizedSelfAdjointEigenSolver()
58 m_isInitialized(false),
59 m_eigenvectorsOk(false),
86 ArpackGeneralizedSelfAdjointEigenSolver(
const MatrixType& A,
const MatrixType& B,
87 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
91 m_isInitialized(false),
92 m_eigenvectorsOk(false),
96 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
121 ArpackGeneralizedSelfAdjointEigenSolver(
const MatrixType& A,
122 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
126 m_isInitialized(false),
127 m_eigenvectorsOk(false),
131 compute(A, nbrEigenvalues, eigs_sigma, options, tol);
158 ArpackGeneralizedSelfAdjointEigenSolver& compute(
const MatrixType& A,
const MatrixType& B,
159 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
184 ArpackGeneralizedSelfAdjointEigenSolver& compute(
const MatrixType& A,
185 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
208 const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors()
const
210 eigen_assert(m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
211 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
230 const Matrix<Scalar, Dynamic, 1>& eigenvalues()
const
232 eigen_assert(m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
254 Matrix<Scalar, Dynamic, Dynamic> operatorSqrt()
const
256 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
257 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
258 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
279 Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt()
const
281 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
282 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
283 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
292 eigen_assert(m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
296 size_t getNbrConvergedEigenValues()
const
297 {
return m_nbrConverged; }
299 size_t getNbrIterations()
const
300 {
return m_nbrIterations; }
303 Matrix<Scalar, Dynamic, Dynamic> m_eivec;
304 Matrix<Scalar, Dynamic, 1> m_eivalues;
306 bool m_isInitialized;
307 bool m_eigenvectorsOk;
309 size_t m_nbrConverged;
310 size_t m_nbrIterations;
317template<
typename MatrixType,
typename MatrixSolver,
bool BisSPD>
318ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>&
319 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>
320::compute(
const MatrixType& A,
Index nbrEigenvalues,
321 std::string eigs_sigma,
int options, RealScalar tol)
324 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
330template<
typename MatrixType,
typename MatrixSolver,
bool BisSPD>
331ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>&
332 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>
333::compute(
const MatrixType& A,
const MatrixType& B,
Index nbrEigenvalues,
334 std::string eigs_sigma,
int options, RealScalar tol)
336 eigen_assert(A.cols() == A.rows());
337 eigen_assert(B.cols() == B.rows());
338 eigen_assert(B.rows() == 0 || A.cols() == B.rows());
339 eigen_assert((options &~ (EigVecMask | GenEigMask)) == 0
340 && (options & EigVecMask) != EigVecMask
341 &&
"invalid option parameter");
343 bool isBempty = (B.rows() == 0) || (B.cols() == 0);
351 int n = (int)A.cols();
359 RealScalar sigma = 0.0;
361 if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
363 eigs_sigma[0] = toupper(eigs_sigma[0]);
364 eigs_sigma[1] = toupper(eigs_sigma[1]);
370 if (eigs_sigma.substr(0,2) !=
"SM")
372 whch[0] = eigs_sigma[0];
373 whch[1] = eigs_sigma[1];
378 eigen_assert(
false &&
"Specifying clustered eigenvalues is not yet supported!");
383 sigma = atof(eigs_sigma.c_str());
392 if (eigs_sigma.substr(0,2) ==
"SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
397 int mode = (bmat[0] ==
'G') + 1;
398 if (eigs_sigma.substr(0,2) ==
"SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
408 int nev = (int)nbrEigenvalues;
412 Scalar *resid =
new Scalar[n];
418 int ncv = std::min(std::max(2*nev, 20), n);
422 Scalar *v =
new Scalar[n*ncv];
427 Scalar *workd =
new Scalar[3*n];
428 int lworkl = ncv*ncv+8*ncv;
429 Scalar *workl =
new Scalar[lworkl];
431 int *iparam=
new int[11];
433 iparam[2] = std::max(300, (
int)std::ceil(2*n/std::max(ncv,1)));
438 int *ipntr =
new int[11];
458 if (mode == 1 || mode == 2)
475 MatrixType AminusSigmaB(A);
476 for (
Index i=0; i<A.rows(); ++i)
477 AminusSigmaB.coeffRef(i,i) -= sigma;
479 OP.compute(AminusSigmaB);
483 MatrixType AminusSigmaB = A - sigma * B;
484 OP.compute(AminusSigmaB);
489 if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() !=
Success)
490 std::cout <<
"Error factoring matrix" << std::endl;
494 internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid,
495 &ncv, v, &ldv, iparam, ipntr, workd, workl,
498 if (ido == -1 || ido == 1)
500 Scalar *in = workd + ipntr[0] - 1;
501 Scalar *out = workd + ipntr[1] - 1;
503 if (ido == 1 && mode != 2)
505 Scalar *out2 = workd + ipntr[2] - 1;
506 if (isBempty || mode == 1)
511 in = workd + ipntr[2] - 1;
526 internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
543 if (ido == 1 || isBempty)
551 Scalar *in = workd + ipntr[0] - 1;
552 Scalar *out = workd + ipntr[1] - 1;
554 if (isBempty || mode == 1)
568 eigen_assert(
false &&
"Unknown ARPACK return value!");
577 char howmny[2] =
"A";
581 int *select =
new int[ncv];
585 m_eivalues.resize(nev, 1);
587 internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv,
588 &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
589 v, &ldv, iparam, ipntr, workd, workl, &lworkl, &info);
599 m_eivec.resize(A.rows(), nev);
600 for (
int i=0; i<nev; i++)
601 for (
int j=0; j<n; j++)
602 m_eivec(j,i) = v[i*n+j] / scale;
604 if (mode == 1 && !isBempty && BisSPD)
605 internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
607 m_eigenvectorsOk =
true;
610 m_nbrIterations = iparam[2];
611 m_nbrConverged = iparam[4];
626 m_isInitialized =
true;
634extern "C" void ssaupd_(
int *ido,
char *bmat,
int *n,
char *which,
635 int *nev,
float *tol,
float *resid,
int *ncv,
636 float *v,
int *ldv,
int *iparam,
int *ipntr,
637 float *workd,
float *workl,
int *lworkl,
640extern "C" void sseupd_(
int *rvec,
char *All,
int *select,
float *d,
641 float *z,
int *ldz,
float *sigma,
642 char *bmat,
int *n,
char *which,
int *nev,
643 float *tol,
float *resid,
int *ncv,
float *v,
644 int *ldv,
int *iparam,
int *ipntr,
float *workd,
645 float *workl,
int *lworkl,
int *ierr);
649extern "C" void dsaupd_(
int *ido,
char *bmat,
int *n,
char *which,
650 int *nev,
double *tol,
double *resid,
int *ncv,
651 double *v,
int *ldv,
int *iparam,
int *ipntr,
652 double *workd,
double *workl,
int *lworkl,
655extern "C" void dseupd_(
int *rvec,
char *All,
int *select,
double *d,
656 double *z,
int *ldz,
double *sigma,
657 char *bmat,
int *n,
char *which,
int *nev,
658 double *tol,
double *resid,
int *ncv,
double *v,
659 int *ldv,
int *iparam,
int *ipntr,
double *workd,
660 double *workl,
int *lworkl,
int *ierr);
665template<
typename Scalar,
typename RealScalar>
struct arpack_wrapper
667 static inline void saupd(
int *ido,
char *bmat,
int *n,
char *which,
668 int *nev, RealScalar *tol, Scalar *resid,
int *ncv,
669 Scalar *v,
int *ldv,
int *iparam,
int *ipntr,
670 Scalar *workd, Scalar *workl,
int *lworkl,
int *info)
672 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
675 static inline void seupd(
int *rvec,
char *All,
int *select, Scalar *d,
676 Scalar *z,
int *ldz, RealScalar *sigma,
677 char *bmat,
int *n,
char *which,
int *nev,
678 RealScalar *tol, Scalar *resid,
int *ncv, Scalar *v,
679 int *ldv,
int *iparam,
int *ipntr, Scalar *workd,
680 Scalar *workl,
int *lworkl,
int *ierr)
682 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
686template <>
struct arpack_wrapper<float, float>
688 static inline void saupd(
int *ido,
char *bmat,
int *n,
char *which,
689 int *nev,
float *tol,
float *resid,
int *ncv,
690 float *v,
int *ldv,
int *iparam,
int *ipntr,
691 float *workd,
float *workl,
int *lworkl,
int *info)
693 ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
696 static inline void seupd(
int *rvec,
char *All,
int *select,
float *d,
697 float *z,
int *ldz,
float *sigma,
698 char *bmat,
int *n,
char *which,
int *nev,
699 float *tol,
float *resid,
int *ncv,
float *v,
700 int *ldv,
int *iparam,
int *ipntr,
float *workd,
701 float *workl,
int *lworkl,
int *ierr)
703 sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
704 workd, workl, lworkl, ierr);
708template <>
struct arpack_wrapper<double, double>
710 static inline void saupd(
int *ido,
char *bmat,
int *n,
char *which,
711 int *nev,
double *tol,
double *resid,
int *ncv,
712 double *v,
int *ldv,
int *iparam,
int *ipntr,
713 double *workd,
double *workl,
int *lworkl,
int *info)
715 dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
718 static inline void seupd(
int *rvec,
char *All,
int *select,
double *d,
719 double *z,
int *ldz,
double *sigma,
720 char *bmat,
int *n,
char *which,
int *nev,
721 double *tol,
double *resid,
int *ncv,
double *v,
722 int *ldv,
int *iparam,
int *ipntr,
double *workd,
723 double *workl,
int *lworkl,
int *ierr)
725 dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
726 workd, workl, lworkl, ierr);
731template<
typename MatrixSolver,
typename MatrixType,
typename Scalar,
bool BisSPD>
734 static inline void applyOP(MatrixSolver &OP,
const MatrixType &A,
int n, Scalar *in, Scalar *out);
735 static inline void project(MatrixSolver &OP,
int n,
int k, Scalar *vecs);
738template<
typename MatrixSolver,
typename MatrixType,
typename Scalar>
739struct OP<MatrixSolver, MatrixType, Scalar, true>
741 static inline void applyOP(MatrixSolver &OP,
const MatrixType &A,
int n, Scalar *in, Scalar *out)
760 static inline void project(MatrixSolver &OP,
int n,
int k, Scalar *vecs)
770template<
typename MatrixSolver,
typename MatrixType,
typename Scalar>
771struct OP<MatrixSolver, MatrixType, Scalar, false>
773 static inline void applyOP(MatrixSolver &OP,
const MatrixType &A,
int n, Scalar *in, Scalar *out)
775 eigen_assert(
false &&
"Should never be in here...");
778 static inline void project(MatrixSolver &OP,
int n,
int k, Scalar *vecs)
780 eigen_assert(
false &&
"Should never be in here...");
friend friend class Eigen::Map
Namespace containing all symbols from the Eigen library.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index