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 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 /*
00011 
00012 NOTE: the _symbolic, and _numeric functions has been adapted from
00013       the LDL library:
00014 
00015 LDL Copyright (c) 2005 by Timothy A. Davis.  All Rights Reserved.
00016 
00017 LDL License:
00018 
00019     Your use or distribution of LDL or any modified version of
00020     LDL implies that you agree to this License.
00021 
00022     This library is free software; you can redistribute it and/or
00023     modify it under the terms of the GNU Lesser General Public
00024     License as published by the Free Software Foundation; either
00025     version 2.1 of the License, or (at your option) any later version.
00026 
00027     This library is distributed in the hope that it will be useful,
00028     but WITHOUT ANY WARRANTY; without even the implied warranty of
00029     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00030     Lesser General Public License for more details.
00031 
00032     You should have received a copy of the GNU Lesser General Public
00033     License along with this library; if not, write to the Free Software
00034     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
00035     USA
00036 
00037     Permission is hereby granted to use or copy this program under the
00038     terms of the GNU LGPL, provided that the Copyright, this License,
00039     and the Availability of the original version is retained on all copies.
00040     User documentation of any code that uses this code or any modified
00041     version of this code must cite the Copyright, this License, the
00042     Availability note, and "Used by permission." Permission to modify
00043     the code and to distribute modified code is granted, provided the
00044     Copyright, this License, and the Availability note are retained,
00045     and a notice that the code was modified is included.
00046  */
00047 
00048 #include "../Core/util/NonMPL2.h"
00049 
00050 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
00051 #define EIGEN_SIMPLICIAL_CHOLESKY_H
00052 
00053 namespace Eigen { 
00054 
00055 enum SimplicialCholeskyMode {
00056   SimplicialCholeskyLLT,
00057   SimplicialCholeskyLDLT
00058 };
00059 
00075 template<typename Derived>
00076 class SimplicialCholeskyBase : internal::noncopyable
00077 {
00078   public:
00079     typedef typename internal::traits<Derived>::MatrixType MatrixType;
00080     enum { UpLo = internal::traits<Derived>::UpLo };
00081     typedef typename MatrixType::Scalar Scalar;
00082     typedef typename MatrixType::RealScalar RealScalar;
00083     typedef typename MatrixType::Index Index;
00084     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00085     typedef Matrix<Scalar,Dynamic,1> VectorType;
00086 
00087   public:
00088 
00090     SimplicialCholeskyBase()
00091       : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
00092     {}
00093 
00094     SimplicialCholeskyBase(const MatrixType& matrix)
00095       : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
00096     {
00097       derived().compute(matrix);
00098     }
00099 
00100     ~SimplicialCholeskyBase()
00101     {
00102     }
00103 
00104     Derived& derived() { return *static_cast<Derived*>(this); }
00105     const Derived& derived() const { return *static_cast<const Derived*>(this); }
00106     
00107     inline Index cols() const { return m_matrix.cols(); }
00108     inline Index rows() const { return m_matrix.rows(); }
00109     
00115     ComputationInfo info() const
00116     {
00117       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00118       return m_info;
00119     }
00120     
00125     template<typename Rhs>
00126     inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
00127     solve(const MatrixBase<Rhs>& b) const
00128     {
00129       eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
00130       eigen_assert(rows()==b.rows()
00131                 && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
00132       return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
00133     }
00134     
00139     template<typename Rhs>
00140     inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
00141     solve(const SparseMatrixBase<Rhs>& b) const
00142     {
00143       eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
00144       eigen_assert(rows()==b.rows()
00145                 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
00146       return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
00147     }
00148     
00151     const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
00152     { return m_P; }
00153     
00156     const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
00157     { return m_Pinv; }
00158 
00168     Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
00169     {
00170       m_shiftOffset = offset;
00171       m_shiftScale = scale;
00172       return derived();
00173     }
00174 
00175 #ifndef EIGEN_PARSED_BY_DOXYGEN
00176 
00177     template<typename Stream>
00178     void dumpMemory(Stream& s)
00179     {
00180       int total = 0;
00181       s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
00182       s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
00183       s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00184       s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00185       s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00186       s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00187       s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
00188     }
00189 
00191     template<typename Rhs,typename Dest>
00192     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00193     {
00194       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00195       eigen_assert(m_matrix.rows()==b.rows());
00196 
00197       if(m_info!=Success)
00198         return;
00199 
00200       if(m_P.size()>0)
00201         dest = m_P * b;
00202       else
00203         dest = b;
00204 
00205       if(m_matrix.nonZeros()>0) // otherwise L==I
00206         derived().matrixL().solveInPlace(dest);
00207 
00208       if(m_diag.size()>0)
00209         dest = m_diag.asDiagonal().inverse() * dest;
00210 
00211       if (m_matrix.nonZeros()>0) // otherwise U==I
00212         derived().matrixU().solveInPlace(dest);
00213 
00214       if(m_P.size()>0)
00215         dest = m_Pinv * dest;
00216     }
00217 
00219     template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00220     void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00221     {
00222       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00223       eigen_assert(m_matrix.rows()==b.rows());
00224       
00225       // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
00226       static const int NbColsAtOnce = 4;
00227       int rhsCols = b.cols();
00228       int size = b.rows();
00229       Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
00230       for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00231       {
00232         int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00233         tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
00234         tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
00235         dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
00236       }
00237     }
00238 
00239 #endif // EIGEN_PARSED_BY_DOXYGEN
00240 
00241   protected:
00242     
00244     template<bool DoLDLT>
00245     void compute(const MatrixType& matrix)
00246     {
00247       eigen_assert(matrix.rows()==matrix.cols());
00248       Index size = matrix.cols();
00249       CholMatrixType ap(size,size);
00250       ordering(matrix, ap);
00251       analyzePattern_preordered(ap, DoLDLT);
00252       factorize_preordered<DoLDLT>(ap);
00253     }
00254     
00255     template<bool DoLDLT>
00256     void factorize(const MatrixType& a)
00257     {
00258       eigen_assert(a.rows()==a.cols());
00259       int size = a.cols();
00260       CholMatrixType ap(size,size);
00261       ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
00262       factorize_preordered<DoLDLT>(ap);
00263     }
00264 
00265     template<bool DoLDLT>
00266     void factorize_preordered(const CholMatrixType& a);
00267 
00268     void analyzePattern(const MatrixType& a, bool doLDLT)
00269     {
00270       eigen_assert(a.rows()==a.cols());
00271       int size = a.cols();
00272       CholMatrixType ap(size,size);
00273       ordering(a, ap);
00274       analyzePattern_preordered(ap,doLDLT);
00275     }
00276     void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
00277     
00278     void ordering(const MatrixType& a, CholMatrixType& ap);
00279 
00281     struct keep_diag {
00282       inline bool operator() (const Index& row, const Index& col, const Scalar&) const
00283       {
00284         return row!=col;
00285       }
00286     };
00287 
00288     mutable ComputationInfo m_info;
00289     bool m_isInitialized;
00290     bool m_factorizationIsOk;
00291     bool m_analysisIsOk;
00292     
00293     CholMatrixType m_matrix;
00294     VectorType m_diag;                                // the diagonal coefficients (LDLT mode)
00295     VectorXi m_parent;                                // elimination tree
00296     VectorXi m_nonZerosPerCol;
00297     PermutationMatrix<Dynamic,Dynamic,Index> m_P;     // the permutation
00298     PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;  // the inverse permutation
00299 
00300     RealScalar m_shiftOffset;
00301     RealScalar m_shiftScale;
00302 };
00303 
00304 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
00305 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
00306 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
00307 
00308 namespace internal {
00309 
00310 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
00311 {
00312   typedef _MatrixType MatrixType;
00313   enum { UpLo = _UpLo };
00314   typedef typename MatrixType::Scalar                         Scalar;
00315   typedef typename MatrixType::Index                          Index;
00316   typedef SparseMatrix<Scalar, ColMajor, Index>               CholMatrixType;
00317   typedef SparseTriangularView<CholMatrixType, Eigen::Lower>  MatrixL;
00318   typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper>   MatrixU;
00319   static inline MatrixL getL(const MatrixType& m) { return m; }
00320   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00321 };
00322 
00323 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
00324 {
00325   typedef _MatrixType MatrixType;
00326   enum { UpLo = _UpLo };
00327   typedef typename MatrixType::Scalar                             Scalar;
00328   typedef typename MatrixType::Index                              Index;
00329   typedef SparseMatrix<Scalar, ColMajor, Index>                   CholMatrixType;
00330   typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower>  MatrixL;
00331   typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
00332   static inline MatrixL getL(const MatrixType& m) { return m; }
00333   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00334 };
00335 
00336 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
00337 {
00338   typedef _MatrixType MatrixType;
00339   enum { UpLo = _UpLo };
00340 };
00341 
00342 }
00343 
00361 template<typename _MatrixType, int _UpLo>
00362     class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
00363 {
00364 public:
00365     typedef _MatrixType MatrixType;
00366     enum { UpLo = _UpLo };
00367     typedef SimplicialCholeskyBase<SimplicialLLT> Base;
00368     typedef typename MatrixType::Scalar Scalar;
00369     typedef typename MatrixType::RealScalar RealScalar;
00370     typedef typename MatrixType::Index Index;
00371     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00372     typedef Matrix<Scalar,Dynamic,1> VectorType;
00373     typedef internal::traits<SimplicialLLT> Traits;
00374     typedef typename Traits::MatrixL  MatrixL;
00375     typedef typename Traits::MatrixU  MatrixU;
00376 public:
00378     SimplicialLLT() : Base() {}
00380     SimplicialLLT(const MatrixType& matrix)
00381         : Base(matrix) {}
00382 
00384     inline const MatrixL matrixL() const {
00385         eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
00386         return Traits::getL(Base::m_matrix);
00387     }
00388 
00390     inline const MatrixU matrixU() const {
00391         eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
00392         return Traits::getU(Base::m_matrix);
00393     }
00394     
00396     SimplicialLLT& compute(const MatrixType& matrix)
00397     {
00398       Base::template compute<false>(matrix);
00399       return *this;
00400     }
00401 
00408     void analyzePattern(const MatrixType& a)
00409     {
00410       Base::analyzePattern(a, false);
00411     }
00412 
00419     void factorize(const MatrixType& a)
00420     {
00421       Base::template factorize<false>(a);
00422     }
00423 
00425     Scalar determinant() const
00426     {
00427       Scalar detL = Base::m_matrix.diagonal().prod();
00428       return internal::abs2(detL);
00429     }
00430 };
00431 
00449 template<typename _MatrixType, int _UpLo>
00450     class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
00451 {
00452 public:
00453     typedef _MatrixType MatrixType;
00454     enum { UpLo = _UpLo };
00455     typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
00456     typedef typename MatrixType::Scalar Scalar;
00457     typedef typename MatrixType::RealScalar RealScalar;
00458     typedef typename MatrixType::Index Index;
00459     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00460     typedef Matrix<Scalar,Dynamic,1> VectorType;
00461     typedef internal::traits<SimplicialLDLT> Traits;
00462     typedef typename Traits::MatrixL  MatrixL;
00463     typedef typename Traits::MatrixU  MatrixU;
00464 public:
00466     SimplicialLDLT() : Base() {}
00467 
00469     SimplicialLDLT(const MatrixType& matrix)
00470         : Base(matrix) {}
00471 
00473     inline const VectorType vectorD() const {
00474         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00475         return Base::m_diag;
00476     }
00478     inline const MatrixL matrixL() const {
00479         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00480         return Traits::getL(Base::m_matrix);
00481     }
00482 
00484     inline const MatrixU matrixU() const {
00485         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00486         return Traits::getU(Base::m_matrix);
00487     }
00488 
00490     SimplicialLDLT& compute(const MatrixType& matrix)
00491     {
00492       Base::template compute<true>(matrix);
00493       return *this;
00494     }
00495     
00502     void analyzePattern(const MatrixType& a)
00503     {
00504       Base::analyzePattern(a, true);
00505     }
00506 
00513     void factorize(const MatrixType& a)
00514     {
00515       Base::template factorize<true>(a);
00516     }
00517 
00519     Scalar determinant() const
00520     {
00521       return Base::m_diag.prod();
00522     }
00523 };
00524 
00531 template<typename _MatrixType, int _UpLo>
00532     class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
00533 {
00534 public:
00535     typedef _MatrixType MatrixType;
00536     enum { UpLo = _UpLo };
00537     typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
00538     typedef typename MatrixType::Scalar Scalar;
00539     typedef typename MatrixType::RealScalar RealScalar;
00540     typedef typename MatrixType::Index Index;
00541     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00542     typedef Matrix<Scalar,Dynamic,1> VectorType;
00543     typedef internal::traits<SimplicialCholesky> Traits;
00544     typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
00545     typedef internal::traits<SimplicialLLT<MatrixType,UpLo>  > LLTTraits;
00546   public:
00547     SimplicialCholesky() : Base(), m_LDLT(true) {}
00548 
00549     SimplicialCholesky(const MatrixType& matrix)
00550       : Base(), m_LDLT(true)
00551     {
00552       compute(matrix);
00553     }
00554 
00555     SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
00556     {
00557       switch(mode)
00558       {
00559       case SimplicialCholeskyLLT:
00560         m_LDLT = false;
00561         break;
00562       case SimplicialCholeskyLDLT:
00563         m_LDLT = true;
00564         break;
00565       default:
00566         break;
00567       }
00568 
00569       return *this;
00570     }
00571 
00572     inline const VectorType vectorD() const {
00573         eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
00574         return Base::m_diag;
00575     }
00576     inline const CholMatrixType rawMatrix() const {
00577         eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
00578         return Base::m_matrix;
00579     }
00580     
00582     SimplicialCholesky& compute(const MatrixType& matrix)
00583     {
00584       if(m_LDLT)
00585         Base::template compute<true>(matrix);
00586       else
00587         Base::template compute<false>(matrix);
00588       return *this;
00589     }
00590 
00597     void analyzePattern(const MatrixType& a)
00598     {
00599       Base::analyzePattern(a, m_LDLT);
00600     }
00601 
00608     void factorize(const MatrixType& a)
00609     {
00610       if(m_LDLT)
00611         Base::template factorize<true>(a);
00612       else
00613         Base::template factorize<false>(a);
00614     }
00615 
00617     template<typename Rhs,typename Dest>
00618     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00619     {
00620       eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00621       eigen_assert(Base::m_matrix.rows()==b.rows());
00622 
00623       if(Base::m_info!=Success)
00624         return;
00625 
00626       if(Base::m_P.size()>0)
00627         dest = Base::m_P * b;
00628       else
00629         dest = b;
00630 
00631       if(Base::m_matrix.nonZeros()>0) // otherwise L==I
00632       {
00633         if(m_LDLT)
00634           LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
00635         else
00636           LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
00637       }
00638 
00639       if(Base::m_diag.size()>0)
00640         dest = Base::m_diag.asDiagonal().inverse() * dest;
00641 
00642       if (Base::m_matrix.nonZeros()>0) // otherwise I==I
00643       {
00644         if(m_LDLT)
00645           LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
00646         else
00647           LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
00648       }
00649 
00650       if(Base::m_P.size()>0)
00651         dest = Base::m_Pinv * dest;
00652     }
00653     
00654     Scalar determinant() const
00655     {
00656       if(m_LDLT)
00657       {
00658         return Base::m_diag.prod();
00659       }
00660       else
00661       {
00662         Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
00663         return internal::abs2(detL);
00664       }
00665     }
00666     
00667   protected:
00668     bool m_LDLT;
00669 };
00670 
00671 template<typename Derived>
00672 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
00673 {
00674   eigen_assert(a.rows()==a.cols());
00675   const Index size = a.rows();
00676   // TODO allows to configure the permutation
00677   // Note that amd compute the inverse permutation
00678   {
00679     CholMatrixType C;
00680     C = a.template selfadjointView<UpLo>();
00681     // remove diagonal entries:
00682     // seems not to be needed
00683     // C.prune(keep_diag());
00684     internal::minimum_degree_ordering(C, m_Pinv);
00685   }
00686 
00687   if(m_Pinv.size()>0)
00688     m_P = m_Pinv.inverse();
00689   else
00690     m_P.resize(0);
00691 
00692   ap.resize(size,size);
00693   ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
00694 }
00695 
00696 template<typename Derived>
00697 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
00698 {
00699   const Index size = ap.rows();
00700   m_matrix.resize(size, size);
00701   m_parent.resize(size);
00702   m_nonZerosPerCol.resize(size);
00703   
00704   ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
00705 
00706   for(Index k = 0; k < size; ++k)
00707   {
00708     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
00709     m_parent[k] = -1;             /* parent of k is not yet known */
00710     tags[k] = k;                  /* mark node k as visited */
00711     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
00712     for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
00713     {
00714       Index i = it.index();
00715       if(i < k)
00716       {
00717         /* follow path from i to root of etree, stop at flagged node */
00718         for(; tags[i] != k; i = m_parent[i])
00719         {
00720           /* find parent of i if not yet determined */
00721           if (m_parent[i] == -1)
00722             m_parent[i] = k;
00723           m_nonZerosPerCol[i]++;        /* L (k,i) is nonzero */
00724           tags[i] = k;                  /* mark i as visited */
00725         }
00726       }
00727     }
00728   }
00729   
00730   /* construct Lp index array from m_nonZerosPerCol column counts */
00731   Index* Lp = m_matrix.outerIndexPtr();
00732   Lp[0] = 0;
00733   for(Index k = 0; k < size; ++k)
00734     Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
00735 
00736   m_matrix.resizeNonZeros(Lp[size]);
00737   
00738   m_isInitialized     = true;
00739   m_info              = Success;
00740   m_analysisIsOk      = true;
00741   m_factorizationIsOk = false;
00742 }
00743 
00744 
00745 template<typename Derived>
00746 template<bool DoLDLT>
00747 void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap)
00748 {
00749   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00750   eigen_assert(ap.rows()==ap.cols());
00751   const Index size = ap.rows();
00752   eigen_assert(m_parent.size()==size);
00753   eigen_assert(m_nonZerosPerCol.size()==size);
00754 
00755   const Index* Lp = m_matrix.outerIndexPtr();
00756   Index* Li = m_matrix.innerIndexPtr();
00757   Scalar* Lx = m_matrix.valuePtr();
00758 
00759   ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
00760   ei_declare_aligned_stack_constructed_variable(Index,  pattern, size, 0);
00761   ei_declare_aligned_stack_constructed_variable(Index,  tags, size, 0);
00762   
00763   bool ok = true;
00764   m_diag.resize(DoLDLT ? size : 0);
00765   
00766   for(Index k = 0; k < size; ++k)
00767   {
00768     // compute nonzero pattern of kth row of L, in topological order
00769     y[k] = 0.0;                     // Y(0:k) is now all zero
00770     Index top = size;               // stack for pattern is empty
00771     tags[k] = k;                    // mark node k as visited
00772     m_nonZerosPerCol[k] = 0;        // count of nonzeros in column k of L
00773     for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
00774     {
00775       Index i = it.index();
00776       if(i <= k)
00777       {
00778         y[i] += internal::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
00779         Index len;
00780         for(len = 0; tags[i] != k; i = m_parent[i])
00781         {
00782           pattern[len++] = i;     /* L(k,i) is nonzero */
00783           tags[i] = k;            /* mark i as visited */
00784         }
00785         while(len > 0)
00786           pattern[--top] = pattern[--len];
00787       }
00788     }
00789 
00790     /* compute numerical values kth row of L (a sparse triangular solve) */
00791 
00792     RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset;    // get D(k,k), apply the shift function, and clear Y(k)
00793     y[k] = 0.0;
00794     for(; top < size; ++top)
00795     {
00796       Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
00797       Scalar yi = y[i];             /* get and clear Y(i) */
00798       y[i] = 0.0;
00799       
00800       /* the nonzero entry L(k,i) */
00801       Scalar l_ki;
00802       if(DoLDLT)
00803         l_ki = yi / m_diag[i];       
00804       else
00805         yi = l_ki = yi / Lx[Lp[i]];
00806       
00807       Index p2 = Lp[i] + m_nonZerosPerCol[i];
00808       Index p;
00809       for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
00810         y[Li[p]] -= internal::conj(Lx[p]) * yi;
00811       d -= internal::real(l_ki * internal::conj(yi));
00812       Li[p] = k;                          /* store L(k,i) in column form of L */
00813       Lx[p] = l_ki;
00814       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
00815     }
00816     if(DoLDLT)
00817     {
00818       m_diag[k] = d;
00819       if(d == RealScalar(0))
00820       {
00821         ok = false;                         /* failure, D(k,k) is zero */
00822         break;
00823       }
00824     }
00825     else
00826     {
00827       Index p = Lp[k] + m_nonZerosPerCol[k]++;
00828       Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
00829       if(d <= RealScalar(0)) {
00830         ok = false;              /* failure, matrix is not positive definite */
00831         break;
00832       }
00833       Lx[p] = internal::sqrt(d) ;
00834     }
00835   }
00836 
00837   m_info = ok ? Success : NumericalIssue;
00838   m_factorizationIsOk = true;
00839 }
00840 
00841 namespace internal {
00842   
00843 template<typename Derived, typename Rhs>
00844 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
00845   : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
00846 {
00847   typedef SimplicialCholeskyBase<Derived> Dec;
00848   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00849 
00850   template<typename Dest> void evalTo(Dest& dst) const
00851   {
00852     dec().derived()._solve(rhs(),dst);
00853   }
00854 };
00855 
00856 template<typename Derived, typename Rhs>
00857 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
00858   : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
00859 {
00860   typedef SimplicialCholeskyBase<Derived> Dec;
00861   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00862 
00863   template<typename Dest> void evalTo(Dest& dst) const
00864   {
00865     dec().derived()._solve_sparse(rhs(),dst);
00866   }
00867 };
00868 
00869 } // end namespace internal
00870 
00871 } // end namespace Eigen
00872 
00873 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:26:05