LLT.h
Go to the documentation of this file.
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 
17 template<typename _MatrixType, int _UpLo> struct traits<LLT<_MatrixType, _UpLo> >
18  : traits<_MatrixType>
19 {
20  typedef MatrixXpr XprKind;
22  typedef int StorageIndex;
23  enum { Flags = 0 };
24 };
25 
26 template<typename MatrixType, int UpLo> struct LLT_Traits;
27 }
28 
66 template<typename _MatrixType, int _UpLo> class LLT
67  : public SolverBase<LLT<_MatrixType, _UpLo> >
68 {
69  public:
70  typedef _MatrixType MatrixType;
72  friend class SolverBase<LLT>;
73 
75  enum {
76  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
77  };
78 
79  enum {
81  AlignmentMask = int(PacketSize)-1,
82  UpLo = _UpLo
83  };
84 
86 
93  LLT() : m_matrix(), m_isInitialized(false) {}
94 
101  explicit LLT(Index size) : m_matrix(size, size),
102  m_isInitialized(false) {}
103 
104  template<typename InputType>
105  explicit LLT(const EigenBase<InputType>& matrix)
106  : m_matrix(matrix.rows(), matrix.cols()),
107  m_isInitialized(false)
108  {
109  compute(matrix.derived());
110  }
111 
119  template<typename InputType>
121  : m_matrix(matrix.derived()),
122  m_isInitialized(false)
123  {
124  compute(matrix.derived());
125  }
126 
128  inline typename Traits::MatrixU matrixU() const
129  {
130  eigen_assert(m_isInitialized && "LLT is not initialized.");
131  return Traits::getU(m_matrix);
132  }
133 
135  inline typename Traits::MatrixL matrixL() const
136  {
137  eigen_assert(m_isInitialized && "LLT is not initialized.");
138  return Traits::getL(m_matrix);
139  }
140 
141  #ifdef EIGEN_PARSED_BY_DOXYGEN
142 
152  template<typename Rhs>
153  inline const Solve<LLT, Rhs>
154  solve(const MatrixBase<Rhs>& b) const;
155  #endif
156 
157  template<typename Derived>
158  void solveInPlace(const MatrixBase<Derived> &bAndX) const;
159 
160  template<typename InputType>
162 
167  {
168  eigen_assert(m_isInitialized && "LLT is not initialized.");
169  eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
170  return internal::rcond_estimate_helper(m_l1_norm, *this);
171  }
172 
177  inline const MatrixType& matrixLLT() const
178  {
179  eigen_assert(m_isInitialized && "LLT is not initialized.");
180  return m_matrix;
181  }
182 
183  MatrixType reconstructedMatrix() const;
184 
185 
192  {
193  eigen_assert(m_isInitialized && "LLT is not initialized.");
194  return m_info;
195  }
196 
202  const LLT& adjoint() const EIGEN_NOEXCEPT { return *this; };
203 
204  inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
205  inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
206 
207  template<typename VectorType>
208  LLT & rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
209 
210  #ifndef EIGEN_PARSED_BY_DOXYGEN
211  template<typename RhsType, typename DstType>
212  void _solve_impl(const RhsType &rhs, DstType &dst) const;
213 
214  template<bool Conjugate, typename RhsType, typename DstType>
215  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
216  #endif
217 
218  protected:
219 
221  {
223  }
224 
229  MatrixType m_matrix;
233 };
234 
235 namespace internal {
236 
237 template<typename Scalar, int UpLo> struct llt_inplace;
238 
239 template<typename MatrixType, typename VectorType>
241 {
242  using std::sqrt;
243  typedef typename MatrixType::Scalar Scalar;
244  typedef typename MatrixType::RealScalar RealScalar;
245  typedef typename MatrixType::ColXpr ColXpr;
246  typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
247  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
248  typedef Matrix<Scalar,Dynamic,1> TempVectorType;
249  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
250 
251  Index n = mat.cols();
252  eigen_assert(mat.rows()==n && vec.size()==n);
253 
254  TempVectorType temp;
255 
256  if(sigma>0)
257  {
258  // This version is based on Givens rotations.
259  // It is faster than the other one below, but only works for updates,
260  // i.e., for sigma > 0
261  temp = sqrt(sigma) * vec;
262 
263  for(Index i=0; i<n; ++i)
264  {
266  g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
267 
268  Index rs = n-i-1;
269  if(rs>0)
270  {
271  ColXprSegment x(mat.col(i).tail(rs));
272  TempVecSegment y(temp.tail(rs));
274  }
275  }
276  }
277  else
278  {
279  temp = vec;
280  RealScalar beta = 1;
281  for(Index j=0; j<n; ++j)
282  {
283  RealScalar Ljj = numext::real(mat.coeff(j,j));
284  RealScalar dj = numext::abs2(Ljj);
285  Scalar wj = temp.coeff(j);
286  RealScalar swj2 = sigma*numext::abs2(wj);
287  RealScalar gamma = dj*beta + swj2;
288 
289  RealScalar x = dj + swj2/beta;
290  if (x<=RealScalar(0))
291  return j;
292  RealScalar nLjj = sqrt(x);
293  mat.coeffRef(j,j) = nLjj;
294  beta += swj2/dj;
295 
296  // Update the terms of L
297  Index rs = n-j-1;
298  if(rs)
299  {
300  temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
301  if(gamma != 0)
302  mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
303  }
304  }
305  }
306  return -1;
307 }
308 
309 template<typename Scalar> struct llt_inplace<Scalar, Lower>
310 {
312  template<typename MatrixType>
314  {
315  using std::sqrt;
316 
317  eigen_assert(mat.rows()==mat.cols());
318  const Index size = mat.rows();
319  for(Index k = 0; k < size; ++k)
320  {
321  Index rs = size-k-1; // remaining size
322 
323  Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
324  Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
325  Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
326 
327  RealScalar x = numext::real(mat.coeff(k,k));
328  if (k>0) x -= A10.squaredNorm();
329  if (x<=RealScalar(0))
330  return k;
331  mat.coeffRef(k,k) = x = sqrt(x);
332  if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
333  if (rs>0) A21 /= x;
334  }
335  return -1;
336  }
337 
338  template<typename MatrixType>
340  {
341  eigen_assert(m.rows()==m.cols());
342  Index size = m.rows();
343  if(size<32)
344  return unblocked(m);
345 
346  Index blockSize = size/8;
347  blockSize = (blockSize/16)*16;
348  blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
349 
350  for (Index k=0; k<size; k+=blockSize)
351  {
352  // partition the matrix:
353  // A00 | - | -
354  // lu = A10 | A11 | -
355  // A20 | A21 | A22
356  Index bs = (std::min)(blockSize, size-k);
357  Index rs = size - k - bs;
358  Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
359  Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
360  Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
361 
362  Index ret;
363  if((ret=unblocked(A11))>=0) return k+ret;
364  if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
365  if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,typename NumTraits<RealScalar>::Literal(-1)); // bottleneck
366  }
367  return -1;
368  }
369 
370  template<typename MatrixType, typename VectorType>
371  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
372  {
373  return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
374  }
375 };
376 
377 template<typename Scalar> struct llt_inplace<Scalar, Upper>
378 {
380 
381  template<typename MatrixType>
383  {
384  Transpose<MatrixType> matt(mat);
386  }
387  template<typename MatrixType>
389  {
390  Transpose<MatrixType> matt(mat);
392  }
393  template<typename MatrixType, typename VectorType>
394  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
395  {
396  Transpose<MatrixType> matt(mat);
397  return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
398  }
399 };
400 
401 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
402 {
405  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
406  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
409 };
410 
411 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
412 {
415  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
416  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
419 };
420 
421 } // end namespace internal
422 
430 template<typename MatrixType, int _UpLo>
431 template<typename InputType>
433 {
434  check_template_parameters();
435 
436  eigen_assert(a.rows()==a.cols());
437  const Index size = a.rows();
438  m_matrix.resize(size, size);
439  if (!internal::is_same_dense(m_matrix, a.derived()))
440  m_matrix = a.derived();
441 
442  // Compute matrix L1 norm = max abs column sum.
443  m_l1_norm = RealScalar(0);
444  // TODO move this code to SelfAdjointView
445  for (Index col = 0; col < size; ++col) {
446  RealScalar abs_col_sum;
447  if (_UpLo == Lower)
448  abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
449  else
450  abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
451  if (abs_col_sum > m_l1_norm)
452  m_l1_norm = abs_col_sum;
453  }
454 
455  m_isInitialized = true;
456  bool ok = Traits::inplace_decomposition(m_matrix);
457  m_info = ok ? Success : NumericalIssue;
458 
459  return *this;
460 }
461 
467 template<typename _MatrixType, int _UpLo>
468 template<typename VectorType>
470 {
472  eigen_assert(v.size()==m_matrix.cols());
473  eigen_assert(m_isInitialized);
475  m_info = NumericalIssue;
476  else
477  m_info = Success;
478 
479  return *this;
480 }
481 
482 #ifndef EIGEN_PARSED_BY_DOXYGEN
483 template<typename _MatrixType,int _UpLo>
484 template<typename RhsType, typename DstType>
485 void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
486 {
487  _solve_impl_transposed<true>(rhs, dst);
488 }
489 
490 template<typename _MatrixType,int _UpLo>
491 template<bool Conjugate, typename RhsType, typename DstType>
492 void LLT<_MatrixType,_UpLo>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
493 {
494  dst = rhs;
495 
496  matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst);
497  matrixU().template conjugateIf<!Conjugate>().solveInPlace(dst);
498 }
499 #endif
500 
514 template<typename MatrixType, int _UpLo>
515 template<typename Derived>
517 {
518  eigen_assert(m_isInitialized && "LLT is not initialized.");
519  eigen_assert(m_matrix.rows()==bAndX.rows());
520  matrixL().solveInPlace(bAndX);
521  matrixU().solveInPlace(bAndX);
522 }
523 
527 template<typename MatrixType, int _UpLo>
529 {
530  eigen_assert(m_isInitialized && "LLT is not initialized.");
531  return matrixL() * matrixL().adjoint().toDenseMatrix();
532 }
533 
538 template<typename Derived>
541 {
542  return LLT<PlainObject>(derived());
543 }
544 
549 template<typename MatrixType, unsigned int UpLo>
552 {
553  return LLT<PlainObject,UpLo>(m_matrix);
554 }
555 
556 } // end namespace Eigen
557 
558 #endif // EIGEN_LLT_H
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define EIGEN_GENERIC_PUBLIC_INTERFACE(Derived)
Definition: Macros.h:1264
#define max(a, b)
Definition: datatypes.h:20
#define EIGEN_STRONG_INLINE
Definition: Macros.h:917
SolverBase< LLT > Base
Definition: LLT.h:71
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: LLT.h:205
void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const
Definition: LLT.h:492
VectorBlock< Derived > SegmentReturnType
Definition: BlockMethods.h:38
float real
Definition: datatypes.h:10
Scalar * b
Definition: benchVecAdd.cpp:17
internal::traits< Derived >::Scalar Scalar
Definition: SolverBase.h:73
const TriangularView< const MatrixType, Lower > MatrixL
Definition: LLT.h:403
#define min(a, b)
Definition: datatypes.h:19
MatrixType reconstructedMatrix() const
Definition: LLT.h:528
Expression of the transpose of a matrix.
Definition: Transpose.h:52
EIGEN_DEVICE_FUNC void apply_rotation_in_the_plane(DenseBase< VectorX > &xpr_x, DenseBase< VectorY > &xpr_y, const JacobiRotation< OtherScalar > &j)
Definition: Jacobi.h:453
int n
const LLT & adjoint() const EIGEN_NOEXCEPT
Definition: LLT.h:202
LLT(Index size)
Default Constructor with memory preallocation.
Definition: LLT.h:101
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Block< Derived, internal::traits< Derived >::RowsAtCompileTime, 1, !IsRowMajor > ColXpr
Definition: BlockMethods.h:14
Rotation given by a cosine-sine pair.
static bool inplace_decomposition(MatrixType &m)
Definition: LLT.h:417
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
const TriangularView< const typename MatrixType::AdjointReturnType, Upper > MatrixU
Definition: LLT.h:404
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
LLT(const EigenBase< InputType > &matrix)
Definition: LLT.h:105
Decomposition::RealScalar rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition &dec)
Reciprocal condition number estimator.
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
void _solve_impl(const RhsType &rhs, DstType &dst) const
Definition: LLT.h:485
void g(const string &key, int i)
Definition: testBTree.cpp:41
LLT(EigenBase< InputType > &matrix)
Constructs a LLT factorization from a given matrix.
Definition: LLT.h:120
AnnoyingScalar conj(const AnnoyingScalar &x)
static EIGEN_STRONG_INLINE Index blocked(MatrixType &mat)
Definition: LLT.h:388
_MatrixType MatrixType
Definition: LLT.h:70
bool m_isInitialized
Definition: LLT.h:231
RealScalar rcond() const
Definition: LLT.h:166
EIGEN_DEVICE_FUNC void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:162
internal::LLT_Traits< MatrixType, UpLo > Traits
Definition: LLT.h:85
#define EIGEN_NOEXCEPT
Definition: Macros.h:1418
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:66
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
Array< int, Dynamic, 1 > v
NumTraits< Scalar >::Real RealScalar
Definition: LLT.h:311
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
#define EIGEN_CONSTEXPR
Definition: Macros.h:787
EIGEN_DEVICE_FUNC bool is_same_dense(const T1 &mat1, const T2 &mat2, typename enable_if< possibly_same_dense< T1, T2 >::value >::type *=0)
Definition: XprHelper.h:695
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
static void check_template_parameters()
Definition: LLT.h:220
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: LLT.h:191
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
static Index rankUpdate(MatrixType &mat, const VectorType &vec, const RealScalar &sigma)
Definition: LLT.h:394
static Index unblocked(MatrixType &mat)
Definition: LLT.h:313
static bool inplace_decomposition(MatrixType &m)
Definition: LLT.h:407
DenseIndex ret
Traits::MatrixL matrixL() const
Definition: LLT.h:135
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
static Index llt_rank_update_lower(MatrixType &mat, const VectorType &vec, const typename MatrixType::RealScalar &sigma)
Definition: LLT.h:240
ComputationInfo m_info
Definition: LLT.h:232
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
static EIGEN_STRONG_INLINE Index unblocked(MatrixType &mat)
Definition: LLT.h:382
static MatrixU getU(const MatrixType &m)
Definition: LLT.h:416
Traits::MatrixU matrixU() const
Definition: LLT.h:128
static MatrixU getU(const MatrixType &m)
Definition: LLT.h:406
static Index rankUpdate(MatrixType &mat, const VectorType &vec, const RealScalar &sigma)
Definition: LLT.h:371
Expression of a triangular part in a matrix.
LLT & rankUpdate(const VectorType &vec, const RealScalar &sigma=1)
const TriangularView< const MatrixType, Upper > MatrixU
Definition: LLT.h:414
const LLT< PlainObject > llt() const
Definition: LLT.h:540
LLT & compute(const EigenBase< InputType > &matrix)
m col(1)
static const double sigma
static Index blocked(MatrixType &m)
Definition: LLT.h:339
const TriangularView< const typename MatrixType::AdjointReturnType, Lower > MatrixL
Definition: LLT.h:413
const MatrixType & matrixLLT() const
Definition: LLT.h:177
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
RealScalar m_l1_norm
Definition: LLT.h:230
Pseudo expression representing a solving operation.
Definition: Solve.h:62
LLT()
Default Constructor.
Definition: LLT.h:93
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
ComputationInfo
Definition: Constants.h:440
static MatrixL getL(const MatrixType &m)
Definition: LLT.h:405
#define EIGEN_STATIC_ASSERT_VECTOR_ONLY(TYPE)
Definition: StaticAssert.h:142
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: LLT.h:204
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:68
EIGEN_DEVICE_FUNC bool abs2(bool x)
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
MatrixType m_matrix
Definition: LLT.h:229
const LLT< PlainObject, UpLo > llt() const
Definition: LLT.h:551
static MatrixL getL(const MatrixType &m)
Definition: LLT.h:415
std::ptrdiff_t j
void solveInPlace(const MatrixBase< Derived > &bAndX) const
Definition: LLT.h:516
NumTraits< Scalar >::Real RealScalar
Definition: LLT.h:379


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:33