SimplicialCholesky.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2010 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 /*
11 
12 NOTE: the _symbolic, and _numeric functions has been adapted from
13  the LDL library:
14 
15 LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved.
16 
17 LDL License:
18 
19  Your use or distribution of LDL or any modified version of
20  LDL implies that you agree to this License.
21 
22  This library is free software; you can redistribute it and/or
23  modify it under the terms of the GNU Lesser General Public
24  License as published by the Free Software Foundation; either
25  version 2.1 of the License, or (at your option) any later version.
26 
27  This library is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
30  Lesser General Public License for more details.
31 
32  You should have received a copy of the GNU Lesser General Public
33  License along with this library; if not, write to the Free Software
34  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
35  USA
36 
37  Permission is hereby granted to use or copy this program under the
38  terms of the GNU LGPL, provided that the Copyright, this License,
39  and the Availability of the original version is retained on all copies.
40  User documentation of any code that uses this code or any modified
41  version of this code must cite the Copyright, this License, the
42  Availability note, and "Used by permission." Permission to modify
43  the code and to distribute modified code is granted, provided the
44  Copyright, this License, and the Availability note are retained,
45  and a notice that the code was modified is included.
46  */
47 
48 #include "../Core/util/NonMPL2.h"
49 
50 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
51 #define EIGEN_SIMPLICIAL_CHOLESKY_H
52 
53 namespace Eigen {
54 
55 enum SimplicialCholeskyMode {
56  SimplicialCholeskyLLT,
57  SimplicialCholeskyLDLT
58 };
59 
75 template<typename Derived>
76 class SimplicialCholeskyBase : internal::noncopyable
77 {
78  public:
79  typedef typename internal::traits<Derived>::MatrixType MatrixType;
80  enum { UpLo = internal::traits<Derived>::UpLo };
81  typedef typename MatrixType::Scalar Scalar;
82  typedef typename MatrixType::RealScalar RealScalar;
83  typedef typename MatrixType::Index Index;
86 
87  public:
88 
91  : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
92  {}
93 
94  SimplicialCholeskyBase(const MatrixType& matrix)
95  : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
96  {
97  derived().compute(matrix);
98  }
99 
101  {
102  }
103 
104  Derived& derived() { return *static_cast<Derived*>(this); }
105  const Derived& derived() const { return *static_cast<const Derived*>(this); }
106 
107  inline Index cols() const { return m_matrix.cols(); }
108  inline Index rows() const { return m_matrix.rows(); }
109 
116  {
117  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
118  return m_info;
119  }
120 
125  template<typename Rhs>
126  inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
127  solve(const MatrixBase<Rhs>& b) const
128  {
129  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
130  eigen_assert(rows()==b.rows()
131  && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
132  return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
133  }
134 
139  template<typename Rhs>
140  inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
141  solve(const SparseMatrixBase<Rhs>& b) const
142  {
143  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
144  eigen_assert(rows()==b.rows()
145  && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
146  return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
147  }
148 
152  { return m_P; }
153 
157  { return m_Pinv; }
158 
168  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
169  {
170  m_shiftOffset = offset;
171  m_shiftScale = scale;
172  return derived();
173  }
174 
175 #ifndef EIGEN_PARSED_BY_DOXYGEN
176 
177  template<typename Stream>
178  void dumpMemory(Stream& s)
179  {
180  int total = 0;
181  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
182  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
183  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
184  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
185  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
186  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
187  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
188  }
189 
191  template<typename Rhs,typename Dest>
192  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
193  {
194  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
195  eigen_assert(m_matrix.rows()==b.rows());
196 
197  if(m_info!=Success)
198  return;
199 
200  if(m_P.size()>0)
201  dest = m_P * b;
202  else
203  dest = b;
204 
205  if(m_matrix.nonZeros()>0) // otherwise L==I
206  derived().matrixL().solveInPlace(dest);
207 
208  if(m_diag.size()>0)
209  dest = m_diag.asDiagonal().inverse() * dest;
210 
211  if (m_matrix.nonZeros()>0) // otherwise U==I
212  derived().matrixU().solveInPlace(dest);
213 
214  if(m_P.size()>0)
215  dest = m_Pinv * dest;
216  }
217 
219  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
220  void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
221  {
222  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
223  eigen_assert(m_matrix.rows()==b.rows());
224 
225  // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
226  static const int NbColsAtOnce = 4;
227  int rhsCols = b.cols();
228  int size = b.rows();
230  for(int k=0; k<rhsCols; k+=NbColsAtOnce)
231  {
232  int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
233  tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
234  tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
235  dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
236  }
237  }
238 
239 #endif // EIGEN_PARSED_BY_DOXYGEN
240 
241  protected:
242 
244  template<bool DoLDLT>
245  void compute(const MatrixType& matrix)
246  {
247  eigen_assert(matrix.rows()==matrix.cols());
248  Index size = matrix.cols();
249  CholMatrixType ap(size,size);
250  ordering(matrix, ap);
251  analyzePattern_preordered(ap, DoLDLT);
252  factorize_preordered<DoLDLT>(ap);
253  }
254 
255  template<bool DoLDLT>
256  void factorize(const MatrixType& a)
257  {
258  eigen_assert(a.rows()==a.cols());
259  int size = a.cols();
260  CholMatrixType ap(size,size);
261  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
262  factorize_preordered<DoLDLT>(ap);
263  }
264 
265  template<bool DoLDLT>
266  void factorize_preordered(const CholMatrixType& a);
267 
268  void analyzePattern(const MatrixType& a, bool doLDLT)
269  {
270  eigen_assert(a.rows()==a.cols());
271  int size = a.cols();
272  CholMatrixType ap(size,size);
273  ordering(a, ap);
274  analyzePattern_preordered(ap,doLDLT);
275  }
276  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
277 
278  void ordering(const MatrixType& a, CholMatrixType& ap);
279 
281  struct keep_diag {
282  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
283  {
284  return row!=col;
285  }
286  };
287 
288  mutable ComputationInfo m_info;
289  bool m_isInitialized;
290  bool m_factorizationIsOk;
291  bool m_analysisIsOk;
292 
293  CholMatrixType m_matrix;
294  VectorType m_diag; // the diagonal coefficients (LDLT mode)
295  VectorXi m_parent; // elimination tree
296  VectorXi m_nonZerosPerCol;
297  PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
298  PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation
299 
300  RealScalar m_shiftOffset;
301  RealScalar m_shiftScale;
302 };
303 
304 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
305 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
306 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
307 
308 namespace internal {
309 
310 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
311 {
312  typedef _MatrixType MatrixType;
313  enum { UpLo = _UpLo };
314  typedef typename MatrixType::Scalar Scalar;
315  typedef typename MatrixType::Index Index;
316  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
317  typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
318  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
319  static inline MatrixL getL(const MatrixType& m) { return m; }
320  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
321 };
322 
323 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
324 {
325  typedef _MatrixType MatrixType;
326  enum { UpLo = _UpLo };
327  typedef typename MatrixType::Scalar Scalar;
328  typedef typename MatrixType::Index Index;
329  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
330  typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
331  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
332  static inline MatrixL getL(const MatrixType& m) { return m; }
333  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
334 };
335 
336 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
337 {
338  typedef _MatrixType MatrixType;
339  enum { UpLo = _UpLo };
340 };
341 
342 }
343 
361 template<typename _MatrixType, int _UpLo>
362  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
363 {
364 public:
365  typedef _MatrixType MatrixType;
366  enum { UpLo = _UpLo };
368  typedef typename MatrixType::Scalar Scalar;
369  typedef typename MatrixType::RealScalar RealScalar;
370  typedef typename MatrixType::Index Index;
373  typedef internal::traits<SimplicialLLT> Traits;
374  typedef typename Traits::MatrixL MatrixL;
375  typedef typename Traits::MatrixU MatrixU;
376 public:
380  SimplicialLLT(const MatrixType& matrix)
381  : Base(matrix) {}
382 
384  inline const MatrixL matrixL() const {
385  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
386  return Traits::getL(Base::m_matrix);
387  }
388 
390  inline const MatrixU matrixU() const {
391  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
392  return Traits::getU(Base::m_matrix);
393  }
394 
396  SimplicialLLT& compute(const MatrixType& matrix)
397  {
398  Base::template compute<false>(matrix);
399  return *this;
400  }
401 
408  void analyzePattern(const MatrixType& a)
409  {
410  Base::analyzePattern(a, false);
411  }
412 
419  void factorize(const MatrixType& a)
420  {
421  Base::template factorize<false>(a);
422  }
423 
425  Scalar determinant() const
426  {
427  Scalar detL = Base::m_matrix.diagonal().prod();
428  return internal::abs2(detL);
429  }
430 };
431 
449 template<typename _MatrixType, int _UpLo>
450  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
451 {
452 public:
453  typedef _MatrixType MatrixType;
454  enum { UpLo = _UpLo };
456  typedef typename MatrixType::Scalar Scalar;
457  typedef typename MatrixType::RealScalar RealScalar;
458  typedef typename MatrixType::Index Index;
461  typedef internal::traits<SimplicialLDLT> Traits;
462  typedef typename Traits::MatrixL MatrixL;
463  typedef typename Traits::MatrixU MatrixU;
464 public:
467 
469  SimplicialLDLT(const MatrixType& matrix)
470  : Base(matrix) {}
471 
473  inline const VectorType vectorD() const {
474  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
475  return Base::m_diag;
476  }
478  inline const MatrixL matrixL() const {
479  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
480  return Traits::getL(Base::m_matrix);
481  }
482 
484  inline const MatrixU matrixU() const {
485  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
486  return Traits::getU(Base::m_matrix);
487  }
488 
490  SimplicialLDLT& compute(const MatrixType& matrix)
491  {
492  Base::template compute<true>(matrix);
493  return *this;
494  }
495 
502  void analyzePattern(const MatrixType& a)
503  {
504  Base::analyzePattern(a, true);
505  }
506 
513  void factorize(const MatrixType& a)
514  {
515  Base::template factorize<true>(a);
516  }
517 
519  Scalar determinant() const
520  {
521  return Base::m_diag.prod();
522  }
523 };
524 
531 template<typename _MatrixType, int _UpLo>
532  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
533 {
534 public:
535  typedef _MatrixType MatrixType;
536  enum { UpLo = _UpLo };
538  typedef typename MatrixType::Scalar Scalar;
539  typedef typename MatrixType::RealScalar RealScalar;
540  typedef typename MatrixType::Index Index;
543  typedef internal::traits<SimplicialCholesky> Traits;
544  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
545  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
546  public:
547  SimplicialCholesky() : Base(), m_LDLT(true) {}
548 
549  SimplicialCholesky(const MatrixType& matrix)
550  : Base(), m_LDLT(true)
551  {
552  compute(matrix);
553  }
554 
555  SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
556  {
557  switch(mode)
558  {
559  case SimplicialCholeskyLLT:
560  m_LDLT = false;
561  break;
562  case SimplicialCholeskyLDLT:
563  m_LDLT = true;
564  break;
565  default:
566  break;
567  }
568 
569  return *this;
570  }
571 
572  inline const VectorType vectorD() const {
573  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
574  return Base::m_diag;
575  }
576  inline const CholMatrixType rawMatrix() const {
577  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
578  return Base::m_matrix;
579  }
580 
582  SimplicialCholesky& compute(const MatrixType& matrix)
583  {
584  if(m_LDLT)
585  Base::template compute<true>(matrix);
586  else
587  Base::template compute<false>(matrix);
588  return *this;
589  }
590 
597  void analyzePattern(const MatrixType& a)
598  {
599  Base::analyzePattern(a, m_LDLT);
600  }
601 
608  void factorize(const MatrixType& a)
609  {
610  if(m_LDLT)
611  Base::template factorize<true>(a);
612  else
613  Base::template factorize<false>(a);
614  }
615 
617  template<typename Rhs,typename Dest>
618  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
619  {
620  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
621  eigen_assert(Base::m_matrix.rows()==b.rows());
622 
623  if(Base::m_info!=Success)
624  return;
625 
626  if(Base::m_P.size()>0)
627  dest = Base::m_P * b;
628  else
629  dest = b;
630 
631  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
632  {
633  if(m_LDLT)
634  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
635  else
636  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
637  }
638 
639  if(Base::m_diag.size()>0)
640  dest = Base::m_diag.asDiagonal().inverse() * dest;
641 
642  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
643  {
644  if(m_LDLT)
645  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
646  else
647  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
648  }
649 
650  if(Base::m_P.size()>0)
651  dest = Base::m_Pinv * dest;
652  }
653 
654  Scalar determinant() const
655  {
656  if(m_LDLT)
657  {
658  return Base::m_diag.prod();
659  }
660  else
661  {
662  Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
663  return internal::abs2(detL);
664  }
665  }
666 
667  protected:
668  bool m_LDLT;
669 };
670 
671 template<typename Derived>
672 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
673 {
674  eigen_assert(a.rows()==a.cols());
675  const Index size = a.rows();
676  // TODO allows to configure the permutation
677  // Note that amd compute the inverse permutation
678  {
679  CholMatrixType C;
680  C = a.template selfadjointView<UpLo>();
681  // remove diagonal entries:
682  // seems not to be needed
683  // C.prune(keep_diag());
684  internal::minimum_degree_ordering(C, m_Pinv);
685  }
686 
687  if(m_Pinv.size()>0)
688  m_P = m_Pinv.inverse();
689  else
690  m_P.resize(0);
691 
692  ap.resize(size,size);
693  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
694 }
695 
696 template<typename Derived>
697 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
698 {
699  const Index size = ap.rows();
700  m_matrix.resize(size, size);
701  m_parent.resize(size);
702  m_nonZerosPerCol.resize(size);
703 
704  ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
705 
706  for(Index k = 0; k < size; ++k)
707  {
708  /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
709  m_parent[k] = -1; /* parent of k is not yet known */
710  tags[k] = k; /* mark node k as visited */
711  m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
712  for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
713  {
714  Index i = it.index();
715  if(i < k)
716  {
717  /* follow path from i to root of etree, stop at flagged node */
718  for(; tags[i] != k; i = m_parent[i])
719  {
720  /* find parent of i if not yet determined */
721  if (m_parent[i] == -1)
722  m_parent[i] = k;
723  m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */
724  tags[i] = k; /* mark i as visited */
725  }
726  }
727  }
728  }
729 
730  /* construct Lp index array from m_nonZerosPerCol column counts */
731  Index* Lp = m_matrix.outerIndexPtr();
732  Lp[0] = 0;
733  for(Index k = 0; k < size; ++k)
734  Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
735 
736  m_matrix.resizeNonZeros(Lp[size]);
737 
738  m_isInitialized = true;
739  m_info = Success;
740  m_analysisIsOk = true;
741  m_factorizationIsOk = false;
742 }
743 
744 
745 template<typename Derived>
746 template<bool DoLDLT>
747 void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap)
748 {
749  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
750  eigen_assert(ap.rows()==ap.cols());
751  const Index size = ap.rows();
752  eigen_assert(m_parent.size()==size);
753  eigen_assert(m_nonZerosPerCol.size()==size);
754 
755  const Index* Lp = m_matrix.outerIndexPtr();
756  Index* Li = m_matrix.innerIndexPtr();
757  Scalar* Lx = m_matrix.valuePtr();
758 
759  ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
760  ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0);
761  ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
762 
763  bool ok = true;
764  m_diag.resize(DoLDLT ? size : 0);
765 
766  for(Index k = 0; k < size; ++k)
767  {
768  // compute nonzero pattern of kth row of L, in topological order
769  y[k] = 0.0; // Y(0:k) is now all zero
770  Index top = size; // stack for pattern is empty
771  tags[k] = k; // mark node k as visited
772  m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L
773  for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
774  {
775  Index i = it.index();
776  if(i <= k)
777  {
778  y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
779  Index len;
780  for(len = 0; tags[i] != k; i = m_parent[i])
781  {
782  pattern[len++] = i; /* L(k,i) is nonzero */
783  tags[i] = k; /* mark i as visited */
784  }
785  while(len > 0)
786  pattern[--top] = pattern[--len];
787  }
788  }
789 
790  /* compute numerical values kth row of L (a sparse triangular solve) */
791 
792  RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
793  y[k] = 0.0;
794  for(; top < size; ++top)
795  {
796  Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
797  Scalar yi = y[i]; /* get and clear Y(i) */
798  y[i] = 0.0;
799 
800  /* the nonzero entry L(k,i) */
801  Scalar l_ki;
802  if(DoLDLT)
803  l_ki = yi / m_diag[i];
804  else
805  yi = l_ki = yi / Lx[Lp[i]];
806 
807  Index p2 = Lp[i] + m_nonZerosPerCol[i];
808  Index p;
809  for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
810  y[Li[p]] -= internal::conj(Lx[p]) * yi;
811  d -= internal::real(l_ki * internal::conj(yi));
812  Li[p] = k; /* store L(k,i) in column form of L */
813  Lx[p] = l_ki;
814  ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
815  }
816  if(DoLDLT)
817  {
818  m_diag[k] = d;
819  if(d == RealScalar(0))
820  {
821  ok = false; /* failure, D(k,k) is zero */
822  break;
823  }
824  }
825  else
826  {
827  Index p = Lp[k] + m_nonZerosPerCol[k]++;
828  Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */
829  if(d <= RealScalar(0)) {
830  ok = false; /* failure, matrix is not positive definite */
831  break;
832  }
833  Lx[p] = internal::sqrt(d) ;
834  }
835  }
836 
837  m_info = ok ? Success : NumericalIssue;
838  m_factorizationIsOk = true;
839 }
840 
841 namespace internal {
842 
843 template<typename Derived, typename Rhs>
844 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
845  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
846 {
847  typedef SimplicialCholeskyBase<Derived> Dec;
848  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
849 
850  template<typename Dest> void evalTo(Dest& dst) const
851  {
852  dec().derived()._solve(rhs(),dst);
853  }
854 };
855 
856 template<typename Derived, typename Rhs>
857 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
858  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
859 {
860  typedef SimplicialCholeskyBase<Derived> Dec;
861  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
862 
863  template<typename Dest> void evalTo(Dest& dst) const
864  {
865  dec().derived()._solve_sparse(rhs(),dst);
866  }
867 };
868 
869 } // end namespace internal
870 
871 } // end namespace Eigen
872 
873 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H