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 
13 namespace Eigen {
14 
15 template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductReturnType
16 {
17   typedef SparseTimeDenseProduct<Lhs,Rhs> Type;
18 };
19 
20 template<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 
28 template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
29 {
30   typedef DenseTimeSparseProduct<Lhs,Rhs> Type;
31 };
32 
33 template<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 
41 namespace internal {
42 
43 template<typename Lhs, typename Rhs, bool Tr>
44 struct 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 
72 template<typename Lhs, typename Rhs, bool Tr>
73 class 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 
116 template<typename Lhs, typename Rhs, bool Transpose>
117 class 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 
151 namespace internal {
152 template<typename Lhs, typename Rhs>
153 struct traits<SparseTimeDenseProduct<Lhs,Rhs> >
154  : traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> >
155 {
156   typedef Dense StorageKind;
157   typedef MatrixXpr XprKind;
158 };
159 
160 template<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>
163 struct sparse_time_dense_product_impl;
164 
165 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
166 struct 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 
189 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
190 struct 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 
211 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
212 struct 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 
230 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
231 struct 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 
249 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,typename AlphaType>
250 inline 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 
257 template<typename Lhs, typename Rhs>
258 class 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
278 namespace internal {
279 template<typename Lhs, typename Rhs>
280 struct traits<DenseTimeSparseProduct<Lhs,Rhs> >
281  : traits<ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> >
282 {
283   typedef Dense StorageKind;
284 };
285 } // end namespace internal
286 
287 template<typename Lhs, typename Rhs>
288 class 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 } // end namespace Eigen
310 
311 #endif // EIGEN_SPARSEDENSEPRODUCT_H
312