00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
00064 #define EIGEN_SIMPLICIAL_CHOLESKY_H
00065
00066 namespace Eigen {
00067
00068 enum SimplicialCholeskyMode {
00069 SimplicialCholeskyLLT,
00070 SimplicialCholeskyLDLT
00071 };
00072
00085 template<typename Derived>
00086 class SimplicialCholeskyBase : internal::noncopyable
00087 {
00088 public:
00089 typedef typename internal::traits<Derived>::MatrixType MatrixType;
00090 enum { UpLo = internal::traits<Derived>::UpLo };
00091 typedef typename MatrixType::Scalar Scalar;
00092 typedef typename MatrixType::RealScalar RealScalar;
00093 typedef typename MatrixType::Index Index;
00094 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00095 typedef Matrix<Scalar,Dynamic,1> VectorType;
00096
00097 public:
00098
00100 SimplicialCholeskyBase()
00101 : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
00102 {}
00103
00104 SimplicialCholeskyBase(const MatrixType& matrix)
00105 : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
00106 {
00107 derived().compute(matrix);
00108 }
00109
00110 ~SimplicialCholeskyBase()
00111 {
00112 }
00113
00114 Derived& derived() { return *static_cast<Derived*>(this); }
00115 const Derived& derived() const { return *static_cast<const Derived*>(this); }
00116
00117 inline Index cols() const { return m_matrix.cols(); }
00118 inline Index rows() const { return m_matrix.rows(); }
00119
00125 ComputationInfo info() const
00126 {
00127 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00128 return m_info;
00129 }
00130
00135 template<typename Rhs>
00136 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
00137 solve(const MatrixBase<Rhs>& b) const
00138 {
00139 eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
00140 eigen_assert(rows()==b.rows()
00141 && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
00142 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
00143 }
00144
00149 template<typename Rhs>
00150 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
00151 solve(const SparseMatrixBase<Rhs>& b) const
00152 {
00153 eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
00154 eigen_assert(rows()==b.rows()
00155 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
00156 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
00157 }
00158
00161 const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
00162 { return m_P; }
00163
00166 const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
00167 { return m_Pinv; }
00168
00178 Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
00179 {
00180 m_shiftOffset = offset;
00181 m_shiftScale = scale;
00182 return derived();
00183 }
00184
00185 #ifndef EIGEN_PARSED_BY_DOXYGEN
00186
00187 template<typename Stream>
00188 void dumpMemory(Stream& s)
00189 {
00190 int total = 0;
00191 s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
00192 s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
00193 s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00194 s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00195 s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00196 s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00197 s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
00198 }
00199
00201 template<typename Rhs,typename Dest>
00202 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00203 {
00204 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00205 eigen_assert(m_matrix.rows()==b.rows());
00206
00207 if(m_info!=Success)
00208 return;
00209
00210 if(m_P.size()>0)
00211 dest = m_Pinv * b;
00212 else
00213 dest = b;
00214
00215 if(m_matrix.nonZeros()>0)
00216 derived().matrixL().solveInPlace(dest);
00217
00218 if(m_diag.size()>0)
00219 dest = m_diag.asDiagonal().inverse() * dest;
00220
00221 if (m_matrix.nonZeros()>0)
00222 derived().matrixU().solveInPlace(dest);
00223
00224 if(m_P.size()>0)
00225 dest = m_P * dest;
00226 }
00227
00229 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00230 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00231 {
00232 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00233 eigen_assert(m_matrix.rows()==b.rows());
00234
00235
00236 static const int NbColsAtOnce = 4;
00237 int rhsCols = b.cols();
00238 int size = b.rows();
00239 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
00240 for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00241 {
00242 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00243 tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
00244 tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
00245 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
00246 }
00247 }
00248
00249 #endif // EIGEN_PARSED_BY_DOXYGEN
00250
00251 protected:
00252
00254 template<bool DoLDLT>
00255 void compute(const MatrixType& matrix)
00256 {
00257 eigen_assert(matrix.rows()==matrix.cols());
00258 Index size = matrix.cols();
00259 CholMatrixType ap(size,size);
00260 ordering(matrix, ap);
00261 analyzePattern_preordered(ap, DoLDLT);
00262 factorize_preordered<DoLDLT>(ap);
00263 }
00264
00265 template<bool DoLDLT>
00266 void factorize(const MatrixType& a)
00267 {
00268 eigen_assert(a.rows()==a.cols());
00269 int size = a.cols();
00270 CholMatrixType ap(size,size);
00271 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
00272 factorize_preordered<DoLDLT>(ap);
00273 }
00274
00275 template<bool DoLDLT>
00276 void factorize_preordered(const CholMatrixType& a);
00277
00278 void analyzePattern(const MatrixType& a, bool doLDLT)
00279 {
00280 eigen_assert(a.rows()==a.cols());
00281 int size = a.cols();
00282 CholMatrixType ap(size,size);
00283 ordering(a, ap);
00284 analyzePattern_preordered(ap,doLDLT);
00285 }
00286 void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
00287
00288 void ordering(const MatrixType& a, CholMatrixType& ap);
00289
00291 struct keep_diag {
00292 inline bool operator() (const Index& row, const Index& col, const Scalar&) const
00293 {
00294 return row!=col;
00295 }
00296 };
00297
00298 mutable ComputationInfo m_info;
00299 bool m_isInitialized;
00300 bool m_factorizationIsOk;
00301 bool m_analysisIsOk;
00302
00303 CholMatrixType m_matrix;
00304 VectorType m_diag;
00305 VectorXi m_parent;
00306 VectorXi m_nonZerosPerCol;
00307 PermutationMatrix<Dynamic,Dynamic,Index> m_P;
00308 PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;
00309
00310 RealScalar m_shiftOffset;
00311 RealScalar m_shiftScale;
00312 };
00313
00314 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
00315 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
00316 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
00317
00318 namespace internal {
00319
00320 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
00321 {
00322 typedef _MatrixType MatrixType;
00323 enum { UpLo = _UpLo };
00324 typedef typename MatrixType::Scalar Scalar;
00325 typedef typename MatrixType::Index Index;
00326 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
00327 typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
00328 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
00329 static inline MatrixL getL(const MatrixType& m) { return m; }
00330 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00331 };
00332
00333 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
00334 {
00335 typedef _MatrixType MatrixType;
00336 enum { UpLo = _UpLo };
00337 typedef typename MatrixType::Scalar Scalar;
00338 typedef typename MatrixType::Index Index;
00339 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
00340 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
00341 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
00342 static inline MatrixL getL(const MatrixType& m) { return m; }
00343 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00344 };
00345
00346 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
00347 {
00348 typedef _MatrixType MatrixType;
00349 enum { UpLo = _UpLo };
00350 };
00351
00352 }
00353
00368 template<typename _MatrixType, int _UpLo>
00369 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
00370 {
00371 public:
00372 typedef _MatrixType MatrixType;
00373 enum { UpLo = _UpLo };
00374 typedef SimplicialCholeskyBase<SimplicialLLT> Base;
00375 typedef typename MatrixType::Scalar Scalar;
00376 typedef typename MatrixType::RealScalar RealScalar;
00377 typedef typename MatrixType::Index Index;
00378 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00379 typedef Matrix<Scalar,Dynamic,1> VectorType;
00380 typedef internal::traits<SimplicialLLT> Traits;
00381 typedef typename Traits::MatrixL MatrixL;
00382 typedef typename Traits::MatrixU MatrixU;
00383 public:
00385 SimplicialLLT() : Base() {}
00387 SimplicialLLT(const MatrixType& matrix)
00388 : Base(matrix) {}
00389
00391 inline const MatrixL matrixL() const {
00392 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
00393 return Traits::getL(Base::m_matrix);
00394 }
00395
00397 inline const MatrixU matrixU() const {
00398 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
00399 return Traits::getU(Base::m_matrix);
00400 }
00401
00403 SimplicialLLT& compute(const MatrixType& matrix)
00404 {
00405 Base::template compute<false>(matrix);
00406 return *this;
00407 }
00408
00415 void analyzePattern(const MatrixType& a)
00416 {
00417 Base::analyzePattern(a, false);
00418 }
00419
00426 void factorize(const MatrixType& a)
00427 {
00428 Base::template factorize<false>(a);
00429 }
00430
00432 Scalar determinant() const
00433 {
00434 Scalar detL = Base::m_matrix.diagonal().prod();
00435 return internal::abs2(detL);
00436 }
00437 };
00438
00453 template<typename _MatrixType, int _UpLo>
00454 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
00455 {
00456 public:
00457 typedef _MatrixType MatrixType;
00458 enum { UpLo = _UpLo };
00459 typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
00460 typedef typename MatrixType::Scalar Scalar;
00461 typedef typename MatrixType::RealScalar RealScalar;
00462 typedef typename MatrixType::Index Index;
00463 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00464 typedef Matrix<Scalar,Dynamic,1> VectorType;
00465 typedef internal::traits<SimplicialLDLT> Traits;
00466 typedef typename Traits::MatrixL MatrixL;
00467 typedef typename Traits::MatrixU MatrixU;
00468 public:
00470 SimplicialLDLT() : Base() {}
00471
00473 SimplicialLDLT(const MatrixType& matrix)
00474 : Base(matrix) {}
00475
00477 inline const VectorType vectorD() const {
00478 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00479 return Base::m_diag;
00480 }
00482 inline const MatrixL matrixL() const {
00483 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00484 return Traits::getL(Base::m_matrix);
00485 }
00486
00488 inline const MatrixU matrixU() const {
00489 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00490 return Traits::getU(Base::m_matrix);
00491 }
00492
00494 SimplicialLDLT& compute(const MatrixType& matrix)
00495 {
00496 Base::template compute<true>(matrix);
00497 return *this;
00498 }
00499
00506 void analyzePattern(const MatrixType& a)
00507 {
00508 Base::analyzePattern(a, true);
00509 }
00510
00517 void factorize(const MatrixType& a)
00518 {
00519 Base::template factorize<true>(a);
00520 }
00521
00523 Scalar determinant() const
00524 {
00525 return Base::m_diag.prod();
00526 }
00527 };
00528
00535 template<typename _MatrixType, int _UpLo>
00536 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
00537 {
00538 public:
00539 typedef _MatrixType MatrixType;
00540 enum { UpLo = _UpLo };
00541 typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
00542 typedef typename MatrixType::Scalar Scalar;
00543 typedef typename MatrixType::RealScalar RealScalar;
00544 typedef typename MatrixType::Index Index;
00545 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00546 typedef Matrix<Scalar,Dynamic,1> VectorType;
00547 typedef internal::traits<SimplicialCholesky> Traits;
00548 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
00549 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
00550 public:
00551 SimplicialCholesky() : Base(), m_LDLT(true) {}
00552
00553 SimplicialCholesky(const MatrixType& matrix)
00554 : Base(), m_LDLT(true)
00555 {
00556 compute(matrix);
00557 }
00558
00559 SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
00560 {
00561 switch(mode)
00562 {
00563 case SimplicialCholeskyLLT:
00564 m_LDLT = false;
00565 break;
00566 case SimplicialCholeskyLDLT:
00567 m_LDLT = true;
00568 break;
00569 default:
00570 break;
00571 }
00572
00573 return *this;
00574 }
00575
00576 inline const VectorType vectorD() const {
00577 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
00578 return Base::m_diag;
00579 }
00580 inline const CholMatrixType rawMatrix() const {
00581 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
00582 return Base::m_matrix;
00583 }
00584
00586 SimplicialCholesky& compute(const MatrixType& matrix)
00587 {
00588 if(m_LDLT)
00589 Base::template compute<true>(matrix);
00590 else
00591 Base::template compute<false>(matrix);
00592 return *this;
00593 }
00594
00601 void analyzePattern(const MatrixType& a)
00602 {
00603 Base::analyzePattern(a, m_LDLT);
00604 }
00605
00612 void factorize(const MatrixType& a)
00613 {
00614 if(m_LDLT)
00615 Base::template factorize<true>(a);
00616 else
00617 Base::template factorize<false>(a);
00618 }
00619
00621 template<typename Rhs,typename Dest>
00622 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00623 {
00624 eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00625 eigen_assert(Base::m_matrix.rows()==b.rows());
00626
00627 if(Base::m_info!=Success)
00628 return;
00629
00630 if(Base::m_P.size()>0)
00631 dest = Base::m_Pinv * b;
00632 else
00633 dest = b;
00634
00635 if(Base::m_matrix.nonZeros()>0)
00636 {
00637 if(m_LDLT)
00638 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
00639 else
00640 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
00641 }
00642
00643 if(Base::m_diag.size()>0)
00644 dest = Base::m_diag.asDiagonal().inverse() * dest;
00645
00646 if (Base::m_matrix.nonZeros()>0)
00647 {
00648 if(m_LDLT)
00649 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
00650 else
00651 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
00652 }
00653
00654 if(Base::m_P.size()>0)
00655 dest = Base::m_P * dest;
00656 }
00657
00658 Scalar determinant() const
00659 {
00660 if(m_LDLT)
00661 {
00662 return Base::m_diag.prod();
00663 }
00664 else
00665 {
00666 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
00667 return internal::abs2(detL);
00668 }
00669 }
00670
00671 protected:
00672 bool m_LDLT;
00673 };
00674
00675 template<typename Derived>
00676 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
00677 {
00678 eigen_assert(a.rows()==a.cols());
00679 const Index size = a.rows();
00680
00681 {
00682 CholMatrixType C;
00683 C = a.template selfadjointView<UpLo>();
00684
00685
00686
00687 internal::minimum_degree_ordering(C, m_P);
00688 }
00689
00690 if(m_P.size()>0)
00691 m_Pinv = m_P.inverse();
00692 else
00693 m_Pinv.resize(0);
00694
00695 ap.resize(size,size);
00696 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
00697 }
00698
00699 template<typename Derived>
00700 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
00701 {
00702 const Index size = ap.rows();
00703 m_matrix.resize(size, size);
00704 m_parent.resize(size);
00705 m_nonZerosPerCol.resize(size);
00706
00707 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
00708
00709 for(Index k = 0; k < size; ++k)
00710 {
00711
00712 m_parent[k] = -1;
00713 tags[k] = k;
00714 m_nonZerosPerCol[k] = 0;
00715 for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
00716 {
00717 Index i = it.index();
00718 if(i < k)
00719 {
00720
00721 for(; tags[i] != k; i = m_parent[i])
00722 {
00723
00724 if (m_parent[i] == -1)
00725 m_parent[i] = k;
00726 m_nonZerosPerCol[i]++;
00727 tags[i] = k;
00728 }
00729 }
00730 }
00731 }
00732
00733
00734 Index* Lp = m_matrix.outerIndexPtr();
00735 Lp[0] = 0;
00736 for(Index k = 0; k < size; ++k)
00737 Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
00738
00739 m_matrix.resizeNonZeros(Lp[size]);
00740
00741 m_isInitialized = true;
00742 m_info = Success;
00743 m_analysisIsOk = true;
00744 m_factorizationIsOk = false;
00745 }
00746
00747
00748 template<typename Derived>
00749 template<bool DoLDLT>
00750 void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap)
00751 {
00752 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00753 eigen_assert(ap.rows()==ap.cols());
00754 const Index size = ap.rows();
00755 eigen_assert(m_parent.size()==size);
00756 eigen_assert(m_nonZerosPerCol.size()==size);
00757
00758 const Index* Lp = m_matrix.outerIndexPtr();
00759 Index* Li = m_matrix.innerIndexPtr();
00760 Scalar* Lx = m_matrix.valuePtr();
00761
00762 ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
00763 ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0);
00764 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
00765
00766 bool ok = true;
00767 m_diag.resize(DoLDLT ? size : 0);
00768
00769 for(Index k = 0; k < size; ++k)
00770 {
00771
00772 y[k] = 0.0;
00773 Index top = size;
00774 tags[k] = k;
00775 m_nonZerosPerCol[k] = 0;
00776 for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
00777 {
00778 Index i = it.index();
00779 if(i <= k)
00780 {
00781 y[i] += internal::conj(it.value());
00782 Index len;
00783 for(len = 0; tags[i] != k; i = m_parent[i])
00784 {
00785 pattern[len++] = i;
00786 tags[i] = k;
00787 }
00788 while(len > 0)
00789 pattern[--top] = pattern[--len];
00790 }
00791 }
00792
00793
00794
00795 RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset;
00796 y[k] = 0.0;
00797 for(; top < size; ++top)
00798 {
00799 Index i = pattern[top];
00800 Scalar yi = y[i];
00801 y[i] = 0.0;
00802
00803
00804 Scalar l_ki;
00805 if(DoLDLT)
00806 l_ki = yi / m_diag[i];
00807 else
00808 yi = l_ki = yi / Lx[Lp[i]];
00809
00810 Index p2 = Lp[i] + m_nonZerosPerCol[i];
00811 Index p;
00812 for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
00813 y[Li[p]] -= internal::conj(Lx[p]) * yi;
00814 d -= internal::real(l_ki * internal::conj(yi));
00815 Li[p] = k;
00816 Lx[p] = l_ki;
00817 ++m_nonZerosPerCol[i];
00818 }
00819 if(DoLDLT)
00820 {
00821 m_diag[k] = d;
00822 if(d == RealScalar(0))
00823 {
00824 ok = false;
00825 break;
00826 }
00827 }
00828 else
00829 {
00830 Index p = Lp[k] + m_nonZerosPerCol[k]++;
00831 Li[p] = k ;
00832 if(d <= RealScalar(0)) {
00833 ok = false;
00834 break;
00835 }
00836 Lx[p] = internal::sqrt(d) ;
00837 }
00838 }
00839
00840 m_info = ok ? Success : NumericalIssue;
00841 m_factorizationIsOk = true;
00842 }
00843
00844 namespace internal {
00845
00846 template<typename Derived, typename Rhs>
00847 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
00848 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
00849 {
00850 typedef SimplicialCholeskyBase<Derived> Dec;
00851 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00852
00853 template<typename Dest> void evalTo(Dest& dst) const
00854 {
00855 dec().derived()._solve(rhs(),dst);
00856 }
00857 };
00858
00859 template<typename Derived, typename Rhs>
00860 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
00861 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
00862 {
00863 typedef SimplicialCholeskyBase<Derived> Dec;
00864 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00865
00866 template<typename Dest> void evalTo(Dest& dst) const
00867 {
00868 dec().derived()._solve_sparse(rhs(),dst);
00869 }
00870 };
00871
00872 }
00873
00874 }
00875
00876 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H