1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2010 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_SPARSEDENSEPRODUCT_H
11#define EIGEN_SPARSEDENSEPRODUCT_H
12
13namespace Eigen {
14
15template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductReturnType
16{
17  typedef SparseTimeDenseProduct<Lhs,Rhs> Type;
18};
19
20template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1>
21{
22  typedef typename internal::conditional<
23    Lhs::IsRowMajor,
24    SparseDenseOuterProduct<Rhs,Lhs,true>,
25    SparseDenseOuterProduct<Lhs,Rhs,false> >::type Type;
26};
27
28template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
29{
30  typedef DenseTimeSparseProduct<Lhs,Rhs> Type;
31};
32
33template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1>
34{
35  typedef typename internal::conditional<
36    Rhs::IsRowMajor,
37    SparseDenseOuterProduct<Rhs,Lhs,true>,
38    SparseDenseOuterProduct<Lhs,Rhs,false> >::type Type;
39};
40
41namespace internal {
42
43template<typename Lhs, typename Rhs, bool Tr>
44struct traits<SparseDenseOuterProduct<Lhs,Rhs,Tr> >
45{
46  typedef Sparse StorageKind;
47  typedef typename scalar_product_traits<typename traits<Lhs>::Scalar,
48                                         typename traits<Rhs>::Scalar>::ReturnType Scalar;
49  typedef typename Lhs::Index Index;
50  typedef typename Lhs::Nested LhsNested;
51  typedef typename Rhs::Nested RhsNested;
52  typedef typename remove_all<LhsNested>::type _LhsNested;
53  typedef typename remove_all<RhsNested>::type _RhsNested;
54
55  enum {
56    LhsCoeffReadCost = traits<_LhsNested>::CoeffReadCost,
57    RhsCoeffReadCost = traits<_RhsNested>::CoeffReadCost,
58
59    RowsAtCompileTime    = Tr ? int(traits<Rhs>::RowsAtCompileTime)     : int(traits<Lhs>::RowsAtCompileTime),
60    ColsAtCompileTime    = Tr ? int(traits<Lhs>::ColsAtCompileTime)     : int(traits<Rhs>::ColsAtCompileTime),
61    MaxRowsAtCompileTime = Tr ? int(traits<Rhs>::MaxRowsAtCompileTime)  : int(traits<Lhs>::MaxRowsAtCompileTime),
62    MaxColsAtCompileTime = Tr ? int(traits<Lhs>::MaxColsAtCompileTime)  : int(traits<Rhs>::MaxColsAtCompileTime),
63
64    Flags = Tr ? RowMajorBit : 0,
65
66    CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + NumTraits<Scalar>::MulCost
67  };
68};
69
70} // end namespace internal
71
72template<typename Lhs, typename Rhs, bool Tr>
73class SparseDenseOuterProduct
74 : public SparseMatrixBase<SparseDenseOuterProduct<Lhs,Rhs,Tr> >
75{
76  public:
77
78    typedef SparseMatrixBase<SparseDenseOuterProduct> Base;
79    EIGEN_DENSE_PUBLIC_INTERFACE(SparseDenseOuterProduct)
80    typedef internal::traits<SparseDenseOuterProduct> Traits;
81
82  private:
83
84    typedef typename Traits::LhsNested LhsNested;
85    typedef typename Traits::RhsNested RhsNested;
86    typedef typename Traits::_LhsNested _LhsNested;
87    typedef typename Traits::_RhsNested _RhsNested;
88
89  public:
90
91    class InnerIterator;
92
93    EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Lhs& lhs, const Rhs& rhs)
94      : m_lhs(lhs), m_rhs(rhs)
95    {
96      EIGEN_STATIC_ASSERT(!Tr,YOU_MADE_A_PROGRAMMING_MISTAKE);
97    }
98
99    EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Rhs& rhs, const Lhs& lhs)
100      : m_lhs(lhs), m_rhs(rhs)
101    {
102      EIGEN_STATIC_ASSERT(Tr,YOU_MADE_A_PROGRAMMING_MISTAKE);
103    }
104
105    EIGEN_STRONG_INLINE Index rows() const { return Tr ? m_rhs.rows() : m_lhs.rows(); }
106    EIGEN_STRONG_INLINE Index cols() const { return Tr ? m_lhs.cols() : m_rhs.cols(); }
107
108    EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
109    EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
110
111  protected:
112    LhsNested m_lhs;
113    RhsNested m_rhs;
114};
115
116template<typename Lhs, typename Rhs, bool Transpose>
117class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNested::InnerIterator
118{
119    typedef typename _LhsNested::InnerIterator Base;
120    typedef typename SparseDenseOuterProduct::Index Index;
121  public:
122    EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer)
123      : Base(prod.lhs(), 0), m_outer(outer), m_factor(get(prod.rhs(), outer, typename internal::traits<Rhs>::StorageKind() ))
124    { }
125
126    inline Index outer() const { return m_outer; }
127    inline Index row() const { return Transpose ? m_outer : Base::index(); }
128    inline Index col() const { return Transpose ? Base::index() : m_outer; }
129
130    inline Scalar value() const { return Base::value() * m_factor; }
131
132  protected:
133    static Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense())
134    {
135      return rhs.coeff(outer);
136    }
137
138    static Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse())
139    {
140      typename Traits::_RhsNested::InnerIterator it(rhs, outer);
141      if (it && it.index()==0)
142        return it.value();
143
144      return Scalar(0);
145    }
146
147    Index m_outer;
148    Scalar m_factor;
149};
150
151namespace internal {
152template<typename Lhs, typename Rhs>
153struct traits<SparseTimeDenseProduct<Lhs,Rhs> >
154 : traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> >
155{
156  typedef Dense StorageKind;
157  typedef MatrixXpr XprKind;
158};
159
160template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,
161         int LhsStorageOrder = ((SparseLhsType::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor,
162         bool ColPerCol = ((DenseRhsType::Flags&RowMajorBit)==0) || DenseRhsType::ColsAtCompileTime==1>
163struct sparse_time_dense_product_impl;
164
165template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
166struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, RowMajor, true>
167{
168  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
169  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
170  typedef typename internal::remove_all<DenseResType>::type Res;
171  typedef typename Lhs::Index Index;
172  typedef typename Lhs::InnerIterator LhsInnerIterator;
173  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
174  {
175    for(Index c=0; c<rhs.cols(); ++c)
176    {
177      Index n = lhs.outerSize();
178      for(Index j=0; j<n; ++j)
179      {
180        typename Res::Scalar tmp(0);
181        for(LhsInnerIterator it(lhs,j); it ;++it)
182          tmp += it.value() * rhs.coeff(it.index(),c);
183        res.coeffRef(j,c) = alpha * tmp;
184      }
185    }
186  }
187};
188
189template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
190struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, ColMajor, true>
191{
192  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
193  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
194  typedef typename internal::remove_all<DenseResType>::type Res;
195  typedef typename Lhs::InnerIterator LhsInnerIterator;
196  typedef typename Lhs::Index Index;
197  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
198  {
199    for(Index c=0; c<rhs.cols(); ++c)
200    {
201      for(Index j=0; j<lhs.outerSize(); ++j)
202      {
203        typename Res::Scalar rhs_j = alpha * rhs.coeff(j,c);
204        for(LhsInnerIterator it(lhs,j); it ;++it)
205          res.coeffRef(it.index(),c) += it.value() * rhs_j;
206      }
207    }
208  }
209};
210
211template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
212struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, RowMajor, false>
213{
214  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
215  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
216  typedef typename internal::remove_all<DenseResType>::type Res;
217  typedef typename Lhs::InnerIterator LhsInnerIterator;
218  typedef typename Lhs::Index Index;
219  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
220  {
221    for(Index j=0; j<lhs.outerSize(); ++j)
222    {
223      typename Res::RowXpr res_j(res.row(j));
224      for(LhsInnerIterator it(lhs,j); it ;++it)
225        res_j += (alpha*it.value()) * rhs.row(it.index());
226    }
227  }
228};
229
230template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
231struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, ColMajor, false>
232{
233  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
234  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
235  typedef typename internal::remove_all<DenseResType>::type Res;
236  typedef typename Lhs::InnerIterator LhsInnerIterator;
237  typedef typename Lhs::Index Index;
238  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
239  {
240    for(Index j=0; j<lhs.outerSize(); ++j)
241    {
242      typename Rhs::ConstRowXpr rhs_j(rhs.row(j));
243      for(LhsInnerIterator it(lhs,j); it ;++it)
244        res.row(it.index()) += (alpha*it.value()) * rhs_j;
245    }
246  }
247};
248
249template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,typename AlphaType>
250inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
251{
252  sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType>::run(lhs, rhs, res, alpha);
253}
254
255} // end namespace internal
256
257template<typename Lhs, typename Rhs>
258class SparseTimeDenseProduct
259  : public ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs>
260{
261  public:
262    EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct)
263
264    SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
265    {}
266
267    template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const
268    {
269      internal::sparse_time_dense_product(m_lhs, m_rhs, dest, alpha);
270    }
271
272  private:
273    SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&);
274};
275
276
277// dense = dense * sparse
278namespace internal {
279template<typename Lhs, typename Rhs>
280struct traits<DenseTimeSparseProduct<Lhs,Rhs> >
281 : traits<ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> >
282{
283  typedef Dense StorageKind;
284};
285} // end namespace internal
286
287template<typename Lhs, typename Rhs>
288class DenseTimeSparseProduct
289  : public ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs>
290{
291  public:
292    EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct)
293
294    DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
295    {}
296
297    template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const
298    {
299      Transpose<const _LhsNested> lhs_t(m_lhs);
300      Transpose<const _RhsNested> rhs_t(m_rhs);
301      Transpose<Dest> dest_t(dest);
302      internal::sparse_time_dense_product(rhs_t, lhs_t, dest_t, alpha);
303    }
304
305  private:
306    DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&);
307};
308
309// sparse * dense
310template<typename Derived>
311template<typename OtherDerived>
312inline const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
313SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
314{
315  return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
316}
317
318} // end namespace Eigen
319
320#endif // EIGEN_SPARSEDENSEPRODUCT_H
321