fix eigenvectors computations :)

This commit is contained in:
Gael Guennebaud 2008-06-03 18:03:55 +00:00
parent 915587d03d
commit a0cff1a295
3 changed files with 8 additions and 5 deletions

View File

@ -281,7 +281,7 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
{ {
if ((Rhs::Flags&Diagonal)==Diagonal) if ((Rhs::Flags&Diagonal)==Diagonal)
{ {
assert(_LhsNested::Flags&RowMajorBit==0); assert((_LhsNested::Flags&RowMajorBit)==0);
return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col))); return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
} }
else else

View File

@ -45,7 +45,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex; typedef std::complex<RealScalar> Complex;
// typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> EigenvalueType;
typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
@ -120,6 +119,9 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
m_eivalues.resize(n,1); m_eivalues.resize(n,1);
m_eivec = matrix; m_eivec = matrix;
// FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?
// the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
// (same for diag and subdiag)
Tridiagonalization<MatrixType> tridiag(m_eivec); Tridiagonalization<MatrixType> tridiag(m_eivec);
RealVectorType& diag = m_eivalues; RealVectorType& diag = m_eivalues;
RealVectorTypeX subdiag(n-1); RealVectorTypeX subdiag(n-1);
@ -128,9 +130,6 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
if (computeEigenvectors) if (computeEigenvectors)
m_eivec = tridiag.matrixQ(); m_eivec = tridiag.matrixQ();
RealVectorTypeX gc(n);
RealVectorTypeX gs(n);
int end = n-1; int end = n-1;
int start = 0; int start = 0;
while (end>0) while (end>0)

View File

@ -197,6 +197,10 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
matA.col(i).coeffRef(i+1) = beta; matA.col(i).coeffRef(i+1) = beta;
hCoeffs.coeffRef(i) = (beta - v0) / beta; hCoeffs.coeffRef(i) = (beta - v0) / beta;
} }
else
{
hCoeffs.coeffRef(n-2) = 0;
}
} }
/** reconstructs and returns the matrix Q */ /** reconstructs and returns the matrix Q */