1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 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_SPARSE_CWISE_BINARY_OP_H
11 #define EIGEN_SPARSE_CWISE_BINARY_OP_H
12 
13 namespace Eigen {
14 
15 // Here we have to handle 3 cases:
16 //  1 - sparse op dense
17 //  2 - dense op sparse
18 //  3 - sparse op sparse
19 // We also need to implement a 4th iterator for:
20 //  4 - dense op dense
21 // Finally, we also need to distinguish between the product and other operations :
22 //                configuration      returned mode
23 //  1 - sparse op dense    product      sparse
24 //                         generic      dense
25 //  2 - dense op sparse    product      sparse
26 //                         generic      dense
27 //  3 - sparse op sparse   product      sparse
28 //                         generic      sparse
29 //  4 - dense op dense     product      dense
30 //                         generic      dense
31 
32 namespace internal {
33 
34 template<> struct promote_storage_type<Dense,Sparse>
35 { typedef Sparse ret; };
36 
37 template<> struct promote_storage_type<Sparse,Dense>
38 { typedef Sparse ret; };
39 
40 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived,
41   typename _LhsStorageMode = typename traits<Lhs>::StorageKind,
42   typename _RhsStorageMode = typename traits<Rhs>::StorageKind>
43 class sparse_cwise_binary_op_inner_iterator_selector;
44 
45 } // end namespace internal
46 
47 template<typename BinaryOp, typename Lhs, typename Rhs>
48 class CwiseBinaryOpImpl<BinaryOp, Lhs, Rhs, Sparse>
49   : public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
50 {
51   public:
52     class InnerIterator;
53     class ReverseInnerIterator;
54     typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> Derived;
55     EIGEN_SPARSE_PUBLIC_INTERFACE(Derived)
56     CwiseBinaryOpImpl()
57     {
58       typedef typename internal::traits<Lhs>::StorageKind LhsStorageKind;
59       typedef typename internal::traits<Rhs>::StorageKind RhsStorageKind;
60       EIGEN_STATIC_ASSERT((
61                 (!internal::is_same<LhsStorageKind,RhsStorageKind>::value)
62             ||  ((Lhs::Flags&RowMajorBit) == (Rhs::Flags&RowMajorBit))),
63             THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH);
64     }
65 };
66 
67 template<typename BinaryOp, typename Lhs, typename Rhs>
68 class CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator
69   : public internal::sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs,typename CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator>
70 {
71   public:
72     typedef typename Lhs::Index Index;
73     typedef internal::sparse_cwise_binary_op_inner_iterator_selector<
74       BinaryOp,Lhs,Rhs, InnerIterator> Base;
75 
76     EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOpImpl& binOp, typename CwiseBinaryOpImpl::Index outer)
77       : Base(binOp.derived(),outer)
78     {}
79 };
80 
81 /***************************************************************************
82 * Implementation of inner-iterators
83 ***************************************************************************/
84 
85 // template<typename T> struct internal::func_is_conjunction { enum { ret = false }; };
86 // template<typename T> struct internal::func_is_conjunction<internal::scalar_product_op<T> > { enum { ret = true }; };
87 
88 // TODO generalize the internal::scalar_product_op specialization to all conjunctions if any !
89 
90 namespace internal {
91 
92 // sparse - sparse  (generic)
93 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
94 class sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived, Sparse, Sparse>
95 {
96     typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
97     typedef typename traits<CwiseBinaryXpr>::Scalar Scalar;
98     typedef typename traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
99     typedef typename traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
100     typedef typename _LhsNested::InnerIterator LhsIterator;
101     typedef typename _RhsNested::InnerIterator RhsIterator;
102     typedef typename Lhs::Index Index;
103 
104   public:
105 
106     EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer)
107       : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
108     {
109       this->operator++();
110     }
111 
112     EIGEN_STRONG_INLINE Derived& operator++()
113     {
114       if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
115       {
116         m_id = m_lhsIter.index();
117         m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
118         ++m_lhsIter;
119         ++m_rhsIter;
120       }
121       else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index())))
122       {
123         m_id = m_lhsIter.index();
124         m_value = m_functor(m_lhsIter.value(), Scalar(0));
125         ++m_lhsIter;
126       }
127       else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index())))
128       {
129         m_id = m_rhsIter.index();
130         m_value = m_functor(Scalar(0), m_rhsIter.value());
131         ++m_rhsIter;
132       }
133       else
134       {
135         m_value = 0; // this is to avoid a compilation warning
136         m_id = -1;
137       }
138       return *static_cast<Derived*>(this);
139     }
140 
141     EIGEN_STRONG_INLINE Scalar value() const { return m_value; }
142 
143     EIGEN_STRONG_INLINE Index index() const { return m_id; }
144     EIGEN_STRONG_INLINE Index row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); }
145     EIGEN_STRONG_INLINE Index col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); }
146 
147     EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }
148 
149   protected:
150     LhsIterator m_lhsIter;
151     RhsIterator m_rhsIter;
152     const BinaryOp& m_functor;
153     Scalar m_value;
154     Index m_id;
155 };
156 
157 // sparse - sparse  (product)
158 template<typename T, typename Lhs, typename Rhs, typename Derived>
159 class sparse_cwise_binary_op_inner_iterator_selector<scalar_product_op<T>, Lhs, Rhs, Derived, Sparse, Sparse>
160 {
161     typedef scalar_product_op<T> BinaryFunc;
162     typedef CwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
163     typedef typename CwiseBinaryXpr::Scalar Scalar;
164     typedef typename traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
165     typedef typename _LhsNested::InnerIterator LhsIterator;
166     typedef typename traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
167     typedef typename _RhsNested::InnerIterator RhsIterator;
168     typedef typename Lhs::Index Index;
169   public:
170 
171     EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer)
172       : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
173     {
174       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
175       {
176         if (m_lhsIter.index() < m_rhsIter.index())
177           ++m_lhsIter;
178         else
179           ++m_rhsIter;
180       }
181     }
182 
183     EIGEN_STRONG_INLINE Derived& operator++()
184     {
185       ++m_lhsIter;
186       ++m_rhsIter;
187       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
188       {
189         if (m_lhsIter.index() < m_rhsIter.index())
190           ++m_lhsIter;
191         else
192           ++m_rhsIter;
193       }
194       return *static_cast<Derived*>(this);
195     }
196 
197     EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); }
198 
199     EIGEN_STRONG_INLINE Index index() const { return m_lhsIter.index(); }
200     EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); }
201     EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); }
202 
203     EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); }
204 
205   protected:
206     LhsIterator m_lhsIter;
207     RhsIterator m_rhsIter;
208     const BinaryFunc& m_functor;
209 };
210 
211 // sparse - dense  (product)
212 template<typename T, typename Lhs, typename Rhs, typename Derived>
213 class sparse_cwise_binary_op_inner_iterator_selector<scalar_product_op<T>, Lhs, Rhs, Derived, Sparse, Dense>
214 {
215     typedef scalar_product_op<T> BinaryFunc;
216     typedef CwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
217     typedef typename CwiseBinaryXpr::Scalar Scalar;
218     typedef typename traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
219     typedef typename traits<CwiseBinaryXpr>::RhsNested RhsNested;
220     typedef typename _LhsNested::InnerIterator LhsIterator;
221     typedef typename Lhs::Index Index;
222     enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
223   public:
224 
225     EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer)
226       : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
227     {}
228 
229     EIGEN_STRONG_INLINE Derived& operator++()
230     {
231       ++m_lhsIter;
232       return *static_cast<Derived*>(this);
233     }
234 
235     EIGEN_STRONG_INLINE Scalar value() const
236     { return m_functor(m_lhsIter.value(),
237                        m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
238 
239     EIGEN_STRONG_INLINE Index index() const { return m_lhsIter.index(); }
240     EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); }
241     EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); }
242 
243     EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }
244 
245   protected:
246     RhsNested m_rhs;
247     LhsIterator m_lhsIter;
248     const BinaryFunc m_functor;
249     const Index m_outer;
250 };
251 
252 // sparse - dense  (product)
253 template<typename T, typename Lhs, typename Rhs, typename Derived>
254 class sparse_cwise_binary_op_inner_iterator_selector<scalar_product_op<T>, Lhs, Rhs, Derived, Dense, Sparse>
255 {
256     typedef scalar_product_op<T> BinaryFunc;
257     typedef CwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
258     typedef typename CwiseBinaryXpr::Scalar Scalar;
259     typedef typename traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
260     typedef typename _RhsNested::InnerIterator RhsIterator;
261     typedef typename Lhs::Index Index;
262 
263     enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit };
264   public:
265 
266     EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer)
267       : m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer)
268     {}
269 
270     EIGEN_STRONG_INLINE Derived& operator++()
271     {
272       ++m_rhsIter;
273       return *static_cast<Derived*>(this);
274     }
275 
276     EIGEN_STRONG_INLINE Scalar value() const
277     { return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); }
278 
279     EIGEN_STRONG_INLINE Index index() const { return m_rhsIter.index(); }
280     EIGEN_STRONG_INLINE Index row() const { return m_rhsIter.row(); }
281     EIGEN_STRONG_INLINE Index col() const { return m_rhsIter.col(); }
282 
283     EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; }
284 
285   protected:
286     const CwiseBinaryXpr& m_xpr;
287     RhsIterator m_rhsIter;
288     const BinaryFunc& m_functor;
289     const Index m_outer;
290 };
291 
292 } // end namespace internal
293 
294 /***************************************************************************
295 * Implementation of SparseMatrixBase and SparseCwise functions/operators
296 ***************************************************************************/
297 
298 template<typename Derived>
299 template<typename OtherDerived>
300 EIGEN_STRONG_INLINE Derived &
301 SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other)
302 {
303   return *this = derived() - other.derived();
304 }
305 
306 template<typename Derived>
307 template<typename OtherDerived>
308 EIGEN_STRONG_INLINE Derived &
309 SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other)
310 {
311   return *this = derived() + other.derived();
312 }
313 
314 template<typename Derived>
315 template<typename OtherDerived>
316 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
317 SparseMatrixBase<Derived>::cwiseProduct(const MatrixBase<OtherDerived> &other) const
318 {
319   return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(derived(), other.derived());
320 }
321 
322 } // end namespace Eigen
323 
324 #endif // EIGEN_SPARSE_CWISE_BINARY_OP_H
325