Go to the documentation of this file.
10 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
11 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
13 #include "../../../../Eigen/Dense"
19 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar,
bool BisSPD>
struct OP;
24 template<
typename MatrixType,
typename MatrixSolver=SimplicialLLT<MatrixType>,
bool BisSPD=false>
87 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
122 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
159 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
185 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
317 template<
typename MatrixType,
typename MatrixSolver,
bool BisSPD>
330 template<
typename MatrixType,
typename MatrixSolver,
bool BisSPD>
341 &&
"invalid option parameter");
343 bool isBempty = (
B.rows() == 0) || (
B.cols() == 0);
351 int n = (
int)
A.cols();
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;
428 int lworkl = ncv*ncv+8*ncv;
431 int *iparam=
new int[11];
438 int *ipntr =
new int[11];
477 AminusSigmaB.coeffRef(
i,
i) -=
sigma;
479 OP.compute(AminusSigmaB);
484 OP.compute(AminusSigmaB);
489 if (!(
mode == 1 && isBempty) && !(
mode == 2 && isBempty) &&
OP.info() !=
Success)
490 std::cout <<
"Error factoring matrix" << std::endl;
495 &ncv,
v, &ldv, iparam, ipntr, workd, workl,
498 if (ido == -1 || ido == 1)
500 Scalar *in = workd + ipntr[0] - 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;
543 if (ido == 1 || isBempty)
551 Scalar *in = workd + ipntr[0] - 1;
554 if (isBempty ||
mode == 1)
577 char howmny[2] =
"A";
581 int *select =
new int[ncv];
585 m_eivalues.resize(nev, 1);
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++)
604 if (
mode == 1 && !isBempty && BisSPD)
607 m_eigenvectorsOk =
true;
610 m_nbrIterations = iparam[2];
611 m_nbrConverged = iparam[4];
626 m_isInitialized =
true;
634 extern "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,
640 extern "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);
649 extern "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,
655 extern "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);
665 template<
typename Scalar,
typename RealScalar>
struct arpack_wrapper
667 static inline void saupd(
int *ido,
char *bmat,
int *
n,
char *which,
669 Scalar *
v,
int *ldv,
int *iparam,
int *ipntr,
677 char *bmat,
int *
n,
char *which,
int *nev,
679 int *ldv,
int *iparam,
int *ipntr,
Scalar *workd,
680 Scalar *workl,
int *lworkl,
int *ierr)
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);
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);
731 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar,
bool BisSPD>
738 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar>
770 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar>
789 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
Namespace containing all symbols from the Eigen library.
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)
static const double d[K][N]
ComputationInfo info() const
Reports whether previous computation was successful.
ArpackGeneralizedSelfAdjointEigenSolver()
Default constructor.
static void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
size_t getNbrConvergedEigenValues() const
static const double sigma
static void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
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)
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)
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
ArpackGeneralizedSelfAdjointEigenSolver & compute(const MatrixType &A, const MatrixType &B, Index nbrEigenvalues, std::string eigs_sigma="LM", int options=ComputeEigenvectors, RealScalar tol=0.0)
Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library.
Matrix< Scalar, Dynamic, 1 > m_eivalues
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
static void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
static const DiscreteKey mode(modeKey, 2)
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)
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
static void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
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)
ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType &A, Index nbrEigenvalues, std::string eigs_sigma="LM", int options=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes eigenvalues of given matrix.
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)
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Matrix< Scalar, Dynamic, Dynamic > operatorSqrt() const
Computes the positive-definite square root of the matrix.
std::ofstream out("Result.txt")
static ConstMapType Map(const Scalar *data)
Matrix< Scalar, Dynamic, Dynamic > operatorInverseSqrt() const
Computes the inverse square root of the matrix.
NumTraits< Scalar >::Real RealScalar
NumTraits< Scalar >::Real RealScalar
Real scalar type for MatrixType.
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)
Matrix< Scalar, Dynamic, Dynamic > m_eivec
const Matrix< Scalar, Dynamic, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Array< int, Dynamic, 1 > v
size_t getNbrIterations() const
const EIGEN_DEVICE_FUNC CeilReturnType ceil() const
ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType &A, const MatrixType &B, Index nbrEigenvalues, std::string eigs_sigma="LM", int options=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes generalized eigenvalues of given matrix with respect to another matrix.
static void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
static void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
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)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
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)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:01:48