Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseCwiseBinaryOp.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 #ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H
26 #define EIGEN_SPARSE_CWISE_BINARY_OP_H
27
28 // Here we have to handle 3 cases:
29 //  1 - sparse op dense
30 //  2 - dense op sparse
31 //  3 - sparse op sparse
32 // We also need to implement a 4th iterator for:
33 //  4 - dense op dense
34 // Finally, we also need to distinguish between the product and other operations :
35 //                configuration      returned mode
36 //  1 - sparse op dense    product      sparse
37 //                         generic      dense
38 //  2 - dense op sparse    product      sparse
39 //                         generic      dense
40 //  3 - sparse op sparse   product      sparse
41 //                         generic      sparse
42 //  4 - dense op dense     product      dense
43 //                         generic      dense
44
45 template<typename BinaryOp, typename Lhs, typename Rhs>
46 struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
47 {
48   typedef typename ei_result_of<
49                      BinaryOp(
50                        typename Lhs::Scalar,
51                        typename Rhs::Scalar
52                      )
53                    >::type Scalar;
54   typedef typename Lhs::Nested LhsNested;
55   typedef typename Rhs::Nested RhsNested;
56   typedef typename ei_unref<LhsNested>::type _LhsNested;
57   typedef typename ei_unref<RhsNested>::type _RhsNested;
58   enum {
59     LhsCoeffReadCost = _LhsNested::CoeffReadCost,
60     RhsCoeffReadCost = _RhsNested::CoeffReadCost,
61     LhsFlags = _LhsNested::Flags,
62     RhsFlags = _RhsNested::Flags,
63     RowsAtCompileTime = Lhs::RowsAtCompileTime,
64     ColsAtCompileTime = Lhs::ColsAtCompileTime,
65     MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
66     MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
67     Flags = (int(LhsFlags) | int(RhsFlags)) & HereditaryBits,
68     CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost
69   };
70 };
71
72 template<typename BinaryOp, typename Lhs, typename Rhs>
73 class SparseCwiseBinaryOp : ei_no_assignment_operator,
74   public SparseMatrixBase<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
75 {
76   public:
77
78     class InnerIterator;
79
80     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseBinaryOp)
81     typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested;
82     typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested;
83     typedef typename ei_unref<LhsNested>::type _LhsNested;
84     typedef typename ei_unref<RhsNested>::type _RhsNested;
85
86     EIGEN_STRONG_INLINE SparseCwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
87       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
88     {
89       EIGEN_STATIC_ASSERT((_LhsNested::Flags&RowMajorBit)==(_RhsNested::Flags&RowMajorBit),
90         BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER)
91       EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret
92                            ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret)
93                            : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)),
94         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
95       // require the sizes to match
96       EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
97       ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
98     }
99
100     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
101     EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }
102
103     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
104     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
105     EIGEN_STRONG_INLINE const BinaryOp& functor() const { return m_functor; }
106
107   protected:
108     const LhsNested m_lhs;
109     const RhsNested m_rhs;
110     const BinaryOp m_functor;
111 };
112
113 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived,
114   int _LhsStorageMode = int(Lhs::Flags) & SparseBit,
115   int _RhsStorageMode = int(Rhs::Flags) & SparseBit>
116 class ei_sparse_cwise_binary_op_inner_iterator_selector;
117
118 template<typename BinaryOp, typename Lhs, typename Rhs>
119 class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
120   : public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator>
121 {
122   public:
123     typedef ei_sparse_cwise_binary_op_inner_iterator_selector<
124       BinaryOp,Lhs,Rhs, InnerIterator> Base;
125
126     EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer)
127       : Base(binOp,outer)
128     {}
129   private:
130     InnerIterator& operator=(const InnerIterator&);
131 };
132
133 /***************************************************************************
134 * Implementation of inner-iterators
135 ***************************************************************************/
136
137 // template<typename T> struct ei_func_is_conjunction { enum { ret = false }; };
138 // template<typename T> struct ei_func_is_conjunction<ei_scalar_product_op<T> > { enum { ret = true }; };
139
140 // TODO generalize the ei_scalar_product_op specialization to all conjunctions if any !
141
142 // sparse - sparse  (generic)
143 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
144 class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived, IsSparse, IsSparse>
145 {
146     typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
147     typedef typename ei_traits<CwiseBinaryXpr>::Scalar Scalar;
148     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
149     typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
150     typedef typename _LhsNested::InnerIterator LhsIterator;
151     typedef typename _RhsNested::InnerIterator RhsIterator;
152   public:
153
154     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
155       : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
156     {
157       this->operator++();
158     }
159
160     EIGEN_STRONG_INLINE Derived& operator++()
161     {
162       if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
163       {
164         m_id = m_lhsIter.index();
165         m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
166         ++m_lhsIter;
167         ++m_rhsIter;
168       }
169       else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index())))
170       {
171         m_id = m_lhsIter.index();
172         m_value = m_functor(m_lhsIter.value(), Scalar(0));
173         ++m_lhsIter;
174       }
175       else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index())))
176       {
177         m_id = m_rhsIter.index();
178         m_value = m_functor(Scalar(0), m_rhsIter.value());
179         ++m_rhsIter;
180       }
181       else
182       {
183         m_id = -1;
184       }
185       return *static_cast<Derived*>(this);
186     }
187
188     EIGEN_STRONG_INLINE Scalar value() const { return m_value; }
189
190     EIGEN_STRONG_INLINE int index() const { return m_id; }
191     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
192     EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }
193
194     EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }
195
196   protected:
197     LhsIterator m_lhsIter;
198     RhsIterator m_rhsIter;
199     const BinaryOp& m_functor;
200     Scalar m_value;
201     int m_id;
202
203   private:
204     ei_sparse_cwise_binary_op_inner_iterator_selector& operator=(const ei_sparse_cwise_binary_op_inner_iterator_selector&);
205 };
206
207 // sparse - sparse  (product)
208 template<typename T, typename Lhs, typename Rhs, typename Derived>
209 class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsSparse>
210 {
211     typedef ei_scalar_product_op<T> BinaryFunc;
212     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
213     typedef typename CwiseBinaryXpr::Scalar Scalar;
214     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
215     typedef typename _LhsNested::InnerIterator LhsIterator;
216     typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
217     typedef typename _RhsNested::InnerIterator RhsIterator;
218   public:
219
220     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
221       : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
222     {
223       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
224       {
225         if (m_lhsIter.index() < m_rhsIter.index())
226           ++m_lhsIter;
227         else
228           ++m_rhsIter;
229       }
230     }
231
232     EIGEN_STRONG_INLINE Derived& operator++()
233     {
234       ++m_lhsIter;
235       ++m_rhsIter;
236       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
237       {
238         if (m_lhsIter.index() < m_rhsIter.index())
239           ++m_lhsIter;
240         else
241           ++m_rhsIter;
242       }
243       return *static_cast<Derived*>(this);
244     }
245
246     EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); }
247
248     EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
249     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
250     EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }
251
252     EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); }
253
254   protected:
255     LhsIterator m_lhsIter;
256     RhsIterator m_rhsIter;
257     const BinaryFunc& m_functor;
258
259   private:
260     ei_sparse_cwise_binary_op_inner_iterator_selector& operator=(const ei_sparse_cwise_binary_op_inner_iterator_selector&);
261 };
262
263 // sparse - dense  (product)
264 template<typename T, typename Lhs, typename Rhs, typename Derived>
265 class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsDense>
266 {
267     typedef ei_scalar_product_op<T> BinaryFunc;
268     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
269     typedef typename CwiseBinaryXpr::Scalar Scalar;
270     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
271     typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested;
272     typedef typename _LhsNested::InnerIterator LhsIterator;
273     enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
274   public:
275
276     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
277       : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
278     {}
279
280     EIGEN_STRONG_INLINE Derived& operator++()
281     {
282       ++m_lhsIter;
283       return *static_cast<Derived*>(this);
284     }
285
286     EIGEN_STRONG_INLINE Scalar value() const
287     { return m_functor(m_lhsIter.value(),
288                        m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
289
290     EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
291     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
292     EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }
293
294     EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }
295
296   protected:
297     const RhsNested m_rhs;
298     LhsIterator m_lhsIter;
299     const BinaryFunc m_functor;
300     const int m_outer;
301
302   private:
303     ei_sparse_cwise_binary_op_inner_iterator_selector& operator=(const ei_sparse_cwise_binary_op_inner_iterator_selector&);
304 };
305
306 // sparse - dense  (product)
307 template<typename T, typename Lhs, typename Rhs, typename Derived>
308 class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsDense, IsSparse>
309 {
310     typedef ei_scalar_product_op<T> BinaryFunc;
311     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
312     typedef typename CwiseBinaryXpr::Scalar Scalar;
313     typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
314     typedef typename _RhsNested::InnerIterator RhsIterator;
315     enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit };
316   public:
317
318     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
319       : m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer)
320     {}
321
322     EIGEN_STRONG_INLINE Derived& operator++()
323     {
324       ++m_rhsIter;
325       return *static_cast<Derived*>(this);
326     }
327
328     EIGEN_STRONG_INLINE Scalar value() const
329     { return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); }
330
331     EIGEN_STRONG_INLINE int index() const { return m_rhsIter.index(); }
332     EIGEN_STRONG_INLINE int row() const { return m_rhsIter.row(); }
333     EIGEN_STRONG_INLINE int col() const { return m_rhsIter.col(); }
334
335     EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; }
336
337   protected:
338     const CwiseBinaryXpr& m_xpr;
339     RhsIterator m_rhsIter;
340     const BinaryFunc& m_functor;
341     const int m_outer;
342 };
343
344
345 /***************************************************************************
346 * Implementation of SparseMatrixBase and SparseCwise functions/operators
347 ***************************************************************************/
348
349 template<typename Derived>
350 template<typename OtherDerived>
351 EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>,
352                                  Derived, OtherDerived>
353 SparseMatrixBase<Derived>::operator-(const SparseMatrixBase<OtherDerived> &other) const
354 {
355   return SparseCwiseBinaryOp<ei_scalar_difference_op<Scalar>,
356                        Derived, OtherDerived>(derived(), other.derived());
357 }
358
359 template<typename Derived>
360 template<typename OtherDerived>
361 EIGEN_STRONG_INLINE Derived &
362 SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other)
363 {
364   return *this = derived() - other.derived();
365 }
366
367 template<typename Derived>
368 template<typename OtherDerived>
369 EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
370 SparseMatrixBase<Derived>::operator+(const SparseMatrixBase<OtherDerived> &other) const
371 {
372   return SparseCwiseBinaryOp<ei_scalar_sum_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
373 }
374
375 template<typename Derived>
376 template<typename OtherDerived>
377 EIGEN_STRONG_INLINE Derived &
378 SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other)
379 {
380   return *this = derived() + other.derived();
381 }
382
383 template<typename ExpressionType>
384 template<typename OtherDerived>
385 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
386 SparseCwise<ExpressionType>::operator*(const SparseMatrixBase<OtherDerived> &other) const
387 {
388   return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
389 }
390
391 template<typename ExpressionType>
392 template<typename OtherDerived>
393 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
394 SparseCwise<ExpressionType>::operator*(const MatrixBase<OtherDerived> &other) const
395 {
396   return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
397 }
398
399 // template<typename ExpressionType>
400 // template<typename OtherDerived>
401 // EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
402 // SparseCwise<ExpressionType>::operator/(const SparseMatrixBase<OtherDerived> &other) const
403 // {
404 //   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
405 // }
406 //
407 // template<typename ExpressionType>
408 // template<typename OtherDerived>
409 // EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
410 // SparseCwise<ExpressionType>::operator/(const MatrixBase<OtherDerived> &other) const
411 // {
412 //   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
413 // }
414
415 template<typename ExpressionType>
416 template<typename OtherDerived>
417 inline ExpressionType& SparseCwise<ExpressionType>::operator*=(const SparseMatrixBase<OtherDerived> &other)
418 {
419   return m_matrix.const_cast_derived() = _expression() * other.derived();
420 }
421
422 // template<typename ExpressionType>
423 // template<typename OtherDerived>
424 // inline ExpressionType& SparseCwise<ExpressionType>::operator/=(const SparseMatrixBase<OtherDerived> &other)
425 // {
426 //   return m_matrix.const_cast_derived() = *this / other;
427 // }
428
429 template<typename ExpressionType>
430 template<typename OtherDerived>
431 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)
432 SparseCwise<ExpressionType>::min(const SparseMatrixBase<OtherDerived> &other) const
433 {
434   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)(_expression(), other.derived());
435 }
436
437 template<typename ExpressionType>
438 template<typename OtherDerived>
439 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)
440 SparseCwise<ExpressionType>::max(const SparseMatrixBase<OtherDerived> &other) const
441 {
442   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)(_expression(), other.derived());
443 }
444
445 // template<typename Derived>
446 // template<typename CustomBinaryOp, typename OtherDerived>
447 // EIGEN_STRONG_INLINE const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
448 // SparseMatrixBase<Derived>::binaryExpr(const SparseMatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const
449 // {
450 //   return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(derived(), other.derived(), func);
451 // }
452
453 #endif // EIGEN_SPARSE_CWISE_BINARY_OP_H