SparseLLT.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 #ifndef EIGEN_SPARSELLT_H
00026 #define EIGEN_SPARSELLT_H
00027 
00038 template<typename _MatrixType, typename Backend = DefaultBackend>
00039 class SparseLLT
00040 {
00041   protected:
00042     typedef typename _MatrixType::Scalar Scalar;
00043     typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
00044 
00045     enum {
00046       SupernodalFactorIsDirty      = 0x10000,
00047       MatrixLIsDirty               = 0x20000
00048     };
00049 
00050   public:
00051     typedef SparseMatrix<Scalar> CholMatrixType;
00052     typedef _MatrixType MatrixType;
00053     typedef typename MatrixType::Index Index;
00054 
00056     SparseLLT(int flags = 0)
00057       : m_flags(flags), m_status(0)
00058     {
00059       m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
00060     }
00061 
00064     SparseLLT(const MatrixType& matrix, int flags = 0)
00065       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
00066     {
00067       m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
00068       compute(matrix);
00069     }
00070 
00084     void setPrecision(RealScalar v) { m_precision = v; }
00085 
00089     RealScalar precision() const { return m_precision; }
00090 
00101     void setFlags(int f) { m_flags = f; }
00103     int flags() const { return m_flags; }
00104 
00106     void compute(const MatrixType& matrix);
00107 
00109     inline const CholMatrixType& matrixL(void) const { return m_matrix; }
00110 
00111     template<typename Derived>
00112     bool solveInPlace(MatrixBase<Derived> &b) const;
00113 
00114     template<typename Rhs>
00115     inline const internal::solve_retval<SparseLLT<MatrixType>, Rhs>
00116     solve(const MatrixBase<Rhs>& b) const
00117     {
00118       eigen_assert(true && "SparseLLT is not initialized.");
00119       return internal::solve_retval<SparseLLT<MatrixType>, Rhs>(*this, b.derived());
00120     }
00121 
00122     inline Index cols() const { return m_matrix.cols(); }
00123     inline Index rows() const { return m_matrix.rows(); }
00124 
00126     inline bool succeeded(void) const { return m_succeeded; }
00127 
00128   protected:
00129     CholMatrixType m_matrix;
00130     RealScalar m_precision;
00131     int m_flags;
00132     mutable int m_status;
00133     bool m_succeeded;
00134 };
00135 
00136 
00137 namespace internal {
00138 
00139 template<typename _MatrixType, typename Rhs>
00140 struct solve_retval<SparseLLT<_MatrixType>, Rhs>
00141   : solve_retval_base<SparseLLT<_MatrixType>, Rhs>
00142 {
00143   typedef SparseLLT<_MatrixType> SpLLTDecType;
00144   EIGEN_MAKE_SOLVE_HELPERS(SpLLTDecType,Rhs)
00145 
00146   template<typename Dest> void evalTo(Dest& dst) const
00147   {
00148     const Index size = dec().matrixL().rows();
00149     eigen_assert(size==rhs().rows());
00150     
00151     Rhs b(rhs().rows(), rhs().cols());
00152     b = rhs();
00153     
00154     dec().matrixL().template triangularView<Lower>().solveInPlace(b);
00155     dec().matrixL().adjoint().template triangularView<Upper>().solveInPlace(b);
00156     
00157     dst = b;
00158     
00159   }
00160     
00161 };
00162 
00163 } // end namespace internal
00164 
00165 
00169 template<typename _MatrixType, typename Backend>
00170 void SparseLLT<_MatrixType,Backend>::compute(const _MatrixType& a)
00171 {
00172   assert(a.rows()==a.cols());
00173   const Index size = a.rows();
00174   m_matrix.resize(size, size);
00175 
00176   // allocate a temporary vector for accumulations
00177   AmbiVector<Scalar,Index> tempVector(size);
00178   RealScalar density = a.nonZeros()/RealScalar(size*size);
00179 
00180   // TODO estimate the number of non zeros
00181   m_matrix.setZero();
00182   m_matrix.reserve(a.nonZeros()*2);
00183   for (Index j = 0; j < size; ++j)
00184   {
00185     Scalar x = internal::real(a.coeff(j,j));
00186 
00187     // TODO better estimate of the density !
00188     tempVector.init(density>0.001? IsDense : IsSparse);
00189     tempVector.setBounds(j+1,size);
00190     tempVector.setZero();
00191     // init with current matrix a
00192     {
00193       typename _MatrixType::InnerIterator it(a,j);
00194       eigen_assert(it.index()==j &&
00195         "matrix must has non zero diagonal entries and only the lower triangular part must be stored");
00196       ++it; // skip diagonal element
00197       for (; it; ++it)
00198         tempVector.coeffRef(it.index()) = it.value();
00199     }
00200     for (Index k=0; k<j+1; ++k)
00201     {
00202       typename CholMatrixType::InnerIterator it(m_matrix, k);
00203       while (it && it.index()<j)
00204         ++it;
00205       if (it && it.index()==j)
00206       {
00207         Scalar y = it.value();
00208         x -= internal::abs2(y);
00209         ++it; // skip j-th element, and process remaining column coefficients
00210         tempVector.restart();
00211         for (; it; ++it)
00212         {
00213           tempVector.coeffRef(it.index()) -= it.value() * y;
00214         }
00215       }
00216     }
00217     // copy the temporary vector to the respective m_matrix.col()
00218     // while scaling the result by 1/real(x)
00219     RealScalar rx = internal::sqrt(internal::real(x));
00220     m_matrix.insert(j,j) = rx; // FIXME use insertBack
00221     Scalar y = Scalar(1)/rx;
00222     for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector, m_precision*rx); it; ++it)
00223     {
00224       // FIXME use insertBack
00225       m_matrix.insert(it.index(), j) = it.value() * y;
00226     }
00227   }
00228   m_matrix.finalize();
00229 }
00230 
00232 template<typename _MatrixType, typename Backend>
00233 template<typename Derived>
00234 bool SparseLLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
00235 {
00236   const Index size = m_matrix.rows();
00237   eigen_assert(size==b.rows());
00238 
00239   m_matrix.template triangularView<Lower>().solveInPlace(b);
00240   m_matrix.adjoint().template triangularView<Upper>().solveInPlace(b);
00241 
00242   return true;
00243 }
00244 
00245 #endif // EIGEN_SPARSELLT_H


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