1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_SPARSESOLVERBASE_H
11#define EIGEN_SPARSESOLVERBASE_H
12
13namespace Eigen {
14
15namespace internal {
16
17 /** \internal
18 * Helper functions to solve with a sparse right-hand-side and result.
19 * The rhs is decomposed into small vertical panels which are solved through dense temporaries.
20 */
21template<typename Decomposition, typename Rhs, typename Dest>
22typename enable_if<Rhs::ColsAtCompileTime!=1 && Dest::ColsAtCompileTime!=1>::type
23solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest)
24{
25 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
26 typedef typename Dest::Scalar DestScalar;
27 // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
28 static const Index NbColsAtOnce = 4;
29 Index rhsCols = rhs.cols();
30 Index size = rhs.rows();
31 // the temporary matrices do not need more columns than NbColsAtOnce:
32 Index tmpCols = (std::min)(rhsCols, NbColsAtOnce);
33 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,tmpCols);
34 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,tmpCols);
35 for(Index k=0; k<rhsCols; k+=NbColsAtOnce)
36 {
37 Index actualCols = std::min<Index>(rhsCols-k, NbColsAtOnce);
38 tmp.leftCols(actualCols) = rhs.middleCols(k,actualCols);
39 tmpX.leftCols(actualCols) = dec.solve(tmp.leftCols(actualCols));
40 dest.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView();
41 }
42}
43
44// Overload for vector as rhs
45template<typename Decomposition, typename Rhs, typename Dest>
46typename enable_if<Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1>::type
47solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest)
48{
49 typedef typename Dest::Scalar DestScalar;
50 Index size = rhs.rows();
51 Eigen::Matrix<DestScalar,Dynamic,1> rhs_dense(rhs);
52 Eigen::Matrix<DestScalar,Dynamic,1> dest_dense(size);
53 dest_dense = dec.solve(rhs_dense);
54 dest = dest_dense.sparseView();
55}
56
57} // end namespace internal
58
59/** \class SparseSolverBase
60 * \ingroup SparseCore_Module
61 * \brief A base class for sparse solvers
62 *
63 * \tparam Derived the actual type of the solver.
64 *
65 */
66template<typename Derived>
67class SparseSolverBase : internal::noncopyable
68{
69 public:
70
71 /** Default constructor */
72 SparseSolverBase()
73 : m_isInitialized(false)
74 {}
75
76 ~SparseSolverBase()
77 {}
78
79 Derived& derived() { return *static_cast<Derived*>(this); }
80 const Derived& derived() const { return *static_cast<const Derived*>(this); }
81
82 /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A.
83 *
84 * \sa compute()
85 */
86 template<typename Rhs>
87 inline const Solve<Derived, Rhs>
88 solve(const MatrixBase<Rhs>& b) const
89 {
90 eigen_assert(m_isInitialized && "Solver is not initialized.");
91 eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
92 return Solve<Derived, Rhs>(derived(), b.derived());
93 }
94
95 /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A.
96 *
97 * \sa compute()
98 */
99 template<typename Rhs>
100 inline const Solve<Derived, Rhs>
101 solve(const SparseMatrixBase<Rhs>& b) const
102 {
103 eigen_assert(m_isInitialized && "Solver is not initialized.");
104 eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
105 return Solve<Derived, Rhs>(derived(), b.derived());
106 }
107
108 #ifndef EIGEN_PARSED_BY_DOXYGEN
109 /** \internal default implementation of solving with a sparse rhs */
110 template<typename Rhs,typename Dest>
111 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
112 {
113 internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived());
114 }
115 #endif // EIGEN_PARSED_BY_DOXYGEN
116
117 protected:
118
119 mutable bool m_isInitialized;
120};
121
122} // end namespace Eigen
123
124#endif // EIGEN_SPARSESOLVERBASE_H
125