From 79877a9917c6482930ef2ac15e2de7ef8a7a757e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 7 Jul 2009 15:32:21 +0200 Subject: [PATCH] * take advantage of new possibilies in LLT (mat -= product) * fix Block::operator+= product which was not optimized * fix some compilation issues --- Eigen/src/Cholesky/LLT.h | 31 +++++++------------ Eigen/src/Core/MapBase.h | 3 +- Eigen/src/Core/SolveTriangular.h | 3 +- Eigen/src/Core/products/GeneralMatrixMatrix.h | 1 - test/cholesky.cpp | 8 +++++ 5 files changed, 23 insertions(+), 23 deletions(-) diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index fccd53459..53e97cd3f 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -26,7 +26,7 @@ #define EIGEN_LLT_H template struct LLT_Traits; - + /** \ingroup cholesky_Module * * \class LLT @@ -70,7 +70,7 @@ template class LLT public: - /** + /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -94,7 +94,7 @@ template class LLT /** \returns a view of the lower triangular matrix L */ inline typename Traits::MatrixL matrixL() const - { + { ei_assert(m_isInitialized && "LLT is not initialized."); return Traits::getL(m_matrix); } @@ -123,7 +123,7 @@ bool ei_inplace_llt_lo(MatrixType& mat) typedef typename MatrixType::RealScalar RealScalar; assert(mat.rows()==mat.cols()); const int size = mat.rows(); - Matrix aux(size); + // The biggest overall is the point of reference to which further diagonals // are compared; if any diagonal is negligible compared // to the largest overall, the algorithm bails. This cutoff is suggested @@ -147,16 +147,10 @@ bool ei_inplace_llt_lo(MatrixType& mat) mat.coeffRef(j,j) = x = ei_sqrt(x); int endSize = size-j-1; - if (endSize>0) { - // Note that when all matrix columns have good alignment, then the following - // product is guaranteed to be optimal with respect to alignment. - aux.end(endSize) = - (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()).lazy(); - - mat.col(j).end(endSize) = (mat.col(j).end(endSize) - aux.end(endSize) ) / x; - // TODO improve the products so that the following is efficient: -// mat.col(j).end(endSize) -= (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()); -// mat.col(j).end(endSize) *= Scalar(1)/x; + if (endSize>0) + { + mat.col(j).end(endSize) -= (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()).lazy(); + mat.col(j).end(endSize) *= RealScalar(1)/x; } } @@ -170,8 +164,7 @@ bool ei_inplace_llt_up(MatrixType& mat) typedef typename MatrixType::RealScalar RealScalar; assert(mat.rows()==mat.cols()); const int size = mat.rows(); - Matrix aux(size); - + const RealScalar cutoff = machine_epsilon() * size * mat.diagonal().cwise().abs().maxCoeff(); RealScalar x; x = ei_real(mat.coeff(0,0)); @@ -190,10 +183,8 @@ bool ei_inplace_llt_up(MatrixType& mat) int endSize = size-j-1; if (endSize>0) { - aux.start(endSize) = - (mat.block(0, j+1, j, endSize).adjoint() * mat.col(j).start(j)).lazy(); - - mat.row(j).end(endSize) = (mat.row(j).end(endSize) - aux.start(endSize).adjoint() ) / x; + mat.row(j).end(endSize) -= (mat.col(j).start(j).adjoint() * mat.block(0, j+1, j, endSize)).lazy(); + mat.row(j).end(endSize) *= RealScalar(1)/x; } } diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index ec9ab3c5e..eddb19b03 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -172,8 +172,9 @@ template class MapBase } using Base::operator=; - using Base::operator*=; + using Base::operator+=; + using Base::operator-=; template Derived& operator+=(const MatrixBase& other) diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index f717a3875..8a6bf3af3 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -176,7 +176,8 @@ struct ei_solve_triangular_selector IsLowerTriangular ? size-endBlock : endBlock+1, &(lhs.const_cast_derived().coeffRef(IsLowerTriangular ? endBlock : 0, IsLowerTriangular ? startBlock : endBlock+1)), lhs.stride(), - btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c))); + btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c)), + Scalar(1)); // if (IsLowerTriangular) // other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock) // * other.col(c).block(startBlock,endBlock-startBlock)).lazy(); diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index e7b8586be..4630e5040 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -151,7 +151,6 @@ static void ei_cache_friendly_product( const Scalar* b3 = &rhs[(j2+3)*rhsStride + k2]; if (hasAlpha) { - std::cerr << "* by " << alpha << "\n"; for(int k=0; k void cholesky(const MatrixType& m) VERIFY_IS_APPROX(symm * vecX, vecB); chol.solve(matB, &matX); VERIFY_IS_APPROX(symm * matX, matB); + + // test the upper mode + LLT cholup(symm); + VERIFY_IS_APPROX(symm, cholup.matrixL().toDense() * chol.matrixL().adjoint().toDense()); + cholup.solve(vecB, &vecX); + VERIFY_IS_APPROX(symm * vecX, vecB); + cholup.solve(matB, &matX); + VERIFY_IS_APPROX(symm * matX, matB); } int sign = ei_random()%2 ? 1 : -1;