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