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 enum SimplicialCholeskyMode {
00067   SimplicialCholeskyLLt,
00068   SimplicialCholeskyLDLt
00069 };
00070 
00082 template<typename _MatrixType, int _UpLo = Lower>
00083 class SimplicialCholesky
00084 {
00085   public:
00086     typedef _MatrixType MatrixType;
00087     enum { UpLo = _UpLo };
00088     typedef typename MatrixType::Scalar Scalar;
00089     typedef typename MatrixType::RealScalar RealScalar;
00090     typedef typename MatrixType::Index Index;
00091     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00092     typedef Matrix<Scalar,MatrixType::ColsAtCompileTime,1> VectorType;
00093 
00094   public:
00095 
00096     SimplicialCholesky()
00097       : m_info(Success), m_isInitialized(false), m_LDLt(true)
00098     {}
00099 
00100     SimplicialCholesky(const MatrixType& matrix)
00101       : m_info(Success), m_isInitialized(false), m_LDLt(true)
00102     {
00103       compute(matrix);
00104     }
00105 
00106     ~SimplicialCholesky()
00107     {
00108     }
00109     
00110     inline Index cols() const { return m_matrix.cols(); }
00111     inline Index rows() const { return m_matrix.rows(); }
00112   
00113     SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
00114     {
00115       switch(mode)
00116       {
00117         case SimplicialCholeskyLLt:
00118           m_LDLt = false;
00119           break;
00120         case SimplicialCholeskyLDLt:
00121           m_LDLt = true;
00122           break;
00123         default:
00124           break;
00125       }
00126       
00127       return *this;
00128     }
00129     
00135     ComputationInfo info() const
00136     {
00137       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00138       return m_info;
00139     }
00140 
00142     SimplicialCholesky& compute(const MatrixType& matrix)
00143     {
00144       analyzePattern(matrix);
00145       factorize(matrix);
00146       return *this;
00147     }
00148     
00153     template<typename Rhs>
00154     inline const internal::solve_retval<SimplicialCholesky, Rhs>
00155     solve(const MatrixBase<Rhs>& b) const
00156     {
00157       eigen_assert(m_isInitialized && "SimplicialCholesky is not initialized.");
00158       eigen_assert(rows()==b.rows()
00159                 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
00160       return internal::solve_retval<SimplicialCholesky, Rhs>(*this, b.derived());
00161     }
00162     
00167 //     template<typename Rhs>
00168 //     inline const internal::sparse_solve_retval<CholmodDecomposition, Rhs>
00169 //     solve(const SparseMatrixBase<Rhs>& b) const
00170 //     {
00171 //       eigen_assert(m_isInitialized && "SimplicialCholesky is not initialized.");
00172 //       eigen_assert(rows()==b.rows()
00173 //                 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
00174 //       return internal::sparse_solve_retval<SimplicialCholesky, Rhs>(*this, b.derived());
00175 //     }
00176     
00183     void analyzePattern(const MatrixType& a);
00184     
00185     
00192     void factorize(const MatrixType& a);
00193     
00196     const PermutationMatrix<Dynamic>& permutationP() const
00197     { return m_P; }
00198     
00201     const PermutationMatrix<Dynamic>& permutationPinv() const
00202     { return m_Pinv; }
00203     
00204     #ifndef EIGEN_PARSED_BY_DOXYGEN
00205 
00206     template<typename Rhs,typename Dest>
00207     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00208     {
00209       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00210       eigen_assert(m_matrix.rows()==b.rows());
00211       
00212       if(m_info!=Success)
00213         return;
00214     
00215       if(m_P.size()>0)
00216         dest = m_Pinv * b;
00217       else
00218         dest = b;
00219       
00220       if(m_LDLt)
00221       {
00222         if(m_matrix.nonZeros()>0) // otherwise L==I
00223           m_matrix.template triangularView<UnitLower>().solveInPlace(dest);
00224       
00225         dest = m_diag.asDiagonal().inverse() * dest;
00226       
00227         if (m_matrix.nonZeros()>0) // otherwise L==I
00228           m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(dest);
00229       }
00230       else
00231       {
00232         if(m_matrix.nonZeros()>0) // otherwise L==I
00233           m_matrix.template triangularView<Lower>().solveInPlace(dest);
00234       
00235         if (m_matrix.nonZeros()>0) // otherwise L==I
00236           m_matrix.adjoint().template triangularView<Upper>().solveInPlace(dest);
00237       }
00238       
00239       if(m_P.size()>0)
00240         dest = m_P * dest;
00241     }
00242     
00244     /*
00245     template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
00246     void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00247     {
00248       // TODO
00249     }
00250     */
00251     #endif // EIGEN_PARSED_BY_DOXYGEN
00252     
00253     template<typename Stream>
00254     void dumpMemory(Stream& s)
00255     {
00256       int total = 0;
00257       s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
00258       s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
00259       s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00260       s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00261       s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00262       s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00263       s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
00264     }
00265 
00266   protected:
00268     struct keep_diag {
00269       inline bool operator() (const Index& row, const Index& col, const Scalar&) const
00270       {
00271         return row!=col;
00272       }
00273     };
00274 
00275     mutable ComputationInfo m_info;
00276     bool m_isInitialized;
00277     bool m_factorizationIsOk;
00278     bool m_analysisIsOk;
00279     bool m_LDLt;
00280     
00281     CholMatrixType m_matrix;
00282     VectorType m_diag;                  // the diagonal coefficients in case of a LDLt decomposition
00283     VectorXi m_parent;                  // elimination tree
00284     VectorXi m_nonZerosPerCol;
00285     PermutationMatrix<Dynamic> m_P;     // the permutation
00286     PermutationMatrix<Dynamic> m_Pinv;  // the inverse permutation
00287 };
00288 
00289 template<typename _MatrixType, int _UpLo>
00290 void SimplicialCholesky<_MatrixType,_UpLo>::analyzePattern(const MatrixType& a)
00291 {
00292   eigen_assert(a.rows()==a.cols());
00293   const Index size = a.rows();
00294   m_matrix.resize(size, size);
00295   m_parent.resize(size);
00296   m_nonZerosPerCol.resize(size);
00297   
00298   ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
00299   
00300   // TODO allows to configure the permutation
00301   {
00302     CholMatrixType C;
00303     C = a.template selfadjointView<UpLo>();
00304     // remove diagonal entries:
00305     C.prune(keep_diag());
00306     internal::minimum_degree_ordering(C, m_P);
00307   }
00308   
00309   if(m_P.size()>0)
00310     m_Pinv  = m_P.inverse();
00311   else
00312     m_Pinv.resize(0);
00313   
00314   SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
00315   ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
00316   
00317   for(Index k = 0; k < size; ++k)
00318   {
00319     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
00320     m_parent[k] = -1;             /* parent of k is not yet known */
00321     tags[k] = k;                  /* mark node k as visited */
00322     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
00323     for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
00324     {
00325       Index i = it.index();
00326       if(i < k)
00327       {
00328         /* follow path from i to root of etree, stop at flagged node */
00329         for(; tags[i] != k; i = m_parent[i])
00330         {
00331           /* find parent of i if not yet determined */
00332           if (m_parent[i] == -1)
00333             m_parent[i] = k;
00334           m_nonZerosPerCol[i]++;        /* L (k,i) is nonzero */
00335           tags[i] = k;                  /* mark i as visited */
00336         }
00337       }
00338     }
00339   }
00340   
00341   /* construct Lp index array from m_nonZerosPerCol column counts */
00342   Index* Lp = m_matrix._outerIndexPtr();
00343   Lp[0] = 0;
00344   for(Index k = 0; k < size; ++k)
00345     Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1);
00346 
00347   m_matrix.resizeNonZeros(Lp[size]);
00348   
00349   m_isInitialized     = true;
00350   m_info              = Success;
00351   m_analysisIsOk      = true;
00352   m_factorizationIsOk = false;
00353 }
00354 
00355 
00356 template<typename _MatrixType, int _UpLo>
00357 void SimplicialCholesky<_MatrixType,_UpLo>::factorize(const MatrixType& a)
00358 {
00359   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00360   eigen_assert(a.rows()==a.cols());
00361   const Index size = a.rows();
00362   eigen_assert(m_parent.size()==size);
00363   eigen_assert(m_nonZerosPerCol.size()==size);
00364 
00365   const Index* Lp = m_matrix._outerIndexPtr();
00366   Index* Li = m_matrix._innerIndexPtr();
00367   Scalar* Lx = m_matrix._valuePtr();
00368 
00369   ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
00370   ei_declare_aligned_stack_constructed_variable(Index,  pattern, size, 0);
00371   ei_declare_aligned_stack_constructed_variable(Index,  tags, size, 0);
00372 
00373   SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
00374   ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
00375   
00376   bool ok = true;
00377   m_diag.resize(m_LDLt ? size : 0);
00378   
00379   for(Index k = 0; k < size; ++k)
00380   {
00381     // compute nonzero pattern of kth row of L, in topological order
00382     y[k] = 0.0;                     // Y(0:k) is now all zero
00383     Index top = size;               // stack for pattern is empty
00384     tags[k] = k;                    // mark node k as visited
00385     m_nonZerosPerCol[k] = 0;        // count of nonzeros in column k of L
00386     for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
00387     {
00388       Index i = it.index();
00389       if(i <= k)
00390       {
00391         y[i] += internal::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
00392         Index len;
00393         for(len = 0; tags[i] != k; i = m_parent[i])
00394         {
00395           pattern[len++] = i;     /* L(k,i) is nonzero */
00396           tags[i] = k;            /* mark i as visited */
00397         }
00398         while(len > 0)
00399           pattern[--top] = pattern[--len];
00400       }
00401     }
00402 
00403     /* compute numerical values kth row of L (a sparse triangular solve) */
00404     Scalar d = y[k];                  // get D(k,k) and clear Y(k)
00405     y[k] = 0.0;
00406     for(; top < size; ++top)
00407     {
00408       Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
00409       Scalar yi = y[i];             /* get and clear Y(i) */
00410       y[i] = 0.0;
00411       
00412       /* the nonzero entry L(k,i) */
00413       Scalar l_ki;
00414       if(m_LDLt)
00415         l_ki = yi / m_diag[i];       
00416       else
00417         yi = l_ki = yi / Lx[Lp[i]];
00418       
00419       Index p2 = Lp[i] + m_nonZerosPerCol[i];
00420       Index p;
00421       for(p = Lp[i] + (m_LDLt ? 0 : 1); p < p2; ++p)
00422         y[Li[p]] -= internal::conj(Lx[p]) * yi;
00423       d -= l_ki * internal::conj(yi);
00424       Li[p] = k;                          /* store L(k,i) in column form of L */
00425       Lx[p] = l_ki;
00426       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
00427     }
00428     if(m_LDLt)
00429       m_diag[k] = d;
00430     else
00431     {
00432       Index p = Lp[k]+m_nonZerosPerCol[k]++;
00433       Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
00434       Lx[p] = internal::sqrt(d) ;
00435     }
00436     if(d == Scalar(0))
00437     {
00438       ok = false;                         /* failure, D(k,k) is zero */
00439       break;
00440     }
00441   }
00442 
00443   m_info = ok ? Success : NumericalIssue;
00444   m_factorizationIsOk = true;
00445 }
00446 
00447 namespace internal {
00448   
00449 template<typename _MatrixType, int _UpLo, typename Rhs>
00450 struct solve_retval<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
00451   : solve_retval_base<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
00452 {
00453   typedef SimplicialCholesky<_MatrixType,_UpLo> Dec;
00454   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00455 
00456   template<typename Dest> void evalTo(Dest& dst) const
00457   {
00458     dec()._solve(rhs(),dst);
00459   }
00460 };
00461 
00462 template<typename _MatrixType, int _UpLo, typename Rhs>
00463 struct sparse_solve_retval<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
00464   : sparse_solve_retval_base<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
00465 {
00466   typedef SimplicialCholesky<_MatrixType,_UpLo> Dec;
00467   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00468 
00469   template<typename Dest> void evalTo(Dest& dst) const
00470   {
00471     dec()._solve(rhs(),dst);
00472   }
00473 };
00474 
00475 }
00476 
00477 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:22