From 5a36f4a8d1c84e0837ec53bbf9c1881b938ef157 Mon Sep 17 00:00:00 2001 From: Adolfo Rodriguez Tsouroukdissian Date: Mon, 8 Mar 2010 19:31:27 +0100 Subject: [PATCH] Propagate all five matrix template parameters to members and temporaries of decomposition classes. One particular advantage of this is that decomposing matrices with max sizes known at compile time will not allocate. NOTE: The ComplexEigenSolver class currently _does_ allocate (line 135 of Eigenvalues/ComplexEigenSolver.h), but the reason appears to be in the implementation of matrix-matrix products, and not in the decomposition itself. The nomalloc unit test has been extended to verify that decompositions do not allocate when max sizes are specified. There are currently two workarounds to prevent the test from failing (see comments in test/nomalloc.cpp), both of which are related to matrix products that allocate on the stack. --- Eigen/src/Cholesky/LDLT.h | 15 +++-- Eigen/src/Cholesky/LLT.h | 8 ++- Eigen/src/Eigenvalues/ComplexEigenSolver.h | 11 +++- Eigen/src/Eigenvalues/ComplexSchur.h | 9 ++- Eigen/src/Eigenvalues/EigenSolver.h | 14 +++-- .../src/Eigenvalues/HessenbergDecomposition.h | 20 +++--- .../src/Eigenvalues/SelfAdjointEigenSolver.h | 10 ++- Eigen/src/Eigenvalues/Tridiagonalization.h | 18 +++--- Eigen/src/LU/FullPivLU.h | 15 +++-- Eigen/src/LU/PartialPivLU.h | 15 +++-- Eigen/src/QR/ColPivHouseholderQR.h | 17 +++--- Eigen/src/QR/FullPivHouseholderQR.h | 19 +++--- Eigen/src/QR/HouseholderQR.h | 11 ++-- Eigen/src/SVD/SVD.h | 18 +++--- test/nomalloc.cpp | 61 +++++++++++++++++++ 15 files changed, 192 insertions(+), 69 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 8699fe7e0..1a4a4a8eb 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -56,11 +56,18 @@ template class LDLT { public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; - typedef Matrix IntColVectorType; - typedef Matrix IntRowVectorType; + typedef Matrix VectorType; + typedef Matrix IntColVectorType; + typedef Matrix IntRowVectorType; /** \brief Default Constructor. * @@ -201,7 +208,7 @@ LDLT& LDLT::compute(const MatrixType& a) // By using a temorary, packet-aligned products are guarenteed. In the LLT // case this is unnecessary because the diagonal is included and will always // have optimal alignment. - Matrix _temporary(size); + Matrix _temporary(size); for (int j = 0; j < size; ++j) { diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 2e8df7661..d552e4e8a 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -57,9 +57,15 @@ template class LLT { public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; + typedef Matrix VectorType; enum { PacketSize = ei_packet_traits::size, diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 27c52b4dc..dae6091fe 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -41,11 +41,18 @@ template class ComplexEigenSolver { public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef std::complex Complex; - typedef Matrix EigenvalueType; - typedef Matrix EigenvectorType; + typedef Matrix EigenvalueType; + typedef Matrix EigenvectorType; /** * \brief Default Constructor. diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index c45151e82..1404af831 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -45,10 +45,17 @@ template class ComplexSchur { public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef std::complex Complex; - typedef Matrix ComplexMatrixType; + typedef Matrix ComplexMatrixType; enum { Size = MatrixType::RowsAtCompileTime }; diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 3f9e30a6e..579585618 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -45,13 +45,19 @@ template class EigenSolver public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef std::complex Complex; - typedef Matrix EigenvalueType; - typedef Matrix EigenvectorType; - typedef Matrix RealVectorType; - typedef Matrix RealVectorTypeX; + typedef Matrix EigenvalueType; + typedef Matrix EigenvectorType; + typedef Matrix RealVectorType; /** * \brief Default Constructor. diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index e3f64b807..f87c0b842 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -44,17 +44,17 @@ template class HessenbergDecomposition public: typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - enum { Size = MatrixType::RowsAtCompileTime, - SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic - ? Dynamic - : MatrixType::RowsAtCompileTime-1 + SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, + Options = MatrixType::Options, + MaxSize = MatrixType::MaxRowsAtCompileTime, + MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 }; - - typedef Matrix CoeffVectorType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix CoeffVectorType; + typedef Matrix VectorType; /** This constructor initializes a HessenbergDecomposition object for * further use with HessenbergDecomposition::compute() @@ -143,7 +143,7 @@ void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVector { assert(matA.rows()==matA.cols()); int n = matA.rows(); - Matrix temp(n); + VectorType temp(n); for (int i = 0; i::matrixQ() const { int n = m_matrix.rows(); MatrixType matQ = MatrixType::Identity(n,n); - Matrix temp(n); + VectorType temp(n); for (int i = n-2; i>=0; i--) { matQ.corner(BottomRight,n-i-1,n-i-1) diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index e8e9bea97..209624c0a 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -42,13 +42,17 @@ template class SelfAdjointEigenSolver { public: - enum {Size = _MatrixType::RowsAtCompileTime }; typedef _MatrixType MatrixType; + enum { + Size = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef std::complex Complex; - typedef Matrix RealVectorType; - typedef Matrix RealVectorTypeX; + typedef Matrix RealVectorType; typedef Tridiagonalization TridiagonalizationType; // typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType; diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index ffd374eb2..3ef493fa7 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -50,16 +50,18 @@ template class Tridiagonalization enum { Size = MatrixType::RowsAtCompileTime, - SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic - ? Dynamic - : MatrixType::RowsAtCompileTime-1, + SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, + Options = MatrixType::Options, + MaxSize = MatrixType::MaxRowsAtCompileTime, + MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1, PacketSize = ei_packet_traits::size }; - typedef Matrix CoeffVectorType; - typedef Matrix DiagonalType; - typedef Matrix SubDiagonalType; - + typedef Matrix CoeffVectorType; + typedef Matrix DiagonalType; + typedef Matrix SubDiagonalType; + typedef Matrix RowVectorType; + typedef typename ei_meta_if::IsComplex, typename Diagonal::RealReturnType, Diagonal @@ -238,7 +240,7 @@ void Tridiagonalization::matrixQInPlace(MatrixBase* q) con QDerived& matQ = q->derived(); int n = m_matrix.rows(); matQ = MatrixType::Identity(n,n); - Matrix aux(n); + RowVectorType aux(n); for (int i = n-2; i>=0; i--) { matQ.corner(BottomRight,n-i-1,n-i-1) diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index dea6ec41c..4092567a6 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -59,12 +59,19 @@ template class FullPivLU { public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix IntRowVectorType; - typedef Matrix IntColVectorType; - typedef PermutationMatrix PermutationQType; - typedef PermutationMatrix PermutationPType; + typedef Matrix IntRowVectorType; + typedef Matrix IntColVectorType; + typedef PermutationMatrix PermutationQType; + typedef PermutationMatrix PermutationPType; /** * \brief Default Constructor. diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 0dc2d0857..8ff9c5b03 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -62,15 +62,18 @@ template class PartialPivLU public: typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix PermutationVectorType; - typedef PermutationMatrix PermutationType; + typedef Matrix PermutationVectorType; + typedef PermutationMatrix PermutationType; - enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN( - MatrixType::MaxColsAtCompileTime, - MatrixType::MaxRowsAtCompileTime) - }; /** * \brief Default Constructor. diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index be1b22979..7ec477153 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -51,16 +51,19 @@ template class ColPivHouseholderQR RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, - DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime) + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime), + MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN(MaxColsAtCompileTime,MaxRowsAtCompileTime) }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix MatrixQType; - typedef Matrix HCoeffsType; - typedef PermutationMatrix PermutationType; - typedef Matrix IntRowVectorType; - typedef Matrix RowVectorType; - typedef Matrix RealRowVectorType; + typedef Matrix MatrixQType; + typedef Matrix HCoeffsType; + typedef PermutationMatrix PermutationType; + typedef Matrix IntRowVectorType; + typedef Matrix RowVectorType; + typedef Matrix RealRowVectorType; typedef typename HouseholderSequence::ConjugateReturnType HouseholderSequenceType; /** diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 6ba0f45f8..cc0b6e4ed 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -51,17 +51,20 @@ template class FullPivHouseholderQR RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, - DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime) + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime), + MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN(MaxColsAtCompileTime,MaxRowsAtCompileTime) }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix MatrixQType; - typedef Matrix HCoeffsType; - typedef Matrix IntRowVectorType; - typedef PermutationMatrix PermutationType; - typedef Matrix IntColVectorType; - typedef Matrix RowVectorType; - typedef Matrix ColVectorType; + typedef Matrix MatrixQType; + typedef Matrix HCoeffsType; + typedef Matrix IntRowVectorType; + typedef PermutationMatrix PermutationType; + typedef Matrix IntColVectorType; + typedef Matrix RowVectorType; + typedef Matrix ColVectorType; /** \brief Default Constructor. * diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 79d1ea6a0..42ad94030 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -55,13 +55,16 @@ template class HouseholderQR RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, - DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime) + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime), + MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN(MaxColsAtCompileTime,MaxRowsAtCompileTime) }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix::Flags&RowMajorBit ? RowMajor : ColMajor> MatrixQType; - typedef Matrix HCoeffsType; - typedef Matrix RowVectorType; + typedef Matrix::Flags&RowMajorBit ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; + typedef Matrix HCoeffsType; + typedef Matrix RowVectorType; typedef typename HouseholderSequence::ConjugateReturnType HouseholderSequenceType; /** diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index fa3b82ce2..bd1ba3cf3 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -52,15 +52,18 @@ template class SVD ColsAtCompileTime = MatrixType::ColsAtCompileTime, PacketSize = ei_packet_traits::size, AlignmentMask = int(PacketSize)-1, - MinSize = EIGEN_SIZE_MIN(RowsAtCompileTime, ColsAtCompileTime) + MinSize = EIGEN_SIZE_MIN(RowsAtCompileTime, ColsAtCompileTime), + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + MatrixOptions = MatrixType::Options }; - typedef Matrix ColVector; - typedef Matrix RowVector; + typedef Matrix ColVector; + typedef Matrix RowVector; - typedef Matrix MatrixUType; - typedef Matrix MatrixVType; - typedef Matrix SingularValuesType; + typedef Matrix MatrixUType; + typedef Matrix MatrixVType; + typedef Matrix SingularValuesType; /** * \brief Default Constructor. @@ -195,7 +198,8 @@ SVD& SVD::compute(const MatrixType& matrix) bool convergence = true; Scalar eps = NumTraits::dummy_precision(); - Matrix rv1(n); + RowVector rv1(n); + g = scale = anorm = 0; // Householder reduction to bidiagonal form. for (i=0; i +#include +#include +#include +#include template void nomalloc(const MatrixType& m) { @@ -73,6 +78,58 @@ template void nomalloc(const MatrixType& m) } } +void ctms_decompositions() +{ + const int maxSize = 16; + const int size = 12; + + typedef Eigen::Matrix Matrix; + + typedef Eigen::Matrix Vector; + + typedef Eigen::Matrix, + Eigen::Dynamic, Eigen::Dynamic, + Eigen::ColMajor | Eigen::AutoAlign, + maxSize, maxSize> ComplexMatrix; + + const Matrix A(Matrix::Random(size, size)); + const ComplexMatrix complexA(ComplexMatrix::Random(size, size)); +// const Matrix saA = A.adjoint() * A; // NOTE: This product allocates on the stack. The two following lines are a kludgy workaround + Matrix saA(Matrix::Constant(size, size, 1.0)); + saA.diagonal().setConstant(2.0); + + // Cholesky module + Eigen::LLT LLT; LLT.compute(A); + Eigen::LDLT LDLT; LDLT.compute(A); + + // Eigenvalues module + Eigen::HessenbergDecomposition hessDecomp; hessDecomp.compute(complexA); + Eigen::ComplexSchur cSchur(size); cSchur.compute(complexA); + Eigen::ComplexEigenSolver cEigSolver; //cEigSolver.compute(complexA); // NOTE: Commented-out because makes test fail (L135 of ComplexEigenSolver.h has a product that allocates on the stack) + Eigen::EigenSolver eigSolver; eigSolver.compute(A); + Eigen::SelfAdjointEigenSolver saEigSolver(size); saEigSolver.compute(saA); + Eigen::Tridiagonalization tridiag; tridiag.compute(saA); + + // LU module + Eigen::PartialPivLU ppLU; ppLU.compute(A); + Eigen::FullPivLU fpLU; fpLU.compute(A); + + // QR module + Eigen::HouseholderQR hQR; hQR.compute(A); + Eigen::ColPivHouseholderQR cpQR; cpQR.compute(A); + Eigen::FullPivHouseholderQR fpQR; fpQR.compute(A); + + // SVD module + Eigen::JacobiSVD jSVD; jSVD.compute(A); + Eigen::SVD svd; svd.compute(A); +} + void test_nomalloc() { // check that our operator new is indeed called: @@ -80,4 +137,8 @@ void test_nomalloc() CALL_SUBTEST(nomalloc(Matrix()) ); CALL_SUBTEST(nomalloc(Matrix4d()) ); CALL_SUBTEST(nomalloc(Matrix()) ); + + // Check decomposition modules with dynamic matrices that have a known compile-time max size (ctms) + CALL_SUBTEST(ctms_decompositions()); + }