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>
58 m_isInitialized(false),
59 m_eigenvectorsOk(false),
87 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
91 m_isInitialized(false),
92 m_eigenvectorsOk(false),
122 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
126 m_isInitialized(false),
127 m_eigenvectorsOk(false),
159 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
185 Index nbrEigenvalues, std::string eigs_sigma=
"LM",
210 eigen_assert(m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
211 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
232 eigen_assert(m_isInitialized &&
"ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
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();
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.");
297 {
return m_nbrConverged; }
300 {
return m_nbrIterations; }
317 template<
typename MatrixType,
typename MatrixSolver,
bool BisSPD>
324 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
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];
458 if (mode == 1 || mode == 2)
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,
675 static inline void seupd(
int *rvec,
char *All,
int *select,
Scalar *d,
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>
735 static inline void project(MatrixSolver &OP,
int n,
int k,
Scalar *vecs);
738 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar>
770 template<
typename MatrixSolver,
typename MatrixType,
typename Scalar>
789 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H 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)
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)
static Point2 project(const Pose3 &pose, const Unit3 &pointAtInfinity, const Cal3_S2::shared_ptr &cal)
Matrix< Scalar, Dynamic, Dynamic > m_eivec
size_t getNbrIterations() const
std::ofstream out("Result.txt")
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...
Namespace containing all symbols from the Eigen library.
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType &A, Index nbrEigenvalues, std::string eigs_sigma="LM", int options=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes eigenvalues of given matrix.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
static void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
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)
Matrix< Scalar, Dynamic, Dynamic > operatorSqrt() const
Computes the positive-definite square root of the matrix.
Matrix< Scalar, Dynamic, Dynamic > operatorInverseSqrt() const
Computes the inverse square root of the matrix.
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
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, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
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, RealScalar *tol, Scalar *resid, int *ncv, Scalar *v, int *ldv, int *iparam, int *ipntr, Scalar *workd, Scalar *workl, int *lworkl, int *info)
ArpackGeneralizedSelfAdjointEigenSolver()
Default constructor.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Array< int, Dynamic, 1 > v
NumTraits< Scalar >::Real RealScalar
Real scalar type for MatrixType.
static void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
NumTraits< Scalar >::Real RealScalar
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
size_t getNbrConvergedEigenValues() const
Matrix< Scalar, Dynamic, 1 > m_eivalues
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 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)
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)
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 ConstMapType Map(const Scalar *data)
ComputationInfo info() const
Reports whether previous computation was successful.
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 const double sigma
static const DiscreteKey mode(modeKey, 2)
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)
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
conditional< is_same< typename traits< ExpressionType >::XprKind, MatrixXpr >::value, MatrixColType, ArrayColType >::type type
EIGEN_DEVICE_FUNC const CeilReturnType ceil() const
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)
static void project(MatrixSolver &OP, int n, int k, Scalar *vecs)