10 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
11 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
19 using Eigen::ComputationInfo;
20 using Eigen::ComputeEigenvectors;
21 using Eigen::NoConvergence;
22 using Eigen::NumericalIssue;
23 #ifndef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
24 using Eigen::InvalidInput;
27 using Eigen::NumTraits;
30 namespace arpack_internal
35 template<
class LMatrixType,
class RMatrixType,
class MatrixOperation,
bool BisSPD=false>
42 typedef typename LMatrixType::Scalar
Scalar;
43 typedef typename LMatrixType::Index
Index;
58 typedef typename Eigen::internal::plain_col_type<LMatrixType, RealScalar>::type
RealVectorType;
94 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
95 int parameters=ComputeEigenvectors,
RealScalar tol=0.0) :
99 compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
124 int parameters=ComputeEigenvectors,
RealScalar tol=0.0) :
128 compute(A, nbrEigenvalues, eigs_sigma, parameters, tol);
155 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
156 int parameters=ComputeEigenvectors,
RealScalar tol=0.0);
181 std::string eigs_sigma=
"LM",
182 int parameters=ComputeEigenvectors,
RealScalar tol=0.0);
202 eigen_assert(
m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
203 eigen_assert(
m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
221 eigen_assert(
m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
241 eigen_assert(
m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
242 eigen_assert(
m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
262 eigen_assert(
m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
263 eigen_assert(
m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
273 eigen_assert(
m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
299 template<
typename LMatrixType,
typename RMatrixType,
class MatrixOperation,
bool BisSPD>
303 std::string eigs_sigma,
int parameters,
RealScalar tol)
306 return compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
309 template<
typename LMatrixType,
typename RMatrixType,
class MatrixOperation,
bool BisSPD>
313 std::string eigs_sigma,
int parameters,
RealScalar tol)
315 eigen_assert(B.cols() == B.rows());
316 eigen_assert(B.rows() == 0 || A.rows() == B.rows() || A.cols() == B.cols());
317 eigen_assert((parameters &~ (Eigen::EigVecMask | Eigen::GenEigMask)) == 0
318 && (parameters & Eigen::EigVecMask) != Eigen::EigVecMask
319 &&
"invalid option parameter");
321 bool isBempty = (B.rows() == 0) || (B.cols() == 0);
326 int n = (int)A.rows();
333 if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
335 eigs_sigma[0] = toupper(eigs_sigma[0]);
336 eigs_sigma[1] = toupper(eigs_sigma[1]);
340 if (eigs_sigma.substr(0,2) !=
"SM")
342 whch[0] = eigs_sigma[0];
343 whch[1] = eigs_sigma[1];
348 eigen_assert(
false &&
"Specifying clustered eigenvalues is not yet supported!");
351 sigma = atof(eigs_sigma.c_str());
358 if (eigs_sigma.substr(0,2) ==
"SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
363 int mode = (bmat[0] ==
'G') + 1;
364 if (eigs_sigma.substr(0,2) ==
"SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
373 int nev = (int)nbrEigenvalues;
379 int ncv = std::min(std::max(2*nev, 20), n);
387 int lworkl = ncv*ncv+8*ncv;
390 int *iparam=
new int[11];
392 iparam[2] = std::max(300, (
int)std::ceil(2*n/std::max(ncv,1)));
396 int *ipntr =
new int[11];
416 MatrixOperation op(A);
424 &ncv, v, &ldv, iparam, ipntr, workd, workl,
426 if (ido == -1 || ido == 1)
428 Scalar *in = workd + ipntr[0] - 1;
429 Scalar *out = workd + ipntr[1] - 1;
431 if (ido == 1 && mode != 2)
433 Scalar *out2 = workd + ipntr[2] - 1;
434 if (isBempty || mode == 1)
435 Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
437 Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
438 in = workd + ipntr[2] - 1;
446 Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
457 Matrix<Scalar, Dynamic, 1>::Map(in, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
459 Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
465 if (ido == 1 || isBempty)
466 Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
468 Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
473 Scalar *in = workd + ipntr[0] - 1;
474 Scalar *out = workd + ipntr[1] - 1;
476 if (isBempty || mode == 1)
477 Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
479 Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
490 m_info = NoConvergence;
495 m_info = NumericalIssue;
500 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
501 m_info = NoConvergence;
503 m_info = InvalidInput;
507 eigen_assert(
false &&
"Unknown ARPACK return value!");
511 int rvec = (parameters & ComputeEigenvectors) == ComputeEigenvectors;
514 char howmny[2] =
"A";
517 int *select =
new int[ncv];
520 m_eivalues.resize(nev, 1);
523 &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
524 v, &ldv, iparam, ipntr, workd, workl, &lworkl, &cinfo);
527 m_info = NoConvergence;
529 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
530 m_info = NoConvergence;
532 m_info = InvalidInput;
538 m_eivec.resize(A.rows(), nev);
539 for (
int i=0; i<nev; i++)
540 for (
int j=0; j<n; j++)
541 m_eivec(j,i) = v[i*n+j] / scale;
546 m_eigenvectorsOk =
true;
549 m_nbrIterations = iparam[2];
550 m_nbrConverged = iparam[4];
562 m_isInitialized =
true;
568 extern "C" void ssaupd_(
int *ido,
char *bmat,
int *n,
char *which,
569 int *nev,
float *tol,
float *resid,
int *ncv,
570 float *v,
int *ldv,
int *iparam,
int *ipntr,
571 float *workd,
float *workl,
int *lworkl,
574 extern "C" void sseupd_(
int *rvec,
char *All,
int *select,
float *d,
575 float *z,
int *ldz,
float *sigma,
576 char *bmat,
int *n,
char *which,
int *nev,
577 float *tol,
float *resid,
int *ncv,
float *v,
578 int *ldv,
int *iparam,
int *ipntr,
float *workd,
579 float *workl,
int *lworkl,
int *ierr);
583 extern "C" void dsaupd_(
int *ido,
char *bmat,
int *n,
char *which,
584 int *nev,
double *tol,
double *resid,
int *ncv,
585 double *v,
int *ldv,
int *iparam,
int *ipntr,
586 double *workd,
double *workl,
int *lworkl,
589 extern "C" void dseupd_(
int *rvec,
char *All,
int *select,
double *d,
590 double *z,
int *ldz,
double *sigma,
591 char *bmat,
int *n,
char *which,
int *nev,
592 double *tol,
double *resid,
int *ncv,
double *v,
593 int *ldv,
int *iparam,
int *ipntr,
double *workd,
594 double *workl,
int *lworkl,
int *ierr);
596 namespace arpack_internal {
598 template<
typename Scalar,
typename RealScalar>
struct arpack_wrapper
600 static inline void saupd(
int *ido,
char *bmat,
int *n,
char *which,
601 int *nev, RealScalar *tol, Scalar *resid,
int *ncv,
602 Scalar *v,
int *ldv,
int *iparam,
int *ipntr,
603 Scalar *workd, Scalar *workl,
int *lworkl,
int *info);
605 static inline void seupd(
int *rvec,
char *All,
int *select, Scalar *d,
606 Scalar *z,
int *ldz, RealScalar *sigma,
607 char *bmat,
int *n,
char *which,
int *nev,
608 RealScalar *tol, Scalar *resid,
int *ncv, Scalar *v,
609 int *ldv,
int *iparam,
int *ipntr, Scalar *workd,
610 Scalar *workl,
int *lworkl,
int *ierr);
615 static inline void saupd(
int *ido,
char *bmat,
int *n,
char *which,
616 int *nev,
float *tol,
float *resid,
int *ncv,
617 float *v,
int *ldv,
int *iparam,
int *ipntr,
618 float *workd,
float *workl,
int *lworkl,
int *info)
620 ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
623 static inline void seupd(
int *rvec,
char *All,
int *select,
float *d,
624 float* z,
int* ldz,
float *sigma,
625 char *bmat,
int *n,
char *which,
int *nev,
626 float *tol,
float *resid,
int *ncv,
float *v,
627 int *ldv,
int *iparam,
int *ipntr,
float *workd,
628 float *workl,
int *lworkl,
int *ierr)
630 sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
631 workd, workl, lworkl, ierr);
637 static inline void saupd(
int *ido,
char *bmat,
int *n,
char *which,
638 int *nev,
double *tol,
double *resid,
int *ncv,
639 double *v,
int *ldv,
int *iparam,
int *ipntr,
640 double *workd,
double *workl,
int *lworkl,
int *info)
642 dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
645 static inline void seupd(
int *rvec,
char *All,
int *select,
double *d,
646 double* z,
int* ldz,
double *sigma,
647 char *bmat,
int *n,
char *which,
int *nev,
648 double *tol,
double *resid,
int *ncv,
double *v,
649 int *ldv,
int *iparam,
int *ipntr,
double *workd,
650 double *workl,
int *lworkl,
int *ierr)
652 dseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
653 workd, workl, lworkl, ierr);
659 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
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)
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)
ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType &A, const RMatrixType &B, Index nbrEigenvalues, std::string eigs_sigma="LM", int parameters=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes generalized eigenvalues of given matrix with respect to another matrix...
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)
ArpackGeneralizedSelfAdjointEigenSolver & compute(const LMatrixType &A, const RMatrixType &B, Index nbrEigenvalues, std::string eigs_sigma="LM", int parameters=ComputeEigenvectors, RealScalar tol=0.0)
Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library...
static 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)
static 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)
Matrix< Scalar, Dynamic, Dynamic > operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Eigen::internal::plain_col_type< LMatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
static 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)
static 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)
size_t getNbrIterations() const
NumTraits< Scalar >::Real RealScalar
Real scalar type for MatrixType.
ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType &A, Index nbrEigenvalues, std::string eigs_sigma="LM", int parameters=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes eigenvalues of given matrix.
Matrix< Scalar, Dynamic, Dynamic > operatorSqrt() const
Computes the positive-definite square root of the matrix.
const Matrix< Scalar, Dynamic, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
static 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)
ComputationInfo info() const
Reports whether previous computation was successful.
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)
LMatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
static 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)
ArpackGeneralizedSelfAdjointEigenSolver()
Default constructor.
Matrix< Scalar, Dynamic, 1 > m_eivalues
Matrix< Scalar, Dynamic, Dynamic > m_eivec
size_t getNbrConvergedEigenValues() const