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