SimplicialCholesky.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 /*
00026 
00027 NOTE: the _symbolic, and _numeric functions has been adapted from
00028       the LDL library:
00029 
00030 LDL Copyright (c) 2005 by Timothy A. Davis.  All Rights Reserved.
00031 
00032 LDL License:
00033 
00034     Your use or distribution of LDL or any modified version of
00035     LDL implies that you agree to this License.
00036 
00037     This library is free software; you can redistribute it and/or
00038     modify it under the terms of the GNU Lesser General Public
00039     License as published by the Free Software Foundation; either
00040     version 2.1 of the License, or (at your option) any later version.
00041 
00042     This library is distributed in the hope that it will be useful,
00043     but WITHOUT ANY WARRANTY; without even the implied warranty of
00044     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00045     Lesser General Public License for more details.
00046 
00047     You should have received a copy of the GNU Lesser General Public
00048     License along with this library; if not, write to the Free Software
00049     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
00050     USA
00051 
00052     Permission is hereby granted to use or copy this program under the
00053     terms of the GNU LGPL, provided that the Copyright, this License,
00054     and the Availability of the original version is retained on all copies.
00055     User documentation of any code that uses this code or any modified
00056     version of this code must cite the Copyright, this License, the
00057     Availability note, and "Used by permission." Permission to modify
00058     the code and to distribute modified code is granted, provided the
00059     Copyright, this License, and the Availability note are retained,
00060     and a notice that the code was modified is included.
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) // otherwise L==I
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) // otherwise I==I
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       // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
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;                                // the diagonal coefficients (LDLT mode)
00305     VectorXi m_parent;                                // elimination tree
00306     VectorXi m_nonZerosPerCol;
00307     PermutationMatrix<Dynamic,Dynamic,Index> m_P;     // the permutation
00308     PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;  // the inverse permutation
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) // otherwise L==I
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) // otherwise I==I
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   // TODO allows to configure the permutation
00681   {
00682     CholMatrixType C;
00683     C = a.template selfadjointView<UpLo>();
00684     // remove diagonal entries:
00685     // seems not to be needed
00686     // C.prune(keep_diag());
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     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
00712     m_parent[k] = -1;             /* parent of k is not yet known */
00713     tags[k] = k;                  /* mark node k as visited */
00714     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
00715     for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
00716     {
00717       Index i = it.index();
00718       if(i < k)
00719       {
00720         /* follow path from i to root of etree, stop at flagged node */
00721         for(; tags[i] != k; i = m_parent[i])
00722         {
00723           /* find parent of i if not yet determined */
00724           if (m_parent[i] == -1)
00725             m_parent[i] = k;
00726           m_nonZerosPerCol[i]++;        /* L (k,i) is nonzero */
00727           tags[i] = k;                  /* mark i as visited */
00728         }
00729       }
00730     }
00731   }
00732   
00733   /* construct Lp index array from m_nonZerosPerCol column counts */
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     // compute nonzero pattern of kth row of L, in topological order
00772     y[k] = 0.0;                     // Y(0:k) is now all zero
00773     Index top = size;               // stack for pattern is empty
00774     tags[k] = k;                    // mark node k as visited
00775     m_nonZerosPerCol[k] = 0;        // count of nonzeros in column k of L
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());            /* scatter A(i,k) into Y (sum duplicates) */
00782         Index len;
00783         for(len = 0; tags[i] != k; i = m_parent[i])
00784         {
00785           pattern[len++] = i;     /* L(k,i) is nonzero */
00786           tags[i] = k;            /* mark i as visited */
00787         }
00788         while(len > 0)
00789           pattern[--top] = pattern[--len];
00790       }
00791     }
00792 
00793     /* compute numerical values kth row of L (a sparse triangular solve) */
00794 
00795     RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset;    // get D(k,k), apply the shift function, and clear Y(k)
00796     y[k] = 0.0;
00797     for(; top < size; ++top)
00798     {
00799       Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
00800       Scalar yi = y[i];             /* get and clear Y(i) */
00801       y[i] = 0.0;
00802       
00803       /* the nonzero entry L(k,i) */
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;                          /* store L(k,i) in column form of L */
00816       Lx[p] = l_ki;
00817       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
00818     }
00819     if(DoLDLT)
00820     {
00821       m_diag[k] = d;
00822       if(d == RealScalar(0))
00823       {
00824         ok = false;                         /* failure, D(k,k) is zero */
00825         break;
00826       }
00827     }
00828     else
00829     {
00830       Index p = Lp[k] + m_nonZerosPerCol[k]++;
00831       Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
00832       if(d <= RealScalar(0)) {
00833         ok = false;              /* failure, matrix is not positive definite */
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 } // end namespace internal
00873 
00874 } // end namespace Eigen
00875 
00876 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H