diff options
Diffstat (limited to 'extern/Eigen3/Eigen/src/Cholesky/LLT.h')
| -rw-r--r-- | extern/Eigen3/Eigen/src/Cholesky/LLT.h | 172 |
1 files changed, 137 insertions, 35 deletions
diff --git a/extern/Eigen3/Eigen/src/Cholesky/LLT.h b/extern/Eigen3/Eigen/src/Cholesky/LLT.h index 3bb76b5787f..41d14e532f1 100644 --- a/extern/Eigen3/Eigen/src/Cholesky/LLT.h +++ b/extern/Eigen3/Eigen/src/Cholesky/LLT.h @@ -3,39 +3,28 @@ // // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> // -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_LLT_H #define EIGEN_LLT_H +namespace Eigen { + namespace internal{ template<typename MatrixType, int UpLo> struct LLT_Traits; } -/** \ingroup cholesky_Module +/** \ingroup Cholesky_Module * * \class LLT * * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features * * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition + * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. + * The other triangular part won't be read. * * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite * matrix A such that A = LL^* = U^*U, where L is lower triangular. @@ -49,6 +38,9 @@ template<typename MatrixType, int UpLo> struct LLT_Traits; * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations * has a solution. * + * Example: \include LLT_example.cpp + * Output: \verbinclude LLT_example.out + * * \sa MatrixBase::llt(), class LDLT */ /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) @@ -178,6 +170,9 @@ template<typename _MatrixType, int _UpLo> class LLT inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } + template<typename VectorType> + LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); + protected: /** \internal * Used to compute and store L @@ -190,16 +185,85 @@ template<typename _MatrixType, int _UpLo> class LLT namespace internal { -template<int UpLo> struct llt_inplace; +template<typename Scalar, int UpLo> struct llt_inplace; + +template<typename MatrixType, typename VectorType> +static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + typedef typename MatrixType::ColXpr ColXpr; + typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; + typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; + typedef Matrix<Scalar,Dynamic,1> TempVectorType; + typedef typename TempVectorType::SegmentReturnType TempVecSegment; + + int n = mat.cols(); + eigen_assert(mat.rows()==n && vec.size()==n); + + TempVectorType temp; + + if(sigma>0) + { + // This version is based on Givens rotations. + // It is faster than the other one below, but only works for updates, + // i.e., for sigma > 0 + temp = sqrt(sigma) * vec; + + for(int i=0; i<n; ++i) + { + JacobiRotation<Scalar> g; + g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); + + int rs = n-i-1; + if(rs>0) + { + ColXprSegment x(mat.col(i).tail(rs)); + TempVecSegment y(temp.tail(rs)); + apply_rotation_in_the_plane(x, y, g); + } + } + } + else + { + temp = vec; + RealScalar beta = 1; + for(int j=0; j<n; ++j) + { + RealScalar Ljj = real(mat.coeff(j,j)); + RealScalar dj = abs2(Ljj); + Scalar wj = temp.coeff(j); + RealScalar swj2 = sigma*abs2(wj); + RealScalar gamma = dj*beta + swj2; + + RealScalar x = dj + swj2/beta; + if (x<=RealScalar(0)) + return j; + RealScalar nLjj = sqrt(x); + mat.coeffRef(j,j) = nLjj; + beta += swj2/dj; + + // Update the terms of L + Index rs = n-j-1; + if(rs) + { + temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs); + if(gamma != 0) + mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs); + } + } + } + return -1; +} -template<> struct llt_inplace<Lower> +template<typename Scalar> struct llt_inplace<Scalar, Lower> { + typedef typename NumTraits<Scalar>::Real RealScalar; template<typename MatrixType> static typename MatrixType::Index unblocked(MatrixType& mat) { typedef typename MatrixType::Index Index; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); @@ -254,21 +318,35 @@ template<> struct llt_inplace<Lower> } return -1; } -}; -template<> struct llt_inplace<Upper> + template<typename MatrixType, typename VectorType> + static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) + { + return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); + } +}; + +template<typename Scalar> struct llt_inplace<Scalar, Upper> { + typedef typename NumTraits<Scalar>::Real RealScalar; + template<typename MatrixType> static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); - return llt_inplace<Lower>::unblocked(matt); + return llt_inplace<Scalar, Lower>::unblocked(matt); } template<typename MatrixType> static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); - return llt_inplace<Lower>::blocked(matt); + return llt_inplace<Scalar, Lower>::blocked(matt); + } + template<typename MatrixType, typename VectorType> + static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) + { + Transpose<MatrixType> matt(mat); + return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma); } }; @@ -276,33 +354,35 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> { typedef const TriangularView<const MatrixType, Lower> MatrixL; typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU; - inline static MatrixL getL(const MatrixType& m) { return m; } - inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } + static inline MatrixL getL(const MatrixType& m) { return m; } + static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } static bool inplace_decomposition(MatrixType& m) - { return llt_inplace<Lower>::blocked(m)==-1; } + { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; } }; template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> { typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL; typedef const TriangularView<const MatrixType, Upper> MatrixU; - inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } - inline static MatrixU getU(const MatrixType& m) { return m; } + static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); } + static inline MatrixU getU(const MatrixType& m) { return m; } static bool inplace_decomposition(MatrixType& m) - { return llt_inplace<Upper>::blocked(m)==-1; } + { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; } }; } // end namespace internal /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix * - * * \returns a reference to *this + * + * Example: \include TutorialLinAlgComputeTwice.cpp + * Output: \verbinclude TutorialLinAlgComputeTwice.out */ template<typename MatrixType, int _UpLo> LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a) { - assert(a.rows()==a.cols()); + eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); m_matrix.resize(size, size); m_matrix = a; @@ -314,6 +394,26 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a) return *this; } +/** Performs a rank one update (or dowdate) of the current decomposition. + * If A = LL^* before the rank one update, + * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector + * of same dimension. + */ +template<typename _MatrixType, int _UpLo> +template<typename VectorType> +LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); + eigen_assert(v.size()==m_matrix.cols()); + eigen_assert(m_isInitialized); + if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0) + m_info = NumericalIssue; + else + m_info = Success; + + return *this; +} + namespace internal { template<typename _MatrixType, int UpLo, typename Rhs> struct solve_retval<LLT<_MatrixType, UpLo>, Rhs> @@ -383,4 +483,6 @@ SelfAdjointView<MatrixType, UpLo>::llt() const return LLT<PlainObject,UpLo>(m_matrix); } +} // end namespace Eigen + #endif // EIGEN_LLT_H |