GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-17 03:43:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

make eigen2 eigensolver test pass

Browse Source
This commit is contained in:
Benoit Jacob 2011-01-31 08:36:14 -05:00
parent e2642ed620
commit 374deaed5f
2 changed files with 38 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -29,6 +29,9 @@
#include "./EigenvaluesCommon.h"
#include "./Tridiagonalization.h"
template<typename _MatrixType>
class GeneralizedSelfAdjointEigenSolver;
/** \eigenvalues_module \ingroup Eigenvalues_Module
*
*
@ -310,6 +313,36 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*/
static const int m_maxIterations = 30;
#ifdef EIGEN2_SUPPORT
SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
: m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()),
m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
m_isInitialized(false)
{
compute(matrix, computeEigenvectors);
}
SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
: m_eivec(matA.cols(), matA.cols()),
m_eivalues(matA.cols()),
m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
m_isInitialized(false)
{
static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
}
void compute(const MatrixType& matrix, bool computeEigenvectors)
{
compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
}
void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
{
compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
}
#endif // EIGEN2_SUPPORT
protected:
MatrixType m_eivec;
RealVectorType m_eivalues;

10
test/eigen2/eigen2_eigensolver.cpp
View File

@ -75,7 +75,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
convert(gEvec, _evec);
// test gsl itself !
VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal().eval(), largerEps));
VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
// compare with eigen
VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
@ -86,7 +86,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
convert(gEval, _eval);
convert(gEvec, _evec);
// test GSL itself:
VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal().eval()), largerEps));
VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
// compare with eigen
// std::cerr << _eval.transpose() << "\n" << eiSymmGen.eigenvalues().transpose() << "\n\n";
@ -102,11 +102,11 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
#endif
VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval(), largerEps));
eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
// generalized eigen problem Ax = lBx
VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal().eval()), largerEps));
symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
MatrixType sqrtSymmA = eiSymm.operatorSqrt();
VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);
@ -141,7 +141,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
EigenSolver<MatrixType> ei1(a);
VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval());
ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
}