Eigen  3.2.2
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LLT.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_LLT_H
11 #define EIGEN_LLT_H
12 
13 namespace Eigen {
14 
15 namespace internal{
16 template<typename MatrixType, int UpLo> struct LLT_Traits;
17 }
18 
46  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
47  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
48  * the strict lower part does not have to store correct values.
49  */
50 template<typename _MatrixType, int _UpLo> class LLT
51 {
52  public:
53  typedef _MatrixType MatrixType;
54  enum {
55  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
56  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
57  Options = MatrixType::Options,
58  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
59  };
60  typedef typename MatrixType::Scalar Scalar;
61  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
62  typedef typename MatrixType::Index Index;
63 
64  enum {
65  PacketSize = internal::packet_traits<Scalar>::size,
66  AlignmentMask = int(PacketSize)-1,
67  UpLo = _UpLo
68  };
69 
70  typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
71 
78  LLT() : m_matrix(), m_isInitialized(false) {}
79 
86  LLT(Index size) : m_matrix(size, size),
87  m_isInitialized(false) {}
88 
89  LLT(const MatrixType& matrix)
90  : m_matrix(matrix.rows(), matrix.cols()),
91  m_isInitialized(false)
92  {
93  compute(matrix);
94  }
95 
97  inline typename Traits::MatrixU matrixU() const
98  {
99  eigen_assert(m_isInitialized && "LLT is not initialized.");
100  return Traits::getU(m_matrix);
101  }
102 
104  inline typename Traits::MatrixL matrixL() const
105  {
106  eigen_assert(m_isInitialized && "LLT is not initialized.");
107  return Traits::getL(m_matrix);
108  }
109 
120  template<typename Rhs>
121  inline const internal::solve_retval<LLT, Rhs>
122  solve(const MatrixBase<Rhs>& b) const
123  {
124  eigen_assert(m_isInitialized && "LLT is not initialized.");
125  eigen_assert(m_matrix.rows()==b.rows()
126  && "LLT::solve(): invalid number of rows of the right hand side matrix b");
127  return internal::solve_retval<LLT, Rhs>(*this, b.derived());
128  }
129 
130  #ifdef EIGEN2_SUPPORT
131  template<typename OtherDerived, typename ResultType>
132  bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
133  {
134  *result = this->solve(b);
135  return true;
136  }
137 
138  bool isPositiveDefinite() const { return true; }
139  #endif
140 
141  template<typename Derived>
142  void solveInPlace(MatrixBase<Derived> &bAndX) const;
143 
144  LLT& compute(const MatrixType& matrix);
145 
150  inline const MatrixType& matrixLLT() const
151  {
152  eigen_assert(m_isInitialized && "LLT is not initialized.");
153  return m_matrix;
154  }
155 
156  MatrixType reconstructedMatrix() const;
157 
158 
165  {
166  eigen_assert(m_isInitialized && "LLT is not initialized.");
167  return m_info;
168  }
169 
170  inline Index rows() const { return m_matrix.rows(); }
171  inline Index cols() const { return m_matrix.cols(); }
172 
173  template<typename VectorType>
174  LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
175 
176  protected:
181  MatrixType m_matrix;
182  bool m_isInitialized;
183  ComputationInfo m_info;
184 };
185 
186 namespace internal {
187 
188 template<typename Scalar, int UpLo> struct llt_inplace;
189 
190 template<typename MatrixType, typename VectorType>
191 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
192 {
193  using std::sqrt;
194  typedef typename MatrixType::Scalar Scalar;
195  typedef typename MatrixType::RealScalar RealScalar;
196  typedef typename MatrixType::Index Index;
197  typedef typename MatrixType::ColXpr ColXpr;
198  typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
199  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
200  typedef Matrix<Scalar,Dynamic,1> TempVectorType;
201  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
202 
203  Index n = mat.cols();
204  eigen_assert(mat.rows()==n && vec.size()==n);
205 
206  TempVectorType temp;
207 
208  if(sigma>0)
209  {
210  // This version is based on Givens rotations.
211  // It is faster than the other one below, but only works for updates,
212  // i.e., for sigma > 0
213  temp = sqrt(sigma) * vec;
214 
215  for(Index i=0; i<n; ++i)
216  {
217  JacobiRotation<Scalar> g;
218  g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
219 
220  Index rs = n-i-1;
221  if(rs>0)
222  {
223  ColXprSegment x(mat.col(i).tail(rs));
224  TempVecSegment y(temp.tail(rs));
225  apply_rotation_in_the_plane(x, y, g);
226  }
227  }
228  }
229  else
230  {
231  temp = vec;
232  RealScalar beta = 1;
233  for(Index j=0; j<n; ++j)
234  {
235  RealScalar Ljj = numext::real(mat.coeff(j,j));
236  RealScalar dj = numext::abs2(Ljj);
237  Scalar wj = temp.coeff(j);
238  RealScalar swj2 = sigma*numext::abs2(wj);
239  RealScalar gamma = dj*beta + swj2;
240 
241  RealScalar x = dj + swj2/beta;
242  if (x<=RealScalar(0))
243  return j;
244  RealScalar nLjj = sqrt(x);
245  mat.coeffRef(j,j) = nLjj;
246  beta += swj2/dj;
247 
248  // Update the terms of L
249  Index rs = n-j-1;
250  if(rs)
251  {
252  temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
253  if(gamma != 0)
254  mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
255  }
256  }
257  }
258  return -1;
259 }
260 
261 template<typename Scalar> struct llt_inplace<Scalar, Lower>
262 {
263  typedef typename NumTraits<Scalar>::Real RealScalar;
264  template<typename MatrixType>
265  static typename MatrixType::Index unblocked(MatrixType& mat)
266  {
267  using std::sqrt;
268  typedef typename MatrixType::Index Index;
269 
270  eigen_assert(mat.rows()==mat.cols());
271  const Index size = mat.rows();
272  for(Index k = 0; k < size; ++k)
273  {
274  Index rs = size-k-1; // remaining size
275 
276  Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
277  Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
278  Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
279 
280  RealScalar x = numext::real(mat.coeff(k,k));
281  if (k>0) x -= A10.squaredNorm();
282  if (x<=RealScalar(0))
283  return k;
284  mat.coeffRef(k,k) = x = sqrt(x);
285  if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
286  if (rs>0) A21 *= RealScalar(1)/x;
287  }
288  return -1;
289  }
290 
291  template<typename MatrixType>
292  static typename MatrixType::Index blocked(MatrixType& m)
293  {
294  typedef typename MatrixType::Index Index;
295  eigen_assert(m.rows()==m.cols());
296  Index size = m.rows();
297  if(size<32)
298  return unblocked(m);
299 
300  Index blockSize = size/8;
301  blockSize = (blockSize/16)*16;
302  blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
303 
304  for (Index k=0; k<size; k+=blockSize)
305  {
306  // partition the matrix:
307  // A00 | - | -
308  // lu = A10 | A11 | -
309  // A20 | A21 | A22
310  Index bs = (std::min)(blockSize, size-k);
311  Index rs = size - k - bs;
312  Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
313  Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
314  Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
315 
316  Index ret;
317  if((ret=unblocked(A11))>=0) return k+ret;
318  if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
319  if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
320  }
321  return -1;
322  }
323 
324  template<typename MatrixType, typename VectorType>
325  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
326  {
327  return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
328  }
329 };
330 
331 template<typename Scalar> struct llt_inplace<Scalar, Upper>
332 {
333  typedef typename NumTraits<Scalar>::Real RealScalar;
334 
335  template<typename MatrixType>
336  static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
337  {
338  Transpose<MatrixType> matt(mat);
339  return llt_inplace<Scalar, Lower>::unblocked(matt);
340  }
341  template<typename MatrixType>
342  static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
343  {
344  Transpose<MatrixType> matt(mat);
345  return llt_inplace<Scalar, Lower>::blocked(matt);
346  }
347  template<typename MatrixType, typename VectorType>
348  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
349  {
350  Transpose<MatrixType> matt(mat);
351  return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
352  }
353 };
354 
355 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
356 {
357  typedef const TriangularView<const MatrixType, Lower> MatrixL;
358  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
359  static inline MatrixL getL(const MatrixType& m) { return m; }
360  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
361  static bool inplace_decomposition(MatrixType& m)
362  { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
363 };
364 
365 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
366 {
367  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
368  typedef const TriangularView<const MatrixType, Upper> MatrixU;
369  static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
370  static inline MatrixU getU(const MatrixType& m) { return m; }
371  static bool inplace_decomposition(MatrixType& m)
372  { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
373 };
374 
375 } // end namespace internal
376 
384 template<typename MatrixType, int _UpLo>
386 {
387  eigen_assert(a.rows()==a.cols());
388  const Index size = a.rows();
389  m_matrix.resize(size, size);
390  m_matrix = a;
391 
392  m_isInitialized = true;
393  bool ok = Traits::inplace_decomposition(m_matrix);
394  m_info = ok ? Success : NumericalIssue;
395 
396  return *this;
397 }
398 
404 template<typename _MatrixType, int _UpLo>
405 template<typename VectorType>
406 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
407 {
408  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
409  eigen_assert(v.size()==m_matrix.cols());
410  eigen_assert(m_isInitialized);
411  if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
412  m_info = NumericalIssue;
413  else
414  m_info = Success;
415 
416  return *this;
417 }
418 
419 namespace internal {
420 template<typename _MatrixType, int UpLo, typename Rhs>
421 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
422  : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
423 {
424  typedef LLT<_MatrixType,UpLo> LLTType;
425  EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
426 
427  template<typename Dest> void evalTo(Dest& dst) const
428  {
429  dst = rhs();
430  dec().solveInPlace(dst);
431  }
432 };
433 }
434 
448 template<typename MatrixType, int _UpLo>
449 template<typename Derived>
450 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
451 {
452  eigen_assert(m_isInitialized && "LLT is not initialized.");
453  eigen_assert(m_matrix.rows()==bAndX.rows());
454  matrixL().solveInPlace(bAndX);
455  matrixU().solveInPlace(bAndX);
456 }
457 
461 template<typename MatrixType, int _UpLo>
463 {
464  eigen_assert(m_isInitialized && "LLT is not initialized.");
465  return matrixL() * matrixL().adjoint().toDenseMatrix();
466 }
467 
471 template<typename Derived>
474 {
475  return LLT<PlainObject>(derived());
476 }
477 
481 template<typename MatrixType, unsigned int UpLo>
484 {
485  return LLT<PlainObject,UpLo>(m_matrix);
486 }
487 
488 } // end namespace Eigen
489 
490 #endif // EIGEN_LLT_H
MatrixType reconstructedMatrix() const
Definition: LLT.h:462
LLT()
Default Constructor.
Definition: LLT.h:78
Definition: Constants.h:378
Traits::MatrixL matrixL() const
Definition: LLT.h:104
LLT(Index size)
Default Constructor with memory preallocation.
Definition: LLT.h:86
Traits::MatrixU matrixU() const
Definition: LLT.h:97
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
LLT & compute(const MatrixType &matrix)
Definition: LLT.h:385
Definition: Constants.h:169
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:50
Definition: Constants.h:167
const LLT< PlainObject, UpLo > llt() const
Definition: LLT.h:483
Definition: Constants.h:376
const internal::solve_retval< LLT, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: LLT.h:122
ComputationInfo
Definition: Constants.h:374
const LLT< PlainObject > llt() const
Definition: LLT.h:473
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: LLT.h:164
const MatrixType & matrixLLT() const
Definition: LLT.h:150