SparseLDLTLegacy.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 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_SPARSELDLT_LEGACY_H
00064 #define EIGEN_SPARSELDLT_LEGACY_H
00065 
00078 template<typename _MatrixType, typename Backend = DefaultBackend>
00079 class SparseLDLT
00080 {
00081   protected:
00082     typedef typename _MatrixType::Scalar Scalar;
00083     typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
00084     
00085     typedef Matrix<Scalar,_MatrixType::ColsAtCompileTime,1> VectorType;
00086 
00087     enum {
00088       SupernodalFactorIsDirty      = 0x10000,
00089       MatrixLIsDirty               = 0x20000
00090     };
00091 
00092   public:
00093     typedef SparseMatrix<Scalar> CholMatrixType;
00094     typedef _MatrixType MatrixType;
00095     typedef typename MatrixType::Index Index;
00096 
00097 
00099     SparseLDLT(int flags = 0)
00100       : m_flags(flags), m_status(0)
00101     {
00102       eigen_assert((MatrixType::Flags&RowMajorBit)==0);
00103       m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
00104     }
00105 
00108     SparseLDLT(const MatrixType& matrix, int flags = 0)
00109       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
00110     {
00111       eigen_assert((MatrixType::Flags&RowMajorBit)==0);
00112       m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
00113       compute(matrix);
00114     }
00115 
00129     void setPrecision(RealScalar v) { m_precision = v; }
00130 
00134     RealScalar precision() const { return m_precision; }
00135 
00146     void settags(int f) { m_flags = f; }
00148     int flags() const { return m_flags; }
00149 
00151     void compute(const MatrixType& matrix);
00152 
00154     void _symbolic(const MatrixType& matrix);
00157     bool _numeric(const MatrixType& matrix);
00158 
00160     inline const CholMatrixType& matrixL(void) const { return m_matrix; }
00161 
00163     inline VectorType vectorD(void) const { return m_diag; }
00164 
00165     template<typename Derived>
00166     bool solveInPlace(MatrixBase<Derived> &b) const;
00167 
00168     template<typename Rhs>
00169     inline const internal::solve_retval<SparseLDLT<MatrixType>, Rhs>
00170     solve(const MatrixBase<Rhs>& b) const
00171     {
00172       eigen_assert(true && "SparseLDLT is not initialized.");
00173       return internal::solve_retval<SparseLDLT<MatrixType>, Rhs>(*this, b.derived());
00174     }
00175 
00176     inline Index cols() const { return m_matrix.cols(); }
00177     inline Index rows() const { return m_matrix.rows(); }
00178 
00179     inline const VectorType& diag() const { return m_diag; }
00180 
00182     inline bool succeeded(void) const { return m_succeeded; }
00183 
00184   protected:
00185     CholMatrixType m_matrix;
00186     VectorType m_diag;
00187     VectorXi m_parent; // elimination tree
00188     VectorXi m_nonZerosPerCol;
00189 //     VectorXi m_w; // workspace
00190     PermutationMatrix<Dynamic> m_P;
00191     PermutationMatrix<Dynamic> m_Pinv;
00192     RealScalar m_precision;
00193     int m_flags;
00194     mutable int m_status;
00195     bool m_succeeded;
00196 };
00197 
00198 namespace internal {
00199 
00200 template<typename _MatrixType, typename Rhs>
00201 struct solve_retval<SparseLDLT<_MatrixType>, Rhs>
00202   : solve_retval_base<SparseLDLT<_MatrixType>, Rhs>
00203 {
00204   typedef SparseLDLT<_MatrixType> SpLDLTDecType;
00205   EIGEN_MAKE_SOLVE_HELPERS(SpLDLTDecType,Rhs)
00206 
00207   template<typename Dest> void evalTo(Dest& dst) const
00208   {
00209     //Index size = dec().matrixL().rows();
00210     eigen_assert(dec().matrixL().rows()==rhs().rows());
00211 
00212     Rhs b(rhs().rows(), rhs().cols());
00213     b = rhs();
00214 
00215     if (dec().matrixL().nonZeros()>0) // otherwise L==I
00216       dec().matrixL().template triangularView<UnitLower>().solveInPlace(b);
00217 
00218     b = b.cwiseQuotient(dec().diag());
00219     if (dec().matrixL().nonZeros()>0) // otherwise L==I
00220       dec().matrixL().adjoint().template triangularView<UnitUpper>().solveInPlace(b);
00221     
00222     dst = b;
00223 
00224   }
00225     
00226 };
00227 
00228 } // end namespace internal
00229 
00233 template<typename _MatrixType, typename Backend>
00234 void SparseLDLT<_MatrixType,Backend>::compute(const _MatrixType& a)
00235 {
00236   _symbolic(a);
00237   m_succeeded = _numeric(a);
00238 }
00239 
00240 template<typename _MatrixType, typename Backend>
00241 void SparseLDLT<_MatrixType,Backend>::_symbolic(const _MatrixType& a)
00242 {
00243   assert(a.rows()==a.cols());
00244   const Index size = a.rows();
00245   m_matrix.resize(size, size);
00246   m_parent.resize(size);
00247   m_nonZerosPerCol.resize(size);
00248 
00249   ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
00250 
00251   const Index* Ap = a._outerIndexPtr();
00252   const Index* Ai = a._innerIndexPtr();
00253   Index* Lp = m_matrix._outerIndexPtr();
00254   
00255   const Index* P = 0;
00256   Index* Pinv = 0;
00257 
00258   if(P)
00259   {
00260     m_P.indices()     = VectorXi::Map(P,size);
00261     m_Pinv = m_P.inverse();
00262     Pinv = m_Pinv.indices().data();
00263   }
00264   else
00265   {
00266     m_P.resize(0);
00267     m_Pinv.resize(0);
00268   }
00269 
00270   for (Index k = 0; k < size; ++k)
00271   {
00272     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
00273     m_parent[k] = -1;             /* parent of k is not yet known */
00274     tags[k] = k;                  /* mark node k as visited */
00275     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
00276     Index kk = P ? P[k] : k;      /* kth original, or permuted, column */
00277     Index p2 = Ap[kk+1];
00278     for (Index p = Ap[kk]; p < p2; ++p)
00279     {
00280       /* A (i,k) is nonzero (original or permuted A) */
00281       Index i = Pinv ? Pinv[Ai[p]] : Ai[p];
00282       if (i < k)
00283       {
00284         /* follow path from i to root of etree, stop at flagged node */
00285         for (; tags[i] != k; i = m_parent[i])
00286         {
00287           /* find parent of i if not yet determined */
00288           if (m_parent[i] == -1)
00289             m_parent[i] = k;
00290           ++m_nonZerosPerCol[i];        /* L (k,i) is nonzero */
00291           tags[i] = k;                  /* mark i as visited */
00292         }
00293       }
00294     }
00295   }
00296   /* construct Lp index array from m_nonZerosPerCol column counts */
00297   Lp[0] = 0;
00298   for (Index k = 0; k < size; ++k)
00299     Lp[k+1] = Lp[k] + m_nonZerosPerCol[k];
00300 
00301   m_matrix.resizeNonZeros(Lp[size]);
00302 }
00303 
00304 template<typename _MatrixType, typename Backend>
00305 bool SparseLDLT<_MatrixType,Backend>::_numeric(const _MatrixType& a)
00306 {
00307   assert(a.rows()==a.cols());
00308   const Index size = a.rows();
00309   assert(m_parent.size()==size);
00310   assert(m_nonZerosPerCol.size()==size);
00311 
00312   const Index* Ap = a._outerIndexPtr();
00313   const Index* Ai = a._innerIndexPtr();
00314   const Scalar* Ax = a._valuePtr();
00315   const Index* Lp = m_matrix._outerIndexPtr();
00316   Index* Li = m_matrix._innerIndexPtr();
00317   Scalar* Lx = m_matrix._valuePtr();
00318   m_diag.resize(size);
00319 
00320   ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
00321   ei_declare_aligned_stack_constructed_variable(Index,  pattern, size, 0);
00322   ei_declare_aligned_stack_constructed_variable(Index,  tags, size, 0);
00323   
00324   Index* P = 0;
00325   Index* Pinv = 0;
00326   
00327   if(m_P.size()==size)
00328   {
00329     P = m_P.indices().data();
00330     Pinv = m_Pinv.indices().data();
00331   }
00332   
00333   bool ok = true;
00334 
00335   for (Index k = 0; k < size; ++k)
00336   {
00337     /* compute nonzero pattern of kth row of L, in topological order */
00338     y[k] = 0.0;                     /* Y(0:k) is now all zero */
00339     Index top = size;               /* stack for pattern is empty */
00340     tags[k] = k;                    /* mark node k as visited */
00341     m_nonZerosPerCol[k] = 0;        /* count of nonzeros in column k of L */
00342     Index kk = (P) ? (P[k]) : (k);  /* kth original, or permuted, column */
00343     Index p2 = Ap[kk+1];
00344     for (Index p = Ap[kk]; p < p2; ++p)
00345     {
00346       Index i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */
00347       if (i <= k)
00348       {
00349         y[i] += internal::conj(Ax[p]);            /* scatter A(i,k) into Y (sum duplicates) */
00350         Index len;
00351         for (len = 0; tags[i] != k; i = m_parent[i])
00352         {
00353           pattern[len++] = i;     /* L(k,i) is nonzero */
00354           tags[i] = k;            /* mark i as visited */
00355         }
00356         while (len > 0)
00357           pattern[--top] = pattern[--len];
00358       }
00359     }
00360 
00361     /* compute numerical values kth row of L (a sparse triangular solve) */
00362     m_diag[k] = y[k];                       /* get D(k,k) and clear Y(k) */
00363     y[k] = 0.0;
00364     for (; top < size; ++top)
00365     {
00366       Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
00367       Scalar yi = (y[i]);             /* get and clear Y(i) */
00368       y[i] = 0.0;
00369       Index p2 = Lp[i] + m_nonZerosPerCol[i];
00370       Index p;
00371       for (p = Lp[i]; p < p2; ++p)
00372         y[Li[p]] -= internal::conj(Lx[p]) * (yi);
00373       Scalar l_ki = yi / m_diag[i];       /* the nonzero entry L(k,i) */
00374       m_diag[k] -= l_ki * internal::conj(yi);
00375       Li[p] = k;                          /* store L(k,i) in column form of L */
00376       Lx[p] = (l_ki);
00377       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
00378     }
00379     if (m_diag[k] == 0.0)
00380     {
00381       ok = false;                         /* failure, D(k,k) is zero */
00382       break;
00383     }
00384   }
00385 
00386   return ok;  /* success, diagonal of D is all nonzero */
00387 }
00388 
00390 template<typename _MatrixType, typename Backend>
00391 template<typename Derived>
00392 bool SparseLDLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
00393 {
00394   //Index size = m_matrix.rows();
00395   eigen_assert(m_matrix.rows()==b.rows());
00396   if (!m_succeeded)
00397     return false;
00398   
00399   if(m_P.size()>0)
00400     b = m_Pinv * b;
00401 
00402   if (m_matrix.nonZeros()>0) // otherwise L==I
00403     m_matrix.template triangularView<UnitLower>().solveInPlace(b);
00404   b = b.cwiseQuotient(m_diag);
00405   if (m_matrix.nonZeros()>0) // otherwise L==I
00406     m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(b);
00407 
00408   if(m_P.size()>0)
00409     b = m_P * b;
00410   
00411   return true;
00412 }
00413 
00414 #endif // EIGEN_SPARSELDLT_LEGACY_H


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