SparseSolve.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 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 #ifndef EIGEN_SPARSE_SOLVE_H
00011 #define EIGEN_SPARSE_SOLVE_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base;
00018 template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval;
00019   
00020 template<typename DecompositionType, typename Rhs>
00021 struct traits<sparse_solve_retval_base<DecompositionType, Rhs> >
00022 {
00023   typedef typename DecompositionType::MatrixType MatrixType;
00024   typedef SparseMatrix<typename Rhs::Scalar, Rhs::Options, typename Rhs::Index> ReturnType;
00025 };
00026 
00027 template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base
00028  : public ReturnByValue<sparse_solve_retval_base<_DecompositionType, Rhs> >
00029 {
00030   typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
00031   typedef _DecompositionType DecompositionType;
00032   typedef ReturnByValue<sparse_solve_retval_base> Base;
00033   typedef typename Base::Index Index;
00034 
00035   sparse_solve_retval_base(const DecompositionType& dec, const Rhs& rhs)
00036     : m_dec(dec), m_rhs(rhs)
00037   {}
00038 
00039   inline Index rows() const { return m_dec.cols(); }
00040   inline Index cols() const { return m_rhs.cols(); }
00041   inline const DecompositionType& dec() const { return m_dec; }
00042   inline const RhsNestedCleaned& rhs() const { return m_rhs; }
00043 
00044   template<typename Dest> inline void evalTo(Dest& dst) const
00045   {
00046     static_cast<const sparse_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
00047   }
00048 
00049   protected:
00050     template<typename DestScalar, int DestOptions, typename DestIndex>
00051     inline void defaultEvalTo(SparseMatrix<DestScalar,DestOptions,DestIndex>& dst) const
00052     {
00053       // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
00054       static const int NbColsAtOnce = 4;
00055       int rhsCols = m_rhs.cols();
00056       int size = m_rhs.rows();
00057       Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
00058       Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,rhsCols);
00059       for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00060       {
00061         int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00062         tmp.leftCols(actualCols) = m_rhs.middleCols(k,actualCols);
00063         tmpX.leftCols(actualCols) = m_dec.solve(tmp.leftCols(actualCols));
00064         dst.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView();
00065       }
00066     }
00067     const DecompositionType& m_dec;
00068     typename Rhs::Nested m_rhs;
00069 };
00070 
00071 #define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType,Rhs) \
00072   typedef typename DecompositionType::MatrixType MatrixType; \
00073   typedef typename MatrixType::Scalar Scalar; \
00074   typedef typename MatrixType::RealScalar RealScalar; \
00075   typedef typename MatrixType::Index Index; \
00076   typedef Eigen::internal::sparse_solve_retval_base<DecompositionType,Rhs> Base; \
00077   using Base::dec; \
00078   using Base::rhs; \
00079   using Base::rows; \
00080   using Base::cols; \
00081   sparse_solve_retval(const DecompositionType& dec, const Rhs& rhs) \
00082     : Base(dec, rhs) {}
00083 
00084 
00085 
00086 template<typename DecompositionType, typename Rhs, typename Guess> struct solve_retval_with_guess;
00087 
00088 template<typename DecompositionType, typename Rhs, typename Guess>
00089 struct traits<solve_retval_with_guess<DecompositionType, Rhs, Guess> >
00090 {
00091   typedef typename DecompositionType::MatrixType MatrixType;
00092   typedef Matrix<typename Rhs::Scalar,
00093                  MatrixType::ColsAtCompileTime,
00094                  Rhs::ColsAtCompileTime,
00095                  Rhs::PlainObject::Options,
00096                  MatrixType::MaxColsAtCompileTime,
00097                  Rhs::MaxColsAtCompileTime> ReturnType;
00098 };
00099 
00100 template<typename DecompositionType, typename Rhs, typename Guess> struct solve_retval_with_guess
00101  : public ReturnByValue<solve_retval_with_guess<DecompositionType, Rhs, Guess> >
00102 {
00103   typedef typename DecompositionType::Index Index;
00104 
00105   solve_retval_with_guess(const DecompositionType& dec, const Rhs& rhs, const Guess& guess)
00106     : m_dec(dec), m_rhs(rhs), m_guess(guess)
00107   {}
00108 
00109   inline Index rows() const { return m_dec.cols(); }
00110   inline Index cols() const { return m_rhs.cols(); }
00111 
00112   template<typename Dest> inline void evalTo(Dest& dst) const
00113   {
00114     dst = m_guess;
00115     m_dec._solveWithGuess(m_rhs,dst);
00116   }
00117 
00118   protected:
00119     const DecompositionType& m_dec;
00120     const typename Rhs::Nested m_rhs;
00121     const typename Guess::Nested m_guess;
00122 };
00123 
00124 } // namepsace internal
00125 
00126 } // end namespace Eigen
00127 
00128 #endif // EIGEN_SPARSE_SOLVE_H