Michael Ernst Peter / Eigen

src/misc/SparseSolve.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 5 //
ykuroda 0:13a5d365ba16 6 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 7 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 9
ykuroda 0:13a5d365ba16 10 #ifndef EIGEN_SPARSE_SOLVE_H
ykuroda 0:13a5d365ba16 11 #define EIGEN_SPARSE_SOLVE_H
ykuroda 0:13a5d365ba16 12
ykuroda 0:13a5d365ba16 13 namespace Eigen {
ykuroda 0:13a5d365ba16 14
ykuroda 0:13a5d365ba16 15 namespace internal {
ykuroda 0:13a5d365ba16 16
ykuroda 0:13a5d365ba16 17 template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base;
ykuroda 0:13a5d365ba16 18 template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval;
ykuroda 0:13a5d365ba16 19
ykuroda 0:13a5d365ba16 20 template<typename DecompositionType, typename Rhs>
ykuroda 0:13a5d365ba16 21 struct traits<sparse_solve_retval_base<DecompositionType, Rhs> >
ykuroda 0:13a5d365ba16 22 {
ykuroda 0:13a5d365ba16 23 typedef typename DecompositionType::MatrixType MatrixType;
ykuroda 0:13a5d365ba16 24 typedef SparseMatrix<typename Rhs::Scalar, Rhs::Options, typename Rhs::Index> ReturnType;
ykuroda 0:13a5d365ba16 25 };
ykuroda 0:13a5d365ba16 26
ykuroda 0:13a5d365ba16 27 template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base
ykuroda 0:13a5d365ba16 28 : public ReturnByValue<sparse_solve_retval_base<_DecompositionType, Rhs> >
ykuroda 0:13a5d365ba16 29 {
ykuroda 0:13a5d365ba16 30 typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
ykuroda 0:13a5d365ba16 31 typedef _DecompositionType DecompositionType;
ykuroda 0:13a5d365ba16 32 typedef ReturnByValue<sparse_solve_retval_base> Base;
ykuroda 0:13a5d365ba16 33 typedef typename Base::Index Index;
ykuroda 0:13a5d365ba16 34
ykuroda 0:13a5d365ba16 35 sparse_solve_retval_base(const DecompositionType& dec, const Rhs& rhs)
ykuroda 0:13a5d365ba16 36 : m_dec(dec), m_rhs(rhs)
ykuroda 0:13a5d365ba16 37 {}
ykuroda 0:13a5d365ba16 38
ykuroda 0:13a5d365ba16 39 inline Index rows() const { return m_dec.cols(); }
ykuroda 0:13a5d365ba16 40 inline Index cols() const { return m_rhs.cols(); }
ykuroda 0:13a5d365ba16 41 inline const DecompositionType& dec() const { return m_dec; }
ykuroda 0:13a5d365ba16 42 inline const RhsNestedCleaned& rhs() const { return m_rhs; }
ykuroda 0:13a5d365ba16 43
ykuroda 0:13a5d365ba16 44 template<typename Dest> inline void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 45 {
ykuroda 0:13a5d365ba16 46 static_cast<const sparse_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
ykuroda 0:13a5d365ba16 47 }
ykuroda 0:13a5d365ba16 48
ykuroda 0:13a5d365ba16 49 protected:
ykuroda 0:13a5d365ba16 50 template<typename DestScalar, int DestOptions, typename DestIndex>
ykuroda 0:13a5d365ba16 51 inline void defaultEvalTo(SparseMatrix<DestScalar,DestOptions,DestIndex>& dst) const
ykuroda 0:13a5d365ba16 52 {
ykuroda 0:13a5d365ba16 53 // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
ykuroda 0:13a5d365ba16 54 static const int NbColsAtOnce = 4;
ykuroda 0:13a5d365ba16 55 int rhsCols = m_rhs.cols();
ykuroda 0:13a5d365ba16 56 int size = m_rhs.rows();
ykuroda 0:13a5d365ba16 57 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
ykuroda 0:13a5d365ba16 58 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,rhsCols);
ykuroda 0:13a5d365ba16 59 for(int k=0; k<rhsCols; k+=NbColsAtOnce)
ykuroda 0:13a5d365ba16 60 {
ykuroda 0:13a5d365ba16 61 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
ykuroda 0:13a5d365ba16 62 tmp.leftCols(actualCols) = m_rhs.middleCols(k,actualCols);
ykuroda 0:13a5d365ba16 63 tmpX.leftCols(actualCols) = m_dec.solve(tmp.leftCols(actualCols));
ykuroda 0:13a5d365ba16 64 dst.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView();
ykuroda 0:13a5d365ba16 65 }
ykuroda 0:13a5d365ba16 66 }
ykuroda 0:13a5d365ba16 67 const DecompositionType& m_dec;
ykuroda 0:13a5d365ba16 68 typename Rhs::Nested m_rhs;
ykuroda 0:13a5d365ba16 69 };
ykuroda 0:13a5d365ba16 70
ykuroda 0:13a5d365ba16 71 #define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType,Rhs) \
ykuroda 0:13a5d365ba16 72 typedef typename DecompositionType::MatrixType MatrixType; \
ykuroda 0:13a5d365ba16 73 typedef typename MatrixType::Scalar Scalar; \
ykuroda 0:13a5d365ba16 74 typedef typename MatrixType::RealScalar RealScalar; \
ykuroda 0:13a5d365ba16 75 typedef typename MatrixType::Index Index; \
ykuroda 0:13a5d365ba16 76 typedef Eigen::internal::sparse_solve_retval_base<DecompositionType,Rhs> Base; \
ykuroda 0:13a5d365ba16 77 using Base::dec; \
ykuroda 0:13a5d365ba16 78 using Base::rhs; \
ykuroda 0:13a5d365ba16 79 using Base::rows; \
ykuroda 0:13a5d365ba16 80 using Base::cols; \
ykuroda 0:13a5d365ba16 81 sparse_solve_retval(const DecompositionType& dec, const Rhs& rhs) \
ykuroda 0:13a5d365ba16 82 : Base(dec, rhs) {}
ykuroda 0:13a5d365ba16 83
ykuroda 0:13a5d365ba16 84
ykuroda 0:13a5d365ba16 85
ykuroda 0:13a5d365ba16 86 template<typename DecompositionType, typename Rhs, typename Guess> struct solve_retval_with_guess;
ykuroda 0:13a5d365ba16 87
ykuroda 0:13a5d365ba16 88 template<typename DecompositionType, typename Rhs, typename Guess>
ykuroda 0:13a5d365ba16 89 struct traits<solve_retval_with_guess<DecompositionType, Rhs, Guess> >
ykuroda 0:13a5d365ba16 90 {
ykuroda 0:13a5d365ba16 91 typedef typename DecompositionType::MatrixType MatrixType;
ykuroda 0:13a5d365ba16 92 typedef Matrix<typename Rhs::Scalar,
ykuroda 0:13a5d365ba16 93 MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 94 Rhs::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 95 Rhs::PlainObject::Options,
ykuroda 0:13a5d365ba16 96 MatrixType::MaxColsAtCompileTime,
ykuroda 0:13a5d365ba16 97 Rhs::MaxColsAtCompileTime> ReturnType;
ykuroda 0:13a5d365ba16 98 };
ykuroda 0:13a5d365ba16 99
ykuroda 0:13a5d365ba16 100 template<typename DecompositionType, typename Rhs, typename Guess> struct solve_retval_with_guess
ykuroda 0:13a5d365ba16 101 : public ReturnByValue<solve_retval_with_guess<DecompositionType, Rhs, Guess> >
ykuroda 0:13a5d365ba16 102 {
ykuroda 0:13a5d365ba16 103 typedef typename DecompositionType::Index Index;
ykuroda 0:13a5d365ba16 104
ykuroda 0:13a5d365ba16 105 solve_retval_with_guess(const DecompositionType& dec, const Rhs& rhs, const Guess& guess)
ykuroda 0:13a5d365ba16 106 : m_dec(dec), m_rhs(rhs), m_guess(guess)
ykuroda 0:13a5d365ba16 107 {}
ykuroda 0:13a5d365ba16 108
ykuroda 0:13a5d365ba16 109 inline Index rows() const { return m_dec.cols(); }
ykuroda 0:13a5d365ba16 110 inline Index cols() const { return m_rhs.cols(); }
ykuroda 0:13a5d365ba16 111
ykuroda 0:13a5d365ba16 112 template<typename Dest> inline void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 113 {
ykuroda 0:13a5d365ba16 114 dst = m_guess;
ykuroda 0:13a5d365ba16 115 m_dec._solveWithGuess(m_rhs,dst);
ykuroda 0:13a5d365ba16 116 }
ykuroda 0:13a5d365ba16 117
ykuroda 0:13a5d365ba16 118 protected:
ykuroda 0:13a5d365ba16 119 const DecompositionType& m_dec;
ykuroda 0:13a5d365ba16 120 const typename Rhs::Nested m_rhs;
ykuroda 0:13a5d365ba16 121 const typename Guess::Nested m_guess;
ykuroda 0:13a5d365ba16 122 };
ykuroda 0:13a5d365ba16 123
ykuroda 0:13a5d365ba16 124 } // namepsace internal
ykuroda 0:13a5d365ba16 125
ykuroda 0:13a5d365ba16 126 } // end namespace Eigen
ykuroda 0:13a5d365ba16 127
ykuroda 0:13a5d365ba16 128 #endif // EIGEN_SPARSE_SOLVE_H