Jaesung Oh / Eigen

Eigne Matrix Class Library

Dependents:   MPC_current_control HydraulicControlBoard_SW AHRS Test_ekf ... more

src/Cholesky/LLT.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 5 //
ykuroda 0:13a5d365ba16 6 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 7 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 9
ykuroda 0:13a5d365ba16 10 #ifndef EIGEN_LLT_H
ykuroda 0:13a5d365ba16 11 #define EIGEN_LLT_H
ykuroda 0:13a5d365ba16 12
ykuroda 0:13a5d365ba16 13 namespace Eigen {
ykuroda 0:13a5d365ba16 14
ykuroda 0:13a5d365ba16 15 namespace internal{
ykuroda 0:13a5d365ba16 16 template<typename MatrixType, int UpLo> struct LLT_Traits;
ykuroda 0:13a5d365ba16 17 }
ykuroda 0:13a5d365ba16 18
ykuroda 0:13a5d365ba16 19 /** \ingroup Cholesky_Module
ykuroda 0:13a5d365ba16 20 *
ykuroda 0:13a5d365ba16 21 * \class LLT
ykuroda 0:13a5d365ba16 22 *
ykuroda 0:13a5d365ba16 23 * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
ykuroda 0:13a5d365ba16 24 *
ykuroda 0:13a5d365ba16 25 * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
ykuroda 0:13a5d365ba16 26 * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
ykuroda 0:13a5d365ba16 27 * The other triangular part won't be read.
ykuroda 0:13a5d365ba16 28 *
ykuroda 0:13a5d365ba16 29 * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
ykuroda 0:13a5d365ba16 30 * matrix A such that A = LL^* = U^*U, where L is lower triangular.
ykuroda 0:13a5d365ba16 31 *
ykuroda 0:13a5d365ba16 32 * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b,
ykuroda 0:13a5d365ba16 33 * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
ykuroda 0:13a5d365ba16 34 * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
ykuroda 0:13a5d365ba16 35 * situations like generalised eigen problems with hermitian matrices.
ykuroda 0:13a5d365ba16 36 *
ykuroda 0:13a5d365ba16 37 * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
ykuroda 0:13a5d365ba16 38 * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
ykuroda 0:13a5d365ba16 39 * has a solution.
ykuroda 0:13a5d365ba16 40 *
ykuroda 0:13a5d365ba16 41 * Example: \include LLT_example.cpp
ykuroda 0:13a5d365ba16 42 * Output: \verbinclude LLT_example.out
ykuroda 0:13a5d365ba16 43 *
ykuroda 0:13a5d365ba16 44 * \sa MatrixBase::llt(), class LDLT
ykuroda 0:13a5d365ba16 45 */
ykuroda 0:13a5d365ba16 46 /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
ykuroda 0:13a5d365ba16 47 * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
ykuroda 0:13a5d365ba16 48 * the strict lower part does not have to store correct values.
ykuroda 0:13a5d365ba16 49 */
ykuroda 0:13a5d365ba16 50 template<typename _MatrixType, int _UpLo> class LLT
ykuroda 0:13a5d365ba16 51 {
ykuroda 0:13a5d365ba16 52 public:
ykuroda 0:13a5d365ba16 53 typedef _MatrixType MatrixType;
ykuroda 0:13a5d365ba16 54 enum {
ykuroda 0:13a5d365ba16 55 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 56 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 57 Options = MatrixType::Options,
ykuroda 0:13a5d365ba16 58 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
ykuroda 0:13a5d365ba16 59 };
ykuroda 0:13a5d365ba16 60 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 61 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 62 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 63
ykuroda 0:13a5d365ba16 64 enum {
ykuroda 0:13a5d365ba16 65 PacketSize = internal::packet_traits<Scalar>::size,
ykuroda 0:13a5d365ba16 66 AlignmentMask = int(PacketSize)-1,
ykuroda 0:13a5d365ba16 67 UpLo = _UpLo
ykuroda 0:13a5d365ba16 68 };
ykuroda 0:13a5d365ba16 69
ykuroda 0:13a5d365ba16 70 typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
ykuroda 0:13a5d365ba16 71
ykuroda 0:13a5d365ba16 72 /**
ykuroda 0:13a5d365ba16 73 * \brief Default Constructor.
ykuroda 0:13a5d365ba16 74 *
ykuroda 0:13a5d365ba16 75 * The default constructor is useful in cases in which the user intends to
ykuroda 0:13a5d365ba16 76 * perform decompositions via LLT::compute(const MatrixType&).
ykuroda 0:13a5d365ba16 77 */
ykuroda 0:13a5d365ba16 78 LLT() : m_matrix(), m_isInitialized(false) {}
ykuroda 0:13a5d365ba16 79
ykuroda 0:13a5d365ba16 80 /** \brief Default Constructor with memory preallocation
ykuroda 0:13a5d365ba16 81 *
ykuroda 0:13a5d365ba16 82 * Like the default constructor but with preallocation of the internal data
ykuroda 0:13a5d365ba16 83 * according to the specified problem \a size.
ykuroda 0:13a5d365ba16 84 * \sa LLT()
ykuroda 0:13a5d365ba16 85 */
ykuroda 0:13a5d365ba16 86 LLT(Index size) : m_matrix(size, size),
ykuroda 0:13a5d365ba16 87 m_isInitialized(false) {}
ykuroda 0:13a5d365ba16 88
ykuroda 0:13a5d365ba16 89 LLT(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 90 : m_matrix(matrix.rows(), matrix.cols()),
ykuroda 0:13a5d365ba16 91 m_isInitialized(false)
ykuroda 0:13a5d365ba16 92 {
ykuroda 0:13a5d365ba16 93 compute(matrix);
ykuroda 0:13a5d365ba16 94 }
ykuroda 0:13a5d365ba16 95
ykuroda 0:13a5d365ba16 96 /** \returns a view of the upper triangular matrix U */
ykuroda 0:13a5d365ba16 97 inline typename Traits::MatrixU matrixU() const
ykuroda 0:13a5d365ba16 98 {
ykuroda 0:13a5d365ba16 99 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 100 return Traits::getU(m_matrix);
ykuroda 0:13a5d365ba16 101 }
ykuroda 0:13a5d365ba16 102
ykuroda 0:13a5d365ba16 103 /** \returns a view of the lower triangular matrix L */
ykuroda 0:13a5d365ba16 104 inline typename Traits::MatrixL matrixL() const
ykuroda 0:13a5d365ba16 105 {
ykuroda 0:13a5d365ba16 106 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 107 return Traits::getL(m_matrix);
ykuroda 0:13a5d365ba16 108 }
ykuroda 0:13a5d365ba16 109
ykuroda 0:13a5d365ba16 110 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
ykuroda 0:13a5d365ba16 111 *
ykuroda 0:13a5d365ba16 112 * Since this LLT class assumes anyway that the matrix A is invertible, the solution
ykuroda 0:13a5d365ba16 113 * theoretically exists and is unique regardless of b.
ykuroda 0:13a5d365ba16 114 *
ykuroda 0:13a5d365ba16 115 * Example: \include LLT_solve.cpp
ykuroda 0:13a5d365ba16 116 * Output: \verbinclude LLT_solve.out
ykuroda 0:13a5d365ba16 117 *
ykuroda 0:13a5d365ba16 118 * \sa solveInPlace(), MatrixBase::llt()
ykuroda 0:13a5d365ba16 119 */
ykuroda 0:13a5d365ba16 120 template<typename Rhs>
ykuroda 0:13a5d365ba16 121 inline const internal::solve_retval<LLT, Rhs>
ykuroda 0:13a5d365ba16 122 solve(const MatrixBase<Rhs>& b) const
ykuroda 0:13a5d365ba16 123 {
ykuroda 0:13a5d365ba16 124 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 125 eigen_assert(m_matrix.rows()==b.rows()
ykuroda 0:13a5d365ba16 126 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
ykuroda 0:13a5d365ba16 127 return internal::solve_retval<LLT, Rhs>(*this, b.derived());
ykuroda 0:13a5d365ba16 128 }
ykuroda 0:13a5d365ba16 129
ykuroda 0:13a5d365ba16 130 #ifdef EIGEN2_SUPPORT
ykuroda 0:13a5d365ba16 131 template<typename OtherDerived, typename ResultType>
ykuroda 0:13a5d365ba16 132 bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
ykuroda 0:13a5d365ba16 133 {
ykuroda 0:13a5d365ba16 134 *result = this->solve(b);
ykuroda 0:13a5d365ba16 135 return true;
ykuroda 0:13a5d365ba16 136 }
ykuroda 0:13a5d365ba16 137
ykuroda 0:13a5d365ba16 138 bool isPositiveDefinite() const { return true; }
ykuroda 0:13a5d365ba16 139 #endif
ykuroda 0:13a5d365ba16 140
ykuroda 0:13a5d365ba16 141 template<typename Derived>
ykuroda 0:13a5d365ba16 142 void solveInPlace(MatrixBase<Derived> &bAndX) const;
ykuroda 0:13a5d365ba16 143
ykuroda 0:13a5d365ba16 144 LLT& compute(const MatrixType& matrix);
ykuroda 0:13a5d365ba16 145
ykuroda 0:13a5d365ba16 146 /** \returns the LLT decomposition matrix
ykuroda 0:13a5d365ba16 147 *
ykuroda 0:13a5d365ba16 148 * TODO: document the storage layout
ykuroda 0:13a5d365ba16 149 */
ykuroda 0:13a5d365ba16 150 inline const MatrixType& matrixLLT() const
ykuroda 0:13a5d365ba16 151 {
ykuroda 0:13a5d365ba16 152 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 153 return m_matrix;
ykuroda 0:13a5d365ba16 154 }
ykuroda 0:13a5d365ba16 155
ykuroda 0:13a5d365ba16 156 MatrixType reconstructedMatrix() const;
ykuroda 0:13a5d365ba16 157
ykuroda 0:13a5d365ba16 158
ykuroda 0:13a5d365ba16 159 /** \brief Reports whether previous computation was successful.
ykuroda 0:13a5d365ba16 160 *
ykuroda 0:13a5d365ba16 161 * \returns \c Success if computation was succesful,
ykuroda 0:13a5d365ba16 162 * \c NumericalIssue if the matrix.appears to be negative.
ykuroda 0:13a5d365ba16 163 */
ykuroda 0:13a5d365ba16 164 ComputationInfo info() const
ykuroda 0:13a5d365ba16 165 {
ykuroda 0:13a5d365ba16 166 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 167 return m_info;
ykuroda 0:13a5d365ba16 168 }
ykuroda 0:13a5d365ba16 169
ykuroda 0:13a5d365ba16 170 inline Index rows() const { return m_matrix.rows(); }
ykuroda 0:13a5d365ba16 171 inline Index cols() const { return m_matrix.cols(); }
ykuroda 0:13a5d365ba16 172
ykuroda 0:13a5d365ba16 173 template<typename VectorType>
ykuroda 0:13a5d365ba16 174 LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
ykuroda 0:13a5d365ba16 175
ykuroda 0:13a5d365ba16 176 protected:
ykuroda 0:13a5d365ba16 177
ykuroda 0:13a5d365ba16 178 static void check_template_parameters()
ykuroda 0:13a5d365ba16 179 {
ykuroda 0:13a5d365ba16 180 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
ykuroda 0:13a5d365ba16 181 }
ykuroda 0:13a5d365ba16 182
ykuroda 0:13a5d365ba16 183 /** \internal
ykuroda 0:13a5d365ba16 184 * Used to compute and store L
ykuroda 0:13a5d365ba16 185 * The strict upper part is not used and even not initialized.
ykuroda 0:13a5d365ba16 186 */
ykuroda 0:13a5d365ba16 187 MatrixType m_matrix;
ykuroda 0:13a5d365ba16 188 bool m_isInitialized;
ykuroda 0:13a5d365ba16 189 ComputationInfo m_info;
ykuroda 0:13a5d365ba16 190 };
ykuroda 0:13a5d365ba16 191
ykuroda 0:13a5d365ba16 192 namespace internal {
ykuroda 0:13a5d365ba16 193
ykuroda 0:13a5d365ba16 194 template<typename Scalar, int UpLo> struct llt_inplace;
ykuroda 0:13a5d365ba16 195
ykuroda 0:13a5d365ba16 196 template<typename MatrixType, typename VectorType>
ykuroda 0:13a5d365ba16 197 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
ykuroda 0:13a5d365ba16 198 {
ykuroda 0:13a5d365ba16 199 using std::sqrt;
ykuroda 0:13a5d365ba16 200 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 201 typedef typename MatrixType::RealScalar RealScalar;
ykuroda 0:13a5d365ba16 202 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 203 typedef typename MatrixType::ColXpr ColXpr;
ykuroda 0:13a5d365ba16 204 typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
ykuroda 0:13a5d365ba16 205 typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
ykuroda 0:13a5d365ba16 206 typedef Matrix<Scalar,Dynamic,1> TempVectorType;
ykuroda 0:13a5d365ba16 207 typedef typename TempVectorType::SegmentReturnType TempVecSegment;
ykuroda 0:13a5d365ba16 208
ykuroda 0:13a5d365ba16 209 Index n = mat.cols();
ykuroda 0:13a5d365ba16 210 eigen_assert(mat.rows()==n && vec.size()==n);
ykuroda 0:13a5d365ba16 211
ykuroda 0:13a5d365ba16 212 TempVectorType temp;
ykuroda 0:13a5d365ba16 213
ykuroda 0:13a5d365ba16 214 if(sigma>0)
ykuroda 0:13a5d365ba16 215 {
ykuroda 0:13a5d365ba16 216 // This version is based on Givens rotations.
ykuroda 0:13a5d365ba16 217 // It is faster than the other one below, but only works for updates,
ykuroda 0:13a5d365ba16 218 // i.e., for sigma > 0
ykuroda 0:13a5d365ba16 219 temp = sqrt(sigma) * vec;
ykuroda 0:13a5d365ba16 220
ykuroda 0:13a5d365ba16 221 for(Index i=0; i<n; ++i)
ykuroda 0:13a5d365ba16 222 {
ykuroda 0:13a5d365ba16 223 JacobiRotation<Scalar> g;
ykuroda 0:13a5d365ba16 224 g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
ykuroda 0:13a5d365ba16 225
ykuroda 0:13a5d365ba16 226 Index rs = n-i-1;
ykuroda 0:13a5d365ba16 227 if(rs>0)
ykuroda 0:13a5d365ba16 228 {
ykuroda 0:13a5d365ba16 229 ColXprSegment x(mat.col(i).tail(rs));
ykuroda 0:13a5d365ba16 230 TempVecSegment y(temp.tail(rs));
ykuroda 0:13a5d365ba16 231 apply_rotation_in_the_plane(x, y, g);
ykuroda 0:13a5d365ba16 232 }
ykuroda 0:13a5d365ba16 233 }
ykuroda 0:13a5d365ba16 234 }
ykuroda 0:13a5d365ba16 235 else
ykuroda 0:13a5d365ba16 236 {
ykuroda 0:13a5d365ba16 237 temp = vec;
ykuroda 0:13a5d365ba16 238 RealScalar beta = 1;
ykuroda 0:13a5d365ba16 239 for(Index j=0; j<n; ++j)
ykuroda 0:13a5d365ba16 240 {
ykuroda 0:13a5d365ba16 241 RealScalar Ljj = numext::real(mat.coeff(j,j));
ykuroda 0:13a5d365ba16 242 RealScalar dj = numext::abs2(Ljj);
ykuroda 0:13a5d365ba16 243 Scalar wj = temp.coeff(j);
ykuroda 0:13a5d365ba16 244 RealScalar swj2 = sigma*numext::abs2(wj);
ykuroda 0:13a5d365ba16 245 RealScalar gamma = dj*beta + swj2;
ykuroda 0:13a5d365ba16 246
ykuroda 0:13a5d365ba16 247 RealScalar x = dj + swj2/beta;
ykuroda 0:13a5d365ba16 248 if (x<=RealScalar(0))
ykuroda 0:13a5d365ba16 249 return j;
ykuroda 0:13a5d365ba16 250 RealScalar nLjj = sqrt(x);
ykuroda 0:13a5d365ba16 251 mat.coeffRef(j,j) = nLjj;
ykuroda 0:13a5d365ba16 252 beta += swj2/dj;
ykuroda 0:13a5d365ba16 253
ykuroda 0:13a5d365ba16 254 // Update the terms of L
ykuroda 0:13a5d365ba16 255 Index rs = n-j-1;
ykuroda 0:13a5d365ba16 256 if(rs)
ykuroda 0:13a5d365ba16 257 {
ykuroda 0:13a5d365ba16 258 temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
ykuroda 0:13a5d365ba16 259 if(gamma != 0)
ykuroda 0:13a5d365ba16 260 mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
ykuroda 0:13a5d365ba16 261 }
ykuroda 0:13a5d365ba16 262 }
ykuroda 0:13a5d365ba16 263 }
ykuroda 0:13a5d365ba16 264 return -1;
ykuroda 0:13a5d365ba16 265 }
ykuroda 0:13a5d365ba16 266
ykuroda 0:13a5d365ba16 267 template<typename Scalar> struct llt_inplace<Scalar, Lower>
ykuroda 0:13a5d365ba16 268 {
ykuroda 0:13a5d365ba16 269 typedef typename NumTraits<Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 270 template<typename MatrixType>
ykuroda 0:13a5d365ba16 271 static typename MatrixType::Index unblocked(MatrixType& mat)
ykuroda 0:13a5d365ba16 272 {
ykuroda 0:13a5d365ba16 273 using std::sqrt;
ykuroda 0:13a5d365ba16 274 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 275
ykuroda 0:13a5d365ba16 276 eigen_assert(mat.rows()==mat.cols());
ykuroda 0:13a5d365ba16 277 const Index size = mat.rows();
ykuroda 0:13a5d365ba16 278 for(Index k = 0; k < size; ++k)
ykuroda 0:13a5d365ba16 279 {
ykuroda 0:13a5d365ba16 280 Index rs = size-k-1; // remaining size
ykuroda 0:13a5d365ba16 281
ykuroda 0:13a5d365ba16 282 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
ykuroda 0:13a5d365ba16 283 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
ykuroda 0:13a5d365ba16 284 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
ykuroda 0:13a5d365ba16 285
ykuroda 0:13a5d365ba16 286 RealScalar x = numext::real(mat.coeff(k,k));
ykuroda 0:13a5d365ba16 287 if (k>0) x -= A10.squaredNorm();
ykuroda 0:13a5d365ba16 288 if (x<=RealScalar(0))
ykuroda 0:13a5d365ba16 289 return k;
ykuroda 0:13a5d365ba16 290 mat.coeffRef(k,k) = x = sqrt(x);
ykuroda 0:13a5d365ba16 291 if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
ykuroda 0:13a5d365ba16 292 if (rs>0) A21 /= x;
ykuroda 0:13a5d365ba16 293 }
ykuroda 0:13a5d365ba16 294 return -1;
ykuroda 0:13a5d365ba16 295 }
ykuroda 0:13a5d365ba16 296
ykuroda 0:13a5d365ba16 297 template<typename MatrixType>
ykuroda 0:13a5d365ba16 298 static typename MatrixType::Index blocked(MatrixType& m)
ykuroda 0:13a5d365ba16 299 {
ykuroda 0:13a5d365ba16 300 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 301 eigen_assert(m.rows()==m.cols());
ykuroda 0:13a5d365ba16 302 Index size = m.rows();
ykuroda 0:13a5d365ba16 303 if(size<32)
ykuroda 0:13a5d365ba16 304 return unblocked(m);
ykuroda 0:13a5d365ba16 305
ykuroda 0:13a5d365ba16 306 Index blockSize = size/8;
ykuroda 0:13a5d365ba16 307 blockSize = (blockSize/16)*16;
ykuroda 0:13a5d365ba16 308 blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
ykuroda 0:13a5d365ba16 309
ykuroda 0:13a5d365ba16 310 for (Index k=0; k<size; k+=blockSize)
ykuroda 0:13a5d365ba16 311 {
ykuroda 0:13a5d365ba16 312 // partition the matrix:
ykuroda 0:13a5d365ba16 313 // A00 | - | -
ykuroda 0:13a5d365ba16 314 // lu = A10 | A11 | -
ykuroda 0:13a5d365ba16 315 // A20 | A21 | A22
ykuroda 0:13a5d365ba16 316 Index bs = (std::min)(blockSize, size-k);
ykuroda 0:13a5d365ba16 317 Index rs = size - k - bs;
ykuroda 0:13a5d365ba16 318 Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
ykuroda 0:13a5d365ba16 319 Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
ykuroda 0:13a5d365ba16 320 Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
ykuroda 0:13a5d365ba16 321
ykuroda 0:13a5d365ba16 322 Index ret;
ykuroda 0:13a5d365ba16 323 if((ret=unblocked(A11))>=0) return k+ret;
ykuroda 0:13a5d365ba16 324 if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
ykuroda 0:13a5d365ba16 325 if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
ykuroda 0:13a5d365ba16 326 }
ykuroda 0:13a5d365ba16 327 return -1;
ykuroda 0:13a5d365ba16 328 }
ykuroda 0:13a5d365ba16 329
ykuroda 0:13a5d365ba16 330 template<typename MatrixType, typename VectorType>
ykuroda 0:13a5d365ba16 331 static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
ykuroda 0:13a5d365ba16 332 {
ykuroda 0:13a5d365ba16 333 return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
ykuroda 0:13a5d365ba16 334 }
ykuroda 0:13a5d365ba16 335 };
ykuroda 0:13a5d365ba16 336
ykuroda 0:13a5d365ba16 337 template<typename Scalar> struct llt_inplace<Scalar, Upper>
ykuroda 0:13a5d365ba16 338 {
ykuroda 0:13a5d365ba16 339 typedef typename NumTraits<Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 340
ykuroda 0:13a5d365ba16 341 template<typename MatrixType>
ykuroda 0:13a5d365ba16 342 static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
ykuroda 0:13a5d365ba16 343 {
ykuroda 0:13a5d365ba16 344 Transpose<MatrixType> matt(mat);
ykuroda 0:13a5d365ba16 345 return llt_inplace<Scalar, Lower>::unblocked(matt);
ykuroda 0:13a5d365ba16 346 }
ykuroda 0:13a5d365ba16 347 template<typename MatrixType>
ykuroda 0:13a5d365ba16 348 static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
ykuroda 0:13a5d365ba16 349 {
ykuroda 0:13a5d365ba16 350 Transpose<MatrixType> matt(mat);
ykuroda 0:13a5d365ba16 351 return llt_inplace<Scalar, Lower>::blocked(matt);
ykuroda 0:13a5d365ba16 352 }
ykuroda 0:13a5d365ba16 353 template<typename MatrixType, typename VectorType>
ykuroda 0:13a5d365ba16 354 static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
ykuroda 0:13a5d365ba16 355 {
ykuroda 0:13a5d365ba16 356 Transpose<MatrixType> matt(mat);
ykuroda 0:13a5d365ba16 357 return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
ykuroda 0:13a5d365ba16 358 }
ykuroda 0:13a5d365ba16 359 };
ykuroda 0:13a5d365ba16 360
ykuroda 0:13a5d365ba16 361 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
ykuroda 0:13a5d365ba16 362 {
ykuroda 0:13a5d365ba16 363 typedef const TriangularView<const MatrixType, Lower> MatrixL;
ykuroda 0:13a5d365ba16 364 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
ykuroda 0:13a5d365ba16 365 static inline MatrixL getL(const MatrixType& m) { return m; }
ykuroda 0:13a5d365ba16 366 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
ykuroda 0:13a5d365ba16 367 static bool inplace_decomposition(MatrixType& m)
ykuroda 0:13a5d365ba16 368 { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
ykuroda 0:13a5d365ba16 369 };
ykuroda 0:13a5d365ba16 370
ykuroda 0:13a5d365ba16 371 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
ykuroda 0:13a5d365ba16 372 {
ykuroda 0:13a5d365ba16 373 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
ykuroda 0:13a5d365ba16 374 typedef const TriangularView<const MatrixType, Upper> MatrixU;
ykuroda 0:13a5d365ba16 375 static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
ykuroda 0:13a5d365ba16 376 static inline MatrixU getU(const MatrixType& m) { return m; }
ykuroda 0:13a5d365ba16 377 static bool inplace_decomposition(MatrixType& m)
ykuroda 0:13a5d365ba16 378 { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
ykuroda 0:13a5d365ba16 379 };
ykuroda 0:13a5d365ba16 380
ykuroda 0:13a5d365ba16 381 } // end namespace internal
ykuroda 0:13a5d365ba16 382
ykuroda 0:13a5d365ba16 383 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
ykuroda 0:13a5d365ba16 384 *
ykuroda 0:13a5d365ba16 385 * \returns a reference to *this
ykuroda 0:13a5d365ba16 386 *
ykuroda 0:13a5d365ba16 387 * Example: \include TutorialLinAlgComputeTwice.cpp
ykuroda 0:13a5d365ba16 388 * Output: \verbinclude TutorialLinAlgComputeTwice.out
ykuroda 0:13a5d365ba16 389 */
ykuroda 0:13a5d365ba16 390 template<typename MatrixType, int _UpLo>
ykuroda 0:13a5d365ba16 391 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
ykuroda 0:13a5d365ba16 392 {
ykuroda 0:13a5d365ba16 393 check_template_parameters();
ykuroda 0:13a5d365ba16 394
ykuroda 0:13a5d365ba16 395 eigen_assert(a.rows()==a.cols());
ykuroda 0:13a5d365ba16 396 const Index size = a.rows();
ykuroda 0:13a5d365ba16 397 m_matrix.resize(size, size);
ykuroda 0:13a5d365ba16 398 m_matrix = a;
ykuroda 0:13a5d365ba16 399
ykuroda 0:13a5d365ba16 400 m_isInitialized = true;
ykuroda 0:13a5d365ba16 401 bool ok = Traits::inplace_decomposition(m_matrix);
ykuroda 0:13a5d365ba16 402 m_info = ok ? Success : NumericalIssue;
ykuroda 0:13a5d365ba16 403
ykuroda 0:13a5d365ba16 404 return *this;
ykuroda 0:13a5d365ba16 405 }
ykuroda 0:13a5d365ba16 406
ykuroda 0:13a5d365ba16 407 /** Performs a rank one update (or dowdate) of the current decomposition.
ykuroda 0:13a5d365ba16 408 * If A = LL^* before the rank one update,
ykuroda 0:13a5d365ba16 409 * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
ykuroda 0:13a5d365ba16 410 * of same dimension.
ykuroda 0:13a5d365ba16 411 */
ykuroda 0:13a5d365ba16 412 template<typename _MatrixType, int _UpLo>
ykuroda 0:13a5d365ba16 413 template<typename VectorType>
ykuroda 0:13a5d365ba16 414 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
ykuroda 0:13a5d365ba16 415 {
ykuroda 0:13a5d365ba16 416 EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
ykuroda 0:13a5d365ba16 417 eigen_assert(v.size()==m_matrix.cols());
ykuroda 0:13a5d365ba16 418 eigen_assert(m_isInitialized);
ykuroda 0:13a5d365ba16 419 if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
ykuroda 0:13a5d365ba16 420 m_info = NumericalIssue;
ykuroda 0:13a5d365ba16 421 else
ykuroda 0:13a5d365ba16 422 m_info = Success;
ykuroda 0:13a5d365ba16 423
ykuroda 0:13a5d365ba16 424 return *this;
ykuroda 0:13a5d365ba16 425 }
ykuroda 0:13a5d365ba16 426
ykuroda 0:13a5d365ba16 427 namespace internal {
ykuroda 0:13a5d365ba16 428 template<typename _MatrixType, int UpLo, typename Rhs>
ykuroda 0:13a5d365ba16 429 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
ykuroda 0:13a5d365ba16 430 : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
ykuroda 0:13a5d365ba16 431 {
ykuroda 0:13a5d365ba16 432 typedef LLT<_MatrixType,UpLo> LLTType;
ykuroda 0:13a5d365ba16 433 EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
ykuroda 0:13a5d365ba16 434
ykuroda 0:13a5d365ba16 435 template<typename Dest> void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 436 {
ykuroda 0:13a5d365ba16 437 dst = rhs();
ykuroda 0:13a5d365ba16 438 dec().solveInPlace(dst);
ykuroda 0:13a5d365ba16 439 }
ykuroda 0:13a5d365ba16 440 };
ykuroda 0:13a5d365ba16 441 }
ykuroda 0:13a5d365ba16 442
ykuroda 0:13a5d365ba16 443 /** \internal use x = llt_object.solve(x);
ykuroda 0:13a5d365ba16 444 *
ykuroda 0:13a5d365ba16 445 * This is the \em in-place version of solve().
ykuroda 0:13a5d365ba16 446 *
ykuroda 0:13a5d365ba16 447 * \param bAndX represents both the right-hand side matrix b and result x.
ykuroda 0:13a5d365ba16 448 *
ykuroda 0:13a5d365ba16 449 * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
ykuroda 0:13a5d365ba16 450 *
ykuroda 0:13a5d365ba16 451 * This version avoids a copy when the right hand side matrix b is not
ykuroda 0:13a5d365ba16 452 * needed anymore.
ykuroda 0:13a5d365ba16 453 *
ykuroda 0:13a5d365ba16 454 * \sa LLT::solve(), MatrixBase::llt()
ykuroda 0:13a5d365ba16 455 */
ykuroda 0:13a5d365ba16 456 template<typename MatrixType, int _UpLo>
ykuroda 0:13a5d365ba16 457 template<typename Derived>
ykuroda 0:13a5d365ba16 458 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
ykuroda 0:13a5d365ba16 459 {
ykuroda 0:13a5d365ba16 460 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 461 eigen_assert(m_matrix.rows()==bAndX.rows());
ykuroda 0:13a5d365ba16 462 matrixL().solveInPlace(bAndX);
ykuroda 0:13a5d365ba16 463 matrixU().solveInPlace(bAndX);
ykuroda 0:13a5d365ba16 464 }
ykuroda 0:13a5d365ba16 465
ykuroda 0:13a5d365ba16 466 /** \returns the matrix represented by the decomposition,
ykuroda 0:13a5d365ba16 467 * i.e., it returns the product: L L^*.
ykuroda 0:13a5d365ba16 468 * This function is provided for debug purpose. */
ykuroda 0:13a5d365ba16 469 template<typename MatrixType, int _UpLo>
ykuroda 0:13a5d365ba16 470 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
ykuroda 0:13a5d365ba16 471 {
ykuroda 0:13a5d365ba16 472 eigen_assert(m_isInitialized && "LLT is not initialized.");
ykuroda 0:13a5d365ba16 473 return matrixL() * matrixL().adjoint().toDenseMatrix();
ykuroda 0:13a5d365ba16 474 }
ykuroda 0:13a5d365ba16 475
ykuroda 0:13a5d365ba16 476 /** \cholesky_module
ykuroda 0:13a5d365ba16 477 * \returns the LLT decomposition of \c *this
ykuroda 0:13a5d365ba16 478 */
ykuroda 0:13a5d365ba16 479 template<typename Derived>
ykuroda 0:13a5d365ba16 480 inline const LLT<typename MatrixBase<Derived>::PlainObject>
ykuroda 0:13a5d365ba16 481 MatrixBase<Derived>::llt() const
ykuroda 0:13a5d365ba16 482 {
ykuroda 0:13a5d365ba16 483 return LLT<PlainObject>(derived());
ykuroda 0:13a5d365ba16 484 }
ykuroda 0:13a5d365ba16 485
ykuroda 0:13a5d365ba16 486 /** \cholesky_module
ykuroda 0:13a5d365ba16 487 * \returns the LLT decomposition of \c *this
ykuroda 0:13a5d365ba16 488 */
ykuroda 0:13a5d365ba16 489 template<typename MatrixType, unsigned int UpLo>
ykuroda 0:13a5d365ba16 490 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
ykuroda 0:13a5d365ba16 491 SelfAdjointView<MatrixType, UpLo>::llt() const
ykuroda 0:13a5d365ba16 492 {
ykuroda 0:13a5d365ba16 493 return LLT<PlainObject,UpLo>(m_matrix);
ykuroda 0:13a5d365ba16 494 }
ykuroda 0:13a5d365ba16 495
ykuroda 0:13a5d365ba16 496 } // end namespace Eigen
ykuroda 0:13a5d365ba16 497
ykuroda 0:13a5d365ba16 498 #endif // EIGEN_LLT_H