src/SparseCore/SparseSolverBase.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SPARSESOLVERBASE_H
0011 #define EIGEN_SPARSESOLVERBASE_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017   /** \internal
0018   * Helper functions to solve with a sparse right-hand-side and result.
0019   * The rhs is decomposed into small vertical panels which are solved through dense temporaries.
0020   */
0021 template<typename Decomposition, typename Rhs, typename Dest>
0022 typename enable_if<Rhs::ColsAtCompileTime!=1 && Dest::ColsAtCompileTime!=1>::type
0023 solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest)
0024 {
0025   EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
0026   typedef typename Dest::Scalar DestScalar;
0027   // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
0028   static const Index NbColsAtOnce = 4;
0029   Index rhsCols = rhs.cols();
0030   Index size = rhs.rows();
0031   // the temporary matrices do not need more columns than NbColsAtOnce:
0032   Index tmpCols = (std::min)(rhsCols, NbColsAtOnce);
0033   Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,tmpCols);
0034   Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,tmpCols);
0035   for(Index k=0; k<rhsCols; k+=NbColsAtOnce)
0036   {
0037     Index actualCols = std::min<Index>(rhsCols-k, NbColsAtOnce);
0038     tmp.leftCols(actualCols) = rhs.middleCols(k,actualCols);
0039     tmpX.leftCols(actualCols) = dec.solve(tmp.leftCols(actualCols));
0040     dest.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView();
0041   }
0042 }
0043
0044 // Overload for vector as rhs
0045 template<typename Decomposition, typename Rhs, typename Dest>
0046 typename enable_if<Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1>::type
0047 solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest)
0048 {
0049   typedef typename Dest::Scalar DestScalar;
0050   Index size = rhs.rows();
0051   Eigen::Matrix<DestScalar,Dynamic,1> rhs_dense(rhs);
0052   Eigen::Matrix<DestScalar,Dynamic,1> dest_dense(size);
0053   dest_dense = dec.solve(rhs_dense);
0054   dest = dest_dense.sparseView();
0055 }
0056
0057 } // end namespace internal
0058
0059 /** \class SparseSolverBase
0060   * \ingroup SparseCore_Module
0061   * \brief A base class for sparse solvers
0062   *
0063   * \tparam Derived the actual type of the solver.
0064   *
0065   */
0066 template<typename Derived>
0067 class SparseSolverBase : internal::noncopyable
0068 {
0069   public:
0070
0071     /** Default constructor */
0072     SparseSolverBase()
0073       : m_isInitialized(false)
0074     {}
0075
0076     ~SparseSolverBase()
0077     {}
0078
0079     Derived& derived() { return *static_cast<Derived*>(this); }
0080     const Derived& derived() const { return *static_cast<const Derived*>(this); }
0081
0082     /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A.
0083       *
0084       * \sa compute()
0085       */
0086     template<typename Rhs>
0087     inline const Solve<Derived, Rhs>
0088     solve(const MatrixBase<Rhs>& b) const
0089     {
0090       eigen_assert(m_isInitialized && "Solver is not initialized.");
0091       eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
0092       return Solve<Derived, Rhs>(derived(), b.derived());
0093     }
0094
0095     /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A.
0096       *
0097       * \sa compute()
0098       */
0099     template<typename Rhs>
0100     inline const Solve<Derived, Rhs>
0101     solve(const SparseMatrixBase<Rhs>& b) const
0102     {
0103       eigen_assert(m_isInitialized && "Solver is not initialized.");
0104       eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
0105       return Solve<Derived, Rhs>(derived(), b.derived());
0106     }
0107
0108     #ifndef EIGEN_PARSED_BY_DOXYGEN
0109     /** \internal default implementation of solving with a sparse rhs */
0110     template<typename Rhs,typename Dest>
0111     void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
0112     {
0113       internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived());
0114     }
0115     #endif // EIGEN_PARSED_BY_DOXYGEN
0116
0117   protected:
0118
0119     mutable bool m_isInitialized;
0120 };
0121
0122 } // end namespace Eigen
0123
0124 #endif // EIGEN_SPARSESOLVERBASE_H