diff options
author | Gael Guennebaud <g.gael@free.fr> | 2011-06-20 15:05:50 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2011-06-20 15:05:50 +0200 |
commit | 2f32e485174c2049a7c364745b2aad091449602b (patch) | |
tree | 7a14c2922c0d6bb88edec8825f1de33ac93e8619 /Eigen/src/Cholesky | |
parent | a55c27a15f8a34a6479f0032ae1fa72b1207a065 (diff) |
New feature: add rank one update in Cholesky decomposition
Diffstat (limited to 'Eigen/src/Cholesky')
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 52 |
1 files changed, 52 insertions, 0 deletions
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index a8fc525e8..33a5c4165 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -178,6 +178,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> + void rankUpdate(const VectorType& vec); + protected: /** \internal * Used to compute and store L @@ -254,6 +257,34 @@ template<> struct llt_inplace<Lower> } return -1; } + + template<typename MatrixType, typename VectorType> + static void rankUpdate(MatrixType& mat, const VectorType& vec) + { + typedef typename MatrixType::ColXpr ColXpr; + typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; + typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; + typedef typename VectorType::SegmentReturnType VecSegment; + + int n = mat.cols(); + eigen_assert(mat.rows()==n && vec.size()==n); + typedef typename MatrixType::Scalar Scalar; + Matrix<Scalar,Dynamic,1> temp(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)); + VecSegment y(temp.tail(rs)); + apply_rotation_in_the_plane(x, y, g); + } + } + } }; template<> struct llt_inplace<Upper> @@ -270,6 +301,12 @@ template<> struct llt_inplace<Upper> Transpose<MatrixType> matt(mat); return llt_inplace<Lower>::blocked(matt); } + template<typename MatrixType, typename VectorType> + static void rankUpdate(MatrixType& mat, const VectorType& vec) + { + Transpose<MatrixType> matt(mat); + return llt_inplace<Lower>::rankUpdate(matt, vec); + } }; template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> @@ -314,6 +351,20 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a) return *this; } +/** Performs a rank one update of the current decomposition. + * If A = LL^* before the rank one update, + * then after it we have LL^* = A + vv^* where \a v must be a vector + * of same dimension. + * + */ +template<typename MatrixType, int _UpLo> +template<typename VectorType> +void LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); + internal::llt_inplace<UpLo>::rankUpdate(m_matrix,v); +} + namespace internal { template<typename _MatrixType, int UpLo, typename Rhs> struct solve_retval<LLT<_MatrixType, UpLo>, Rhs> @@ -384,3 +435,4 @@ SelfAdjointView<MatrixType, UpLo>::llt() const } #endif // EIGEN_LLT_H + |