![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2014 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_SPARSESOLVERBASE_H 00011 #define EIGEN_SPARSESOLVERBASE_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 00021 template<typename Decomposition, typename Rhs, typename Dest> 00022 typename enable_if<Rhs::ColsAtCompileTime!=1 && Dest::ColsAtCompileTime!=1>::type 00023 solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) 00024 { 00025 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 00026 typedef typename Dest::Scalar DestScalar; 00027 // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. 00028 static const Index NbColsAtOnce = 4; 00029 Index rhsCols = rhs.cols(); 00030 Index size = rhs.rows(); 00031 // the temporary matrices do not need more columns than NbColsAtOnce: 00032 Index tmpCols = (std::min)(rhsCols, NbColsAtOnce); 00033 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,tmpCols); 00034 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,tmpCols); 00035 for(Index k=0; k<rhsCols; k+=NbColsAtOnce) 00036 { 00037 Index actualCols = std::min<Index>(rhsCols-k, NbColsAtOnce); 00038 tmp.leftCols(actualCols) = rhs.middleCols(k,actualCols); 00039 tmpX.leftCols(actualCols) = dec.solve(tmp.leftCols(actualCols)); 00040 dest.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView(); 00041 } 00042 } 00043 00044 // Overload for vector as rhs 00045 template<typename Decomposition, typename Rhs, typename Dest> 00046 typename enable_if<Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1>::type 00047 solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) 00048 { 00049 typedef typename Dest::Scalar DestScalar; 00050 Index size = rhs.rows(); 00051 Eigen::Matrix<DestScalar,Dynamic,1> rhs_dense(rhs); 00052 Eigen::Matrix<DestScalar,Dynamic,1> dest_dense(size); 00053 dest_dense = dec.solve(rhs_dense); 00054 dest = dest_dense.sparseView(); 00055 } 00056 00057 } // end namespace internal 00058 00066 template<typename Derived> 00067 class SparseSolverBase : internal::noncopyable 00068 { 00069 public: 00070 00072 SparseSolverBase() 00073 : m_isInitialized(false) 00074 {} 00075 00076 ~SparseSolverBase() 00077 {} 00078 00079 Derived& derived() { return *static_cast<Derived*>(this); } 00080 const Derived& derived() const { return *static_cast<const Derived*>(this); } 00081 00086 template<typename Rhs> 00087 inline const Solve<Derived, Rhs> 00088 solve(const MatrixBase<Rhs>& b) const 00089 { 00090 eigen_assert(m_isInitialized && "Solver is not initialized."); 00091 eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); 00092 return Solve<Derived, Rhs>(derived(), b.derived()); 00093 } 00094 00099 template<typename Rhs> 00100 inline const Solve<Derived, Rhs> 00101 solve(const SparseMatrixBase<Rhs>& b) const 00102 { 00103 eigen_assert(m_isInitialized && "Solver is not initialized."); 00104 eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); 00105 return Solve<Derived, Rhs>(derived(), b.derived()); 00106 } 00107 00108 #ifndef EIGEN_PARSED_BY_DOXYGEN 00109 00110 template<typename Rhs,typename Dest> 00111 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 00112 { 00113 internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived()); 00114 } 00115 #endif // EIGEN_PARSED_BY_DOXYGEN 00116 00117 protected: 00118 00119 mutable bool m_isInitialized; 00120 }; 00121 00122 } // end namespace Eigen 00123 00124 #endif // EIGEN_SPARSESOLVERBASE_H