Replace local variables by member variables in compute() methods. 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().
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 8e5f131..f6b90d7 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -3,6 +3,7 @@ // // Copyright (C) 2009 Claire Maurice // 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 // modify it under the terms of the GNU Lesser General Public @@ -99,7 +100,8 @@ : m_eivec(), m_eivalues(), m_schur(), - m_isInitialized(false) + m_isInitialized(false), + m_matX() {} /** \brief Default Constructor with memory preallocation @@ -112,7 +114,8 @@ : m_eivec(size, size), m_eivalues(size), m_schur(size), - m_isInitialized(false) + m_isInitialized(false), + m_matX(size, size) {} /** \brief Constructor; computes eigendecomposition of given matrix. @@ -125,7 +128,8 @@ : m_eivec(matrix.rows(),matrix.cols()), m_eivalues(matrix.cols()), m_schur(matrix.rows()), - m_isInitialized(false) + m_isInitialized(false), + m_matX(matrix.rows(),matrix.cols()) { compute(matrix); } @@ -199,6 +203,7 @@ EigenvalueType m_eivalues; ComplexSchur<MatrixType> m_schur; bool m_isInitialized; + EigenvectorType m_matX; }; @@ -217,16 +222,16 @@ // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T. // 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--) { - 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 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) - 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); if(z==ComplexScalar(0)) { @@ -234,12 +239,12 @@ // Use a small value instead, to prevent division by zero. 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) - m_eivec = m_schur.matrixU() * X; + m_eivec.noalias() = m_schur.matrixU() * m_matX; // .. and normalize the eigenvectors for(int k=0 ; k<n ; k++) {
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index f8b953d..7713e04 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -120,7 +120,7 @@ * * \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 * @@ -131,7 +131,11 @@ EigenSolver(int size) : m_eivec(size, 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. * @@ -148,7 +152,10 @@ EigenSolver(const MatrixType& matrix) : m_eivec(matrix.rows(), 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); } @@ -261,12 +268,17 @@ EigenSolver& compute(const MatrixType& matrix); private: - void computeEigenvectors(MatrixType& matH); + void computeEigenvectors(); protected: MatrixType m_eivec; EigenvalueType m_eivalues; 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> @@ -324,32 +336,32 @@ assert(matrix.cols() == matrix.rows()); // Reduce to real Schur form. - RealSchur<MatrixType> rs(matrix); - MatrixType matT = rs.matrixT(); - m_eivec = rs.matrixU(); + m_realSchur.compute(matrix); + m_matT = m_realSchur.matrixT(); + m_eivec = m_realSchur.matrixU(); // Compute eigenvalues from matT m_eivalues.resize(matrix.cols()); int i = 0; 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; } else { - Scalar p = Scalar(0.5) * (matT.coeff(i, i) - matT.coeff(i+1, i+1)); - Scalar z = ei_sqrt(ei_abs(p * p + matT.coeff(i+1, i) * matT.coeff(i, i+1))); - m_eivalues.coeffRef(i) = ComplexScalar(matT.coeff(i+1, i+1) + p, z); - m_eivalues.coeffRef(i+1) = ComplexScalar(matT.coeff(i+1, i+1) + p, -z); + 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 + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); + m_eivalues.coeffRef(i) = ComplexScalar(m_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; } } // Compute eigenvectors. - computeEigenvectors(matT); + computeEigenvectors(); m_isInitialized = true; return *this; @@ -376,7 +388,7 @@ template<typename MatrixType> -void EigenSolver<MatrixType>::computeEigenvectors(MatrixType& matH) +void EigenSolver<MatrixType>::computeEigenvectors() { const int size = m_eivec.cols(); const Scalar eps = NumTraits<Scalar>::epsilon(); @@ -385,7 +397,7 @@ Scalar norm = 0.0; 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 @@ -405,11 +417,11 @@ Scalar lastr=0, lastw=0; int l = n; - matH.coeffRef(n,n) = 1.0; + m_matT.coeffRef(n,n) = 1.0; for (int i = n-1; i >= 0; i--) { - Scalar w = matH.coeff(i,i) - p; - Scalar r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1)); + Scalar w = m_matT.coeff(i,i) - p; + 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) { @@ -422,27 +434,27 @@ if (m_eivalues.coeff(i).imag() == 0.0) { if (w != 0.0) - matH.coeffRef(i,n) = -r / w; + m_matT.coeffRef(i,n) = -r / w; else - matH.coeffRef(i,n) = -r / (eps * norm); + m_matT.coeffRef(i,n) = -r / (eps * norm); } else // Solve real equations { - Scalar x = matH.coeff(i,i+1); - Scalar y = matH.coeff(i+1,i); + Scalar x = m_matT.coeff(i,i+1); + 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 t = (x * lastr - lastw * r) / denom; - matH.coeffRef(i,n) = t; + m_matT.coeffRef(i,n) = t; 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 - matH.coeffRef(i+1,n) = (-lastr - y * t) / lastw; + m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw; } // Overflow control - Scalar t = ei_abs(matH.coeff(i,n)); + Scalar t = ei_abs(m_matT.coeff(i,n)); if ((eps * t) * t > 1) - matH.col(n).tail(size-i) /= t; + m_matT.col(n).tail(size-i) /= t; } } } @@ -452,24 +464,24 @@ int l = n-1; // 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); - matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1); + m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1); + m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1); } else { - std::complex<Scalar> cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q); - matH.coeffRef(n-1,n-1) = ei_real(cc); - matH.coeffRef(n-1,n) = ei_imag(cc); + std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q); + m_matT.coeffRef(n-1,n-1) = ei_real(cc); + m_matT.coeffRef(n-1,n) = ei_imag(cc); } - matH.coeffRef(n,n-1) = 0.0; - matH.coeffRef(n,n) = 1.0; + m_matT.coeffRef(n,n-1) = 0.0; + m_matT.coeffRef(n,n) = 1.0; 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 sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1)); - Scalar w = matH.coeff(i,i) - p; + Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).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 = m_matT.coeff(i,i) - p; if (m_eivalues.coeff(i).imag() < 0.0) { @@ -483,39 +495,39 @@ if (m_eivalues.coeff(i).imag() == 0) { std::complex<Scalar> cc = cdiv(-ra,-sa,w,q); - matH.coeffRef(i,n-1) = ei_real(cc); - matH.coeffRef(i,n) = ei_imag(cc); + m_matT.coeffRef(i,n-1) = ei_real(cc); + m_matT.coeffRef(i,n) = ei_imag(cc); } else { // Solve complex equations - Scalar x = matH.coeff(i,i+1); - Scalar y = matH.coeff(i+1,i); + Scalar x = m_matT.coeff(i,i+1); + 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 vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; 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)); 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); - matH.coeffRef(i,n) = ei_imag(cc); + m_matT.coeffRef(i,n-1) = ei_real(cc); + m_matT.coeffRef(i,n) = ei_imag(cc); 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; - matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x; + m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; + m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x; } else { - cc = cdiv(-lastra-y*matH.coeff(i,n-1),-lastsa-y*matH.coeff(i,n),lastw,q); - matH.coeffRef(i+1,n-1) = ei_real(cc); - matH.coeffRef(i+1,n) = ei_imag(cc); + cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q); + m_matT.coeffRef(i+1,n-1) = ei_real(cc); + m_matT.coeffRef(i+1,n) = ei_imag(cc); } } // 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) - matH.block(i, n-1, size-i, 2) /= t; + m_matT.block(i, n-1, size-i, 2) /= t; } } @@ -525,7 +537,8 @@ // Back transformation to get eigenvectors of original matrix 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; } }
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 00a9d36..25b18dd 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -2,6 +2,7 @@ // for linear algebra. // // 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 // modify it under the terms of the GNU Lesser General Public @@ -110,9 +111,10 @@ * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out */ SelfAdjointEigenSolver() - : m_eivec(int(Size), int(Size)), - m_eivalues(int(Size)), - m_subdiag(int(TridiagonalizationType::SizeMinusOne)) + : m_eivec(), + m_eivalues(), + m_tridiag(), + m_subdiag() { ei_assert(Size!=Dynamic); } @@ -133,7 +135,8 @@ SelfAdjointEigenSolver(int size) : m_eivec(size, 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. @@ -157,9 +160,9 @@ SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) : m_eivec(matrix.rows(), 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); } @@ -187,9 +190,9 @@ SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) : m_eivec(matA.rows(), 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); } @@ -351,6 +354,7 @@ protected: MatrixType m_eivec; RealVectorType m_eivalues; + TridiagonalizationType m_tridiag; typename TridiagonalizationType::SubDiagonalType m_subdiag; #ifndef NDEBUG bool m_eigenvectorsOk; @@ -396,14 +400,12 @@ return *this; } - 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) + m_tridiag.compute(matrix); RealVectorType& diag = m_eivalues; - m_subdiag.resize(n-1); - TridiagonalizationType::decomposeInPlace(m_eivec, diag, m_subdiag, computeEigenvectors); + diag = m_tridiag.diagonal(); + m_subdiag = m_tridiag.subDiagonal(); + if (computeEigenvectors) + m_eivec = m_tridiag.matrixQ(); int end = n-1; int start = 0;
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index 4350986..6ea852a 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -70,10 +70,10 @@ enum { Size = MatrixType::RowsAtCompileTime, - SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, + SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1), Options = MatrixType::Options, 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; @@ -108,7 +108,7 @@ * \sa compute() for an example. */ 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. @@ -122,7 +122,7 @@ * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out */ 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); }