mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-09 07:00:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Replace local variables by member variables in compute() methods.

Browse Source 
This is to avoid dynamic memory allocations in the compute() methods of
ComplexEigenSolver, EigenSolver, and SelfAdjointEigenSolver where possible.
As a result, Tridiagonalization::decomposeInPlace() is no longer used.
Biggest remaining issue is the allocation in HouseholderSequence::evalTo().
This commit is contained in:
Jitse Niesen 2010-05-24 17:43:06 +01:00
parent 68820fd4e8
commit 7a43a4408b
4 changed files with 101 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@ -3,6 +3,7 @@
// //
// Copyright (C) 2009 Claire Maurice // Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
// //
// Eigen is free software; you can redistribute it and/or // Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public // modify it under the terms of the GNU Lesser General Public
@ -99,7 +100,8 @@ template<typename _MatrixType> class ComplexEigenSolver
: m_eivec(), : m_eivec(),
m_eivalues(), m_eivalues(),
m_schur(), m_schur(),
m_isInitialized(false) m_isInitialized(false),
m_matX()
{} {}
/** \brief Default Constructor with memory preallocation /** \brief Default Constructor with memory preallocation
@ -112,7 +114,8 @@ template<typename _MatrixType> class ComplexEigenSolver
: m_eivec(size, size), : m_eivec(size, size),
m_eivalues(size), m_eivalues(size),
m_schur(size), m_schur(size),
m_isInitialized(false) m_isInitialized(false),
m_matX(size, size)
{} {}
/** \brief Constructor; computes eigendecomposition of given matrix. /** \brief Constructor; computes eigendecomposition of given matrix.
@ -125,7 +128,8 @@ template<typename _MatrixType> class ComplexEigenSolver
: m_eivec(matrix.rows(),matrix.cols()), : m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()), m_eivalues(matrix.cols()),
m_schur(matrix.rows()), m_schur(matrix.rows()),
m_isInitialized(false) m_isInitialized(false),
m_matX(matrix.rows(),matrix.cols())
{ {
compute(matrix); compute(matrix);
} }
@ -199,6 +203,7 @@ template<typename _MatrixType> class ComplexEigenSolver
EigenvalueType m_eivalues; EigenvalueType m_eivalues;
ComplexSchur<MatrixType> m_schur; ComplexSchur<MatrixType> m_schur;
bool m_isInitialized; bool m_isInitialized;
EigenvectorType m_matX;
}; };
@ -217,16 +222,16 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
// Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T. // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular. // The matrix X is unit triangular.
EigenvectorType X = EigenvectorType::Zero(n, n); m_matX = EigenvectorType::Zero(n, n);
for(int k=n-1 ; k>=0 ; k--) for(int k=n-1 ; k>=0 ; k--)
{ {
X.coeffRef(k,k) = ComplexScalar(1.0,0.0); m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
// Compute X(i,k) using the (i,k) entry of the equation X T = D X // Compute X(i,k) using the (i,k) entry of the equation X T = D X
for(int i=k-1 ; i>=0 ; i--) for(int i=k-1 ; i>=0 ; i--)
{ {
X.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
if(k-i-1>0) if(k-i-1>0)
X.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value(); m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
if(z==ComplexScalar(0)) if(z==ComplexScalar(0))
{ {
@ -234,12 +239,12 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
// Use a small value instead, to prevent division by zero. // Use a small value instead, to prevent division by zero.
ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
} }
X.coeffRef(i,k) = X.coeff(i,k) / z; m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
} }
} }
// Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) // Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
m_eivec = m_schur.matrixU() * X; m_eivec.noalias() = m_schur.matrixU() * m_matX;
// .. and normalize the eigenvectors // .. and normalize the eigenvectors
for(int k=0 ; k<n ; k++) for(int k=0 ; k<n ; k++)
{ {

119
Eigen/src/Eigenvalues/EigenSolver.h
View File

@ -120,7 +120,7 @@ template<typename _MatrixType> class EigenSolver
* *
* \sa compute() for an example. * \sa compute() for an example.
*/ */
EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {} EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
/** \brief Default Constructor with memory preallocation /** \brief Default Constructor with memory preallocation
* *
@ -131,7 +131,11 @@ template<typename _MatrixType> class EigenSolver
EigenSolver(int size) EigenSolver(int size)
: m_eivec(size, size), : m_eivec(size, size),
m_eivalues(size), m_eivalues(size),
m_isInitialized(false) {} m_isInitialized(false),
m_realSchur(size),
m_matT(size, size),
m_tmp(size)
{}
/** \brief Constructor; computes eigendecomposition of given matrix. /** \brief Constructor; computes eigendecomposition of given matrix.
* *
@ -148,7 +152,10 @@ template<typename _MatrixType> class EigenSolver
EigenSolver(const MatrixType& matrix) EigenSolver(const MatrixType& matrix)
: m_eivec(matrix.rows(), matrix.cols()), : m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()), m_eivalues(matrix.cols()),
m_isInitialized(false) m_isInitialized(false),
m_realSchur(matrix.cols()),
m_matT(matrix.rows(), matrix.cols()),
m_tmp(matrix.cols())
{ {
compute(matrix); compute(matrix);
} }
@ -261,12 +268,17 @@ template<typename _MatrixType> class EigenSolver
EigenSolver& compute(const MatrixType& matrix); EigenSolver& compute(const MatrixType& matrix);
private: private:
void computeEigenvectors(MatrixType& matH); void computeEigenvectors();
protected: protected:
MatrixType m_eivec; MatrixType m_eivec;
EigenvalueType m_eivalues; EigenvalueType m_eivalues;
bool m_isInitialized; bool m_isInitialized;
RealSchur<MatrixType> m_realSchur;
MatrixType m_matT;
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
ColumnVectorType m_tmp;
}; };
template<typename MatrixType> template<typename MatrixType>
@ -324,32 +336,32 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
assert(matrix.cols() == matrix.rows()); assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form. // Reduce to real Schur form.
RealSchur<MatrixType> rs(matrix); m_realSchur.compute(matrix);
MatrixType matT = rs.matrixT(); m_matT = m_realSchur.matrixT();
m_eivec = rs.matrixU(); m_eivec = m_realSchur.matrixU();
// Compute eigenvalues from matT // Compute eigenvalues from matT
m_eivalues.resize(matrix.cols()); m_eivalues.resize(matrix.cols());
int i = 0; int i = 0;
while (i < matrix.cols()) while (i < matrix.cols())
{ {
if (i == matrix.cols() - 1 || matT.coeff(i+1, i) == Scalar(0)) if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
{ {
m_eivalues.coeffRef(i) = matT.coeff(i, i); m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
++i; ++i;
} }
else else
{ {
Scalar p = Scalar(0.5) * (matT.coeff(i, i) - matT.coeff(i+1, i+1)); Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
Scalar z = ei_sqrt(ei_abs(p * p + matT.coeff(i+1, i) * matT.coeff(i, i+1))); Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
m_eivalues.coeffRef(i) = ComplexScalar(matT.coeff(i+1, i+1) + p, z); m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
m_eivalues.coeffRef(i+1) = ComplexScalar(matT.coeff(i+1, i+1) + p, -z); m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
i += 2; i += 2;
} }
} }
// Compute eigenvectors. // Compute eigenvectors.
computeEigenvectors(matT); computeEigenvectors();
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;
@ -376,7 +388,7 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
template<typename MatrixType> template<typename MatrixType>
void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH) void EigenSolver<MatrixType>::computeEigenvectors()
{ {
const int size = m_eivec.cols(); const int size = m_eivec.cols();
const Scalar eps = NumTraits<Scalar>::epsilon(); const Scalar eps = NumTraits<Scalar>::epsilon();
@ -385,7 +397,7 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
Scalar norm = 0.0; Scalar norm = 0.0;
for (int j = 0; j < size; ++j) for (int j = 0; j < size; ++j)
{ {
norm += matH.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum(); norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum();
} }
// Backsubstitute to find vectors of upper triangular form // Backsubstitute to find vectors of upper triangular form
@ -405,11 +417,11 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
Scalar lastr=0, lastw=0; Scalar lastr=0, lastw=0;
int l = n; int l = n;
matH.coeffRef(n,n) = 1.0; m_matT.coeffRef(n,n) = 1.0;
for (int i = n-1; i >= 0; i--) for (int i = n-1; i >= 0; i--)
{ {
Scalar w = matH.coeff(i,i) - p; Scalar w = m_matT.coeff(i,i) - p;
Scalar r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1)); Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
if (m_eivalues.coeff(i).imag() < 0.0) if (m_eivalues.coeff(i).imag() < 0.0)
{ {
@ -422,27 +434,27 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
if (m_eivalues.coeff(i).imag() == 0.0) if (m_eivalues.coeff(i).imag() == 0.0)
{ {
if (w != 0.0) if (w != 0.0)
matH.coeffRef(i,n) = -r / w; m_matT.coeffRef(i,n) = -r / w;
else else
matH.coeffRef(i,n) = -r / (eps * norm); m_matT.coeffRef(i,n) = -r / (eps * norm);
} }
else // Solve real equations else // Solve real equations
{ {
Scalar x = matH.coeff(i,i+1); Scalar x = m_matT.coeff(i,i+1);
Scalar y = matH.coeff(i+1,i); Scalar y = m_matT.coeff(i+1,i);
Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
Scalar t = (x * lastr - lastw * r) / denom; Scalar t = (x * lastr - lastw * r) / denom;
matH.coeffRef(i,n) = t; m_matT.coeffRef(i,n) = t;
if (ei_abs(x) > ei_abs(lastw)) if (ei_abs(x) > ei_abs(lastw))
matH.coeffRef(i+1,n) = (-r - w * t) / x; m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
else else
matH.coeffRef(i+1,n) = (-lastr - y * t) / lastw; m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
} }
// Overflow control // Overflow control
Scalar t = ei_abs(matH.coeff(i,n)); Scalar t = ei_abs(m_matT.coeff(i,n));
if ((eps * t) * t > 1) if ((eps * t) * t > 1)
matH.col(n).tail(size-i) /= t; m_matT.col(n).tail(size-i) /= t;
} }
} }
} }
@ -452,24 +464,24 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
int l = n-1; int l = n-1;
// Last vector component imaginary so matrix is triangular // Last vector component imaginary so matrix is triangular
if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n))) if (ei_abs(m_matT.coeff(n,n-1)) > ei_abs(m_matT.coeff(n-1,n)))
{ {
matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1); m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1); m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
} }
else else
{ {
std::complex<Scalar> cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q); std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
matH.coeffRef(n-1,n-1) = ei_real(cc); m_matT.coeffRef(n-1,n-1) = ei_real(cc);
matH.coeffRef(n-1,n) = ei_imag(cc); m_matT.coeffRef(n-1,n) = ei_imag(cc);
} }
matH.coeffRef(n,n-1) = 0.0; m_matT.coeffRef(n,n-1) = 0.0;
matH.coeffRef(n,n) = 1.0; m_matT.coeffRef(n,n) = 1.0;
for (int i = n-2; i >= 0; i--) for (int i = n-2; i >= 0; i--)
{ {
Scalar ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1)); Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
Scalar sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1)); Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
Scalar w = matH.coeff(i,i) - p; Scalar w = m_matT.coeff(i,i) - p;
if (m_eivalues.coeff(i).imag() < 0.0) if (m_eivalues.coeff(i).imag() < 0.0)
{ {
@ -483,39 +495,39 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
if (m_eivalues.coeff(i).imag() == 0) if (m_eivalues.coeff(i).imag() == 0)
{ {
std::complex<Scalar> cc = cdiv(-ra,-sa,w,q); std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
matH.coeffRef(i,n-1) = ei_real(cc); m_matT.coeffRef(i,n-1) = ei_real(cc);
matH.coeffRef(i,n) = ei_imag(cc); m_matT.coeffRef(i,n) = ei_imag(cc);
} }
else else
{ {
// Solve complex equations // Solve complex equations
Scalar x = matH.coeff(i,i+1); Scalar x = m_matT.coeff(i,i+1);
Scalar y = matH.coeff(i+1,i); Scalar y = m_matT.coeff(i+1,i);
Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q; Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
if ((vr == 0.0) && (vi == 0.0)) if ((vr == 0.0) && (vi == 0.0))
vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(lastw)); vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(lastw));
std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
matH.coeffRef(i,n-1) = ei_real(cc); m_matT.coeffRef(i,n-1) = ei_real(cc);
matH.coeffRef(i,n) = ei_imag(cc); m_matT.coeffRef(i,n) = ei_imag(cc);
if (ei_abs(x) > (ei_abs(lastw) + ei_abs(q))) if (ei_abs(x) > (ei_abs(lastw) + ei_abs(q)))
{ {
matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x; m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x; m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
} }
else else
{ {
cc = cdiv(-lastra-y*matH.coeff(i,n-1),-lastsa-y*matH.coeff(i,n),lastw,q); cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
matH.coeffRef(i+1,n-1) = ei_real(cc); m_matT.coeffRef(i+1,n-1) = ei_real(cc);
matH.coeffRef(i+1,n) = ei_imag(cc); m_matT.coeffRef(i+1,n) = ei_imag(cc);
} }
} }
// Overflow control // Overflow control
Scalar t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n))); Scalar t = std::max(ei_abs(m_matT.coeff(i,n-1)),ei_abs(m_matT.coeff(i,n)));
if ((eps * t) * t > 1) if ((eps * t) * t > 1)
matH.block(i, n-1, size-i, 2) /= t; m_matT.block(i, n-1, size-i, 2) /= t;
} }
} }
@ -525,7 +537,8 @@ void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH)
// Back transformation to get eigenvectors of original matrix // Back transformation to get eigenvectors of original matrix
for (int j = size-1; j >= 0; j--) for (int j = size-1; j >= 0; j--)
{ {
m_eivec.col(j).segment(0, size) = m_eivec.leftCols(j+1) * matH.col(j).segment(0, j+1); m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
m_eivec.col(j) = m_tmp;
} }
} }

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

@ -2,6 +2,7 @@
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
// //
// Eigen is free software; you can redistribute it and/or // Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public // modify it under the terms of the GNU Lesser General Public
@ -110,9 +111,10 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
* Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
*/ */
SelfAdjointEigenSolver() SelfAdjointEigenSolver()
: m_eivec(int(Size), int(Size)), : m_eivec(),
m_eivalues(int(Size)), m_eivalues(),
m_subdiag(int(TridiagonalizationType::SizeMinusOne)) m_tridiag(),
m_subdiag()
{ {
ei_assert(Size!=Dynamic); ei_assert(Size!=Dynamic);
} }
@ -133,7 +135,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
SelfAdjointEigenSolver(int size) SelfAdjointEigenSolver(int size)
: m_eivec(size, size), : m_eivec(size, size),
m_eivalues(size), m_eivalues(size),
m_subdiag(TridiagonalizationType::SizeMinusOne) m_tridiag(size),
m_subdiag(size > 1 ? size - 1 : 1)
{} {}
/** \brief Constructor; computes eigendecomposition of given matrix. /** \brief Constructor; computes eigendecomposition of given matrix.
@ -157,9 +160,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
: m_eivec(matrix.rows(), matrix.cols()), : m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()), m_eivalues(matrix.cols()),
m_subdiag() m_tridiag(matrix.rows()),
m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1)
{ {
if (matrix.rows() > 1) m_subdiag.resize(matrix.rows() - 1);
compute(matrix, computeEigenvectors); compute(matrix, computeEigenvectors);
} }
@ -187,9 +190,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
: m_eivec(matA.rows(), matA.cols()), : m_eivec(matA.rows(), matA.cols()),
m_eivalues(matA.cols()), m_eivalues(matA.cols()),
m_subdiag() m_tridiag(matA.rows()),
m_subdiag(matA.rows() > 1 ? matA.rows() - 1 : 1)
{ {
if (matA.rows() > 1) m_subdiag.resize(matA.rows() - 1);
compute(matA, matB, computeEigenvectors); compute(matA, matB, computeEigenvectors);
} }
@ -351,6 +354,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
protected: protected:
MatrixType m_eivec; MatrixType m_eivec;
RealVectorType m_eivalues; RealVectorType m_eivalues;
TridiagonalizationType m_tridiag;
typename TridiagonalizationType::SubDiagonalType m_subdiag; typename TridiagonalizationType::SubDiagonalType m_subdiag;
#ifndef NDEBUG #ifndef NDEBUG
bool m_eigenvectorsOk; bool m_eigenvectorsOk;
@ -396,14 +400,12 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
return *this; return *this;
} }
m_eivec = matrix; m_tridiag.compute(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)
RealVectorType& diag = m_eivalues; RealVectorType& diag = m_eivalues;
m_subdiag.resize(n-1); diag = m_tridiag.diagonal();
TridiagonalizationType::decomposeInPlace(m_eivec, diag, m_subdiag, computeEigenvectors); m_subdiag = m_tridiag.subDiagonal();
if (computeEigenvectors)
m_eivec = m_tridiag.matrixQ();
int end = n-1; int end = n-1;
int start = 0; int start = 0;

8
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -70,10 +70,10 @@ template<typename _MatrixType> class Tridiagonalization
enum { enum {
Size = MatrixType::RowsAtCompileTime, Size = MatrixType::RowsAtCompileTime,
SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
Options = MatrixType::Options, Options = MatrixType::Options,
MaxSize = MatrixType::MaxRowsAtCompileTime, MaxSize = MatrixType::MaxRowsAtCompileTime,
MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
}; };
typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
@ -108,7 +108,7 @@ template<typename _MatrixType> class Tridiagonalization
* \sa compute() for an example. * \sa compute() for an example.
*/ */
Tridiagonalization(int size = Size==Dynamic ? 2 : Size) Tridiagonalization(int size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size), m_hCoeffs(size-1) : m_matrix(size,size), m_hCoeffs(size > 1 ? size-1 : 1)
{} {}
/** \brief Constructor; computes tridiagonal decomposition of given matrix. /** \brief Constructor; computes tridiagonal decomposition of given matrix.
@ -122,7 +122,7 @@ template<typename _MatrixType> class Tridiagonalization
* Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
*/ */
Tridiagonalization(const MatrixType& matrix) Tridiagonalization(const MatrixType& matrix)
: m_matrix(matrix), m_hCoeffs(matrix.cols()-1) : m_matrix(matrix), m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1)
{ {
_compute(m_matrix, m_hCoeffs); _compute(m_matrix, m_hCoeffs);
} }