Eigen  3.2.91
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 Eigen::Index Index;
63  typedef typename MatrixType::StorageIndex StorageIndex;
64 
65  enum {
66  PacketSize = internal::packet_traits<Scalar>::size,
67  AlignmentMask = int(PacketSize)-1,
68  UpLo = _UpLo
69  };
70 
71  typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
72 
79  LLT() : m_matrix(), m_isInitialized(false) {}
80 
87  explicit LLT(Index size) : m_matrix(size, size),
88  m_isInitialized(false) {}
89 
90  explicit LLT(const MatrixType& matrix)
91  : m_matrix(matrix.rows(), matrix.cols()),
92  m_isInitialized(false)
93  {
94  compute(matrix);
95  }
96 
98  inline typename Traits::MatrixU matrixU() const
99  {
100  eigen_assert(m_isInitialized && "LLT is not initialized.");
101  return Traits::getU(m_matrix);
102  }
103 
105  inline typename Traits::MatrixL matrixL() const
106  {
107  eigen_assert(m_isInitialized && "LLT is not initialized.");
108  return Traits::getL(m_matrix);
109  }
110 
121  template<typename Rhs>
122  inline const Solve<LLT, Rhs>
123  solve(const MatrixBase<Rhs>& b) const
124  {
125  eigen_assert(m_isInitialized && "LLT is not initialized.");
126  eigen_assert(m_matrix.rows()==b.rows()
127  && "LLT::solve(): invalid number of rows of the right hand side matrix b");
128  return Solve<LLT, Rhs>(*this, b.derived());
129  }
130 
131  template<typename Derived>
132  void solveInPlace(MatrixBase<Derived> &bAndX) const;
133 
134  LLT& compute(const MatrixType& matrix);
135 
140  inline const MatrixType& matrixLLT() const
141  {
142  eigen_assert(m_isInitialized && "LLT is not initialized.");
143  return m_matrix;
144  }
145 
146  MatrixType reconstructedMatrix() const;
147 
148 
155  {
156  eigen_assert(m_isInitialized && "LLT is not initialized.");
157  return m_info;
158  }
159 
160  inline Index rows() const { return m_matrix.rows(); }
161  inline Index cols() const { return m_matrix.cols(); }
162 
163  template<typename VectorType>
164  LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
165 
166  #ifndef EIGEN_PARSED_BY_DOXYGEN
167  template<typename RhsType, typename DstType>
168  EIGEN_DEVICE_FUNC
169  void _solve_impl(const RhsType &rhs, DstType &dst) const;
170  #endif
171 
172  protected:
173 
174  static void check_template_parameters()
175  {
176  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
177  }
178 
183  MatrixType m_matrix;
184  bool m_isInitialized;
185  ComputationInfo m_info;
186 };
187 
188 namespace internal {
189 
190 template<typename Scalar, int UpLo> struct llt_inplace;
191 
192 template<typename MatrixType, typename VectorType>
193 static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
194 {
195  using std::sqrt;
196  typedef typename MatrixType::Scalar Scalar;
197  typedef typename MatrixType::RealScalar RealScalar;
198  typedef typename MatrixType::ColXpr ColXpr;
199  typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
200  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
201  typedef Matrix<Scalar,Dynamic,1> TempVectorType;
202  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
203 
204  Index n = mat.cols();
205  eigen_assert(mat.rows()==n && vec.size()==n);
206 
207  TempVectorType temp;
208 
209  if(sigma>0)
210  {
211  // This version is based on Givens rotations.
212  // It is faster than the other one below, but only works for updates,
213  // i.e., for sigma > 0
214  temp = sqrt(sigma) * vec;
215 
216  for(Index i=0; i<n; ++i)
217  {
218  JacobiRotation<Scalar> g;
219  g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
220 
221  Index rs = n-i-1;
222  if(rs>0)
223  {
224  ColXprSegment x(mat.col(i).tail(rs));
225  TempVecSegment y(temp.tail(rs));
226  apply_rotation_in_the_plane(x, y, g);
227  }
228  }
229  }
230  else
231  {
232  temp = vec;
233  RealScalar beta = 1;
234  for(Index j=0; j<n; ++j)
235  {
236  RealScalar Ljj = numext::real(mat.coeff(j,j));
237  RealScalar dj = numext::abs2(Ljj);
238  Scalar wj = temp.coeff(j);
239  RealScalar swj2 = sigma*numext::abs2(wj);
240  RealScalar gamma = dj*beta + swj2;
241 
242  RealScalar x = dj + swj2/beta;
243  if (x<=RealScalar(0))
244  return j;
245  RealScalar nLjj = sqrt(x);
246  mat.coeffRef(j,j) = nLjj;
247  beta += swj2/dj;
248 
249  // Update the terms of L
250  Index rs = n-j-1;
251  if(rs)
252  {
253  temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
254  if(gamma != 0)
255  mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
256  }
257  }
258  }
259  return -1;
260 }
261 
262 template<typename Scalar> struct llt_inplace<Scalar, Lower>
263 {
264  typedef typename NumTraits<Scalar>::Real RealScalar;
265  template<typename MatrixType>
266  static Index unblocked(MatrixType& mat)
267  {
268  using std::sqrt;
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 Index blocked(MatrixType& m)
293  {
294  eigen_assert(m.rows()==m.cols());
295  Index size = m.rows();
296  if(size<32)
297  return unblocked(m);
298 
299  Index blockSize = size/8;
300  blockSize = (blockSize/16)*16;
301  blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
302 
303  for (Index k=0; k<size; k+=blockSize)
304  {
305  // partition the matrix:
306  // A00 | - | -
307  // lu = A10 | A11 | -
308  // A20 | A21 | A22
309  Index bs = (std::min)(blockSize, size-k);
310  Index rs = size - k - bs;
311  Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
312  Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
313  Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
314 
315  Index ret;
316  if((ret=unblocked(A11))>=0) return k+ret;
317  if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
318  if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
319  }
320  return -1;
321  }
322 
323  template<typename MatrixType, typename VectorType>
324  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
325  {
326  return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
327  }
328 };
329 
330 template<typename Scalar> struct llt_inplace<Scalar, Upper>
331 {
332  typedef typename NumTraits<Scalar>::Real RealScalar;
333 
334  template<typename MatrixType>
335  static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat)
336  {
337  Transpose<MatrixType> matt(mat);
338  return llt_inplace<Scalar, Lower>::unblocked(matt);
339  }
340  template<typename MatrixType>
341  static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat)
342  {
343  Transpose<MatrixType> matt(mat);
344  return llt_inplace<Scalar, Lower>::blocked(matt);
345  }
346  template<typename MatrixType, typename VectorType>
347  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
348  {
349  Transpose<MatrixType> matt(mat);
350  return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
351  }
352 };
353 
354 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
355 {
356  typedef const TriangularView<const MatrixType, Lower> MatrixL;
357  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
358  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
359  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
360  static bool inplace_decomposition(MatrixType& m)
361  { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
362 };
363 
364 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
365 {
366  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
367  typedef const TriangularView<const MatrixType, Upper> MatrixU;
368  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
369  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
370  static bool inplace_decomposition(MatrixType& m)
371  { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
372 };
373 
374 } // end namespace internal
375 
383 template<typename MatrixType, int _UpLo>
385 {
386  check_template_parameters();
387 
388  eigen_assert(a.rows()==a.cols());
389  const Index size = a.rows();
390  m_matrix.resize(size, size);
391  m_matrix = a;
392 
393  m_isInitialized = true;
394  bool ok = Traits::inplace_decomposition(m_matrix);
395  m_info = ok ? Success : NumericalIssue;
396 
397  return *this;
398 }
399 
405 template<typename _MatrixType, int _UpLo>
406 template<typename VectorType>
407 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
408 {
409  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
410  eigen_assert(v.size()==m_matrix.cols());
411  eigen_assert(m_isInitialized);
412  if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
413  m_info = NumericalIssue;
414  else
415  m_info = Success;
416 
417  return *this;
418 }
419 
420 #ifndef EIGEN_PARSED_BY_DOXYGEN
421 template<typename _MatrixType,int _UpLo>
422 template<typename RhsType, typename DstType>
423 void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
424 {
425  dst = rhs;
426  solveInPlace(dst);
427 }
428 #endif
429 
443 template<typename MatrixType, int _UpLo>
444 template<typename Derived>
445 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
446 {
447  eigen_assert(m_isInitialized && "LLT is not initialized.");
448  eigen_assert(m_matrix.rows()==bAndX.rows());
449  matrixL().solveInPlace(bAndX);
450  matrixU().solveInPlace(bAndX);
451 }
452 
456 template<typename MatrixType, int _UpLo>
458 {
459  eigen_assert(m_isInitialized && "LLT is not initialized.");
460  return matrixL() * matrixL().adjoint().toDenseMatrix();
461 }
462 
463 #ifndef __CUDACC__
464 
468 template<typename Derived>
471 {
472  return LLT<PlainObject>(derived());
473 }
474 
479 template<typename MatrixType, unsigned int UpLo>
482 {
483  return LLT<PlainObject,UpLo>(m_matrix);
484 }
485 #endif // __CUDACC__
486 
487 } // end namespace Eigen
488 
489 #endif // EIGEN_LLT_H
const LLT< PlainObject, UpLo > llt() const
Definition: LLT.h:481
MatrixType reconstructedMatrix() const
Definition: LLT.h:457
Definition: Constants.h:196
Traits::MatrixU matrixU() const
Definition: LLT.h:98
LLT(Index size)
Default Constructor with memory preallocation.
Definition: LLT.h:87
Definition: LDLT.h:16
const MatrixType & matrixLLT() const
Definition: LLT.h:140
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
const Solve< LLT, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: LLT.h:123
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: LLT.h:154
Definition: Constants.h:198
Eigen::Index Index
Definition: LLT.h:62
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:50
Definition: Constants.h:426
Definition: Constants.h:424
Definition: Eigen_Colamd.h:54
LLT & compute(const MatrixType &matrix)
Definition: LLT.h:384
const LLT< PlainObject > llt() const
Definition: LLT.h:470
Pseudo expression representing a solving operation.
Definition: Solve.h:63
LLT()
Default Constructor.
Definition: LLT.h:79
Traits::MatrixL matrixL() const
Definition: LLT.h:105
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48