src/SparseCore/SparseCwiseBinaryOp.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H
0011 #define EIGEN_SPARSE_CWISE_BINARY_OP_H
0012
0013 namespace Eigen {
0014
0015 // Here we have to handle 3 cases:
0016 //  1 - sparse op dense
0017 //  2 - dense op sparse
0018 //  3 - sparse op sparse
0019 // We also need to implement a 4th iterator for:
0020 //  4 - dense op dense
0021 // Finally, we also need to distinguish between the product and other operations :
0022 //                configuration      returned mode
0023 //  1 - sparse op dense    product      sparse
0024 //                         generic      dense
0025 //  2 - dense op sparse    product      sparse
0026 //                         generic      dense
0027 //  3 - sparse op sparse   product      sparse
0028 //                         generic      sparse
0029 //  4 - dense op dense     product      dense
0030 //                         generic      dense
0031 //
0032 // TODO to ease compiler job, we could specialize product/quotient with a scalar
0033 //      and fallback to cwise-unary evaluator using bind1st_op and bind2nd_op.
0034
0035 template<typename BinaryOp, typename Lhs, typename Rhs>
0036 class CwiseBinaryOpImpl<BinaryOp, Lhs, Rhs, Sparse>
0037   : public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
0038 {
0039   public:
0040     typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> Derived;
0041     typedef SparseMatrixBase<Derived> Base;
0042     EIGEN_SPARSE_PUBLIC_INTERFACE(Derived)
0043     CwiseBinaryOpImpl()
0044     {
0045       EIGEN_STATIC_ASSERT((
0046                 (!internal::is_same<typename internal::traits<Lhs>::StorageKind,
0047                                     typename internal::traits<Rhs>::StorageKind>::value)
0048             ||  ((internal::evaluator<Lhs>::Flags&RowMajorBit) == (internal::evaluator<Rhs>::Flags&RowMajorBit))),
0049             THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH);
0050     }
0051 };
0052
0053 namespace internal {
0054
0055
0056 // Generic "sparse OP sparse"
0057 template<typename XprType> struct binary_sparse_evaluator;
0058
0059 template<typename BinaryOp, typename Lhs, typename Rhs>
0060 struct binary_evaluator<CwiseBinaryOp<BinaryOp, Lhs, Rhs>, IteratorBased, IteratorBased>
0061   : evaluator_base<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
0062 {
0063 protected:
0064   typedef typename evaluator<Lhs>::InnerIterator  LhsIterator;
0065   typedef typename evaluator<Rhs>::InnerIterator  RhsIterator;
0066   typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> XprType;
0067   typedef typename traits<XprType>::Scalar Scalar;
0068   typedef typename XprType::StorageIndex StorageIndex;
0069 public:
0070
0071   class InnerIterator
0072   {
0073   public:
0074
0075     EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer)
0076       : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor)
0077     {
0078       this->operator++();
0079     }
0080
0081     EIGEN_STRONG_INLINE InnerIterator& operator++()
0082     {
0083       if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
0084       {
0085         m_id = m_lhsIter.index();
0086         m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
0087         ++m_lhsIter;
0088         ++m_rhsIter;
0089       }
0090       else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index())))
0091       {
0092         m_id = m_lhsIter.index();
0093         m_value = m_functor(m_lhsIter.value(), Scalar(0));
0094         ++m_lhsIter;
0095       }
0096       else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index())))
0097       {
0098         m_id = m_rhsIter.index();
0099         m_value = m_functor(Scalar(0), m_rhsIter.value());
0100         ++m_rhsIter;
0101       }
0102       else
0103       {
0104         m_value = Scalar(0); // this is to avoid a compilation warning
0105         m_id = -1;
0106       }
0107       return *this;
0108     }
0109
0110     EIGEN_STRONG_INLINE Scalar value() const { return m_value; }
0111
0112     EIGEN_STRONG_INLINE StorageIndex index() const { return m_id; }
0113     EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); }
0114     EIGEN_STRONG_INLINE Index row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); }
0115     EIGEN_STRONG_INLINE Index col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); }
0116
0117     EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }
0118
0119   protected:
0120     LhsIterator m_lhsIter;
0121     RhsIterator m_rhsIter;
0122     const BinaryOp& m_functor;
0123     Scalar m_value;
0124     StorageIndex m_id;
0125   };
0126
0127
0128   enum {
0129     CoeffReadCost = int(evaluator<Lhs>::CoeffReadCost) + int(evaluator<Rhs>::CoeffReadCost) + int(functor_traits<BinaryOp>::Cost),
0130     Flags = XprType::Flags
0131   };
0132
0133   explicit binary_evaluator(const XprType& xpr)
0134     : m_functor(xpr.functor()),
0135       m_lhsImpl(xpr.lhs()),
0136       m_rhsImpl(xpr.rhs())
0137   {
0138     EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits<BinaryOp>::Cost);
0139     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0140   }
0141
0142   inline Index nonZerosEstimate() const {
0143     return m_lhsImpl.nonZerosEstimate() + m_rhsImpl.nonZerosEstimate();
0144   }
0145
0146 protected:
0147   const BinaryOp m_functor;
0148   evaluator<Lhs> m_lhsImpl;
0149   evaluator<Rhs> m_rhsImpl;
0150 };
0151
0152 // dense op sparse
0153 template<typename BinaryOp, typename Lhs, typename Rhs>
0154 struct binary_evaluator<CwiseBinaryOp<BinaryOp, Lhs, Rhs>, IndexBased, IteratorBased>
0155   : evaluator_base<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
0156 {
0157 protected:
0158   typedef typename evaluator<Rhs>::InnerIterator  RhsIterator;
0159   typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> XprType;
0160   typedef typename traits<XprType>::Scalar Scalar;
0161   typedef typename XprType::StorageIndex StorageIndex;
0162 public:
0163
0164   class InnerIterator
0165   {
0166     enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit };
0167   public:
0168
0169     EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer)
0170       : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_value(0), m_id(-1), m_innerSize(aEval.m_expr.rhs().innerSize())
0171     {
0172       this->operator++();
0173     }
0174
0175     EIGEN_STRONG_INLINE InnerIterator& operator++()
0176     {
0177       ++m_id;
0178       if(m_id<m_innerSize)
0179       {
0180         Scalar lhsVal = m_lhsEval.coeff(IsRowMajor?m_rhsIter.outer():m_id,
0181                                         IsRowMajor?m_id:m_rhsIter.outer());
0182         if(m_rhsIter && m_rhsIter.index()==m_id)
0183         {
0184           m_value = m_functor(lhsVal, m_rhsIter.value());
0185           ++m_rhsIter;
0186         }
0187         else
0188           m_value = m_functor(lhsVal, Scalar(0));
0189       }
0190
0191       return *this;
0192     }
0193
0194     EIGEN_STRONG_INLINE Scalar value() const { eigen_internal_assert(m_id<m_innerSize); return m_value; }
0195
0196     EIGEN_STRONG_INLINE StorageIndex index() const { return m_id; }
0197     EIGEN_STRONG_INLINE Index outer() const { return m_rhsIter.outer(); }
0198     EIGEN_STRONG_INLINE Index row() const { return IsRowMajor ? m_rhsIter.outer() : m_id; }
0199     EIGEN_STRONG_INLINE Index col() const { return IsRowMajor ? m_id : m_rhsIter.outer(); }
0200
0201     EIGEN_STRONG_INLINE operator bool() const { return m_id<m_innerSize; }
0202
0203   protected:
0204     const evaluator<Lhs> &m_lhsEval;
0205     RhsIterator m_rhsIter;
0206     const BinaryOp& m_functor;
0207     Scalar m_value;
0208     StorageIndex m_id;
0209     StorageIndex m_innerSize;
0210   };
0211
0212
0213   enum {
0214     CoeffReadCost = int(evaluator<Lhs>::CoeffReadCost) + int(evaluator<Rhs>::CoeffReadCost) + int(functor_traits<BinaryOp>::Cost),
0215     Flags = XprType::Flags
0216   };
0217
0218   explicit binary_evaluator(const XprType& xpr)
0219     : m_functor(xpr.functor()),
0220       m_lhsImpl(xpr.lhs()),
0221       m_rhsImpl(xpr.rhs()),
0222       m_expr(xpr)
0223   {
0224     EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits<BinaryOp>::Cost);
0225     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0226   }
0227
0228   inline Index nonZerosEstimate() const {
0229     return m_expr.size();
0230   }
0231
0232 protected:
0233   const BinaryOp m_functor;
0234   evaluator<Lhs> m_lhsImpl;
0235   evaluator<Rhs> m_rhsImpl;
0236   const XprType &m_expr;
0237 };
0238
0239 // sparse op dense
0240 template<typename BinaryOp, typename Lhs, typename Rhs>
0241 struct binary_evaluator<CwiseBinaryOp<BinaryOp, Lhs, Rhs>, IteratorBased, IndexBased>
0242   : evaluator_base<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
0243 {
0244 protected:
0245   typedef typename evaluator<Lhs>::InnerIterator  LhsIterator;
0246   typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> XprType;
0247   typedef typename traits<XprType>::Scalar Scalar;
0248   typedef typename XprType::StorageIndex StorageIndex;
0249 public:
0250
0251   class InnerIterator
0252   {
0253     enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
0254   public:
0255
0256     EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer)
0257       : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_value(0), m_id(-1), m_innerSize(aEval.m_expr.lhs().innerSize())
0258     {
0259       this->operator++();
0260     }
0261
0262     EIGEN_STRONG_INLINE InnerIterator& operator++()
0263     {
0264       ++m_id;
0265       if(m_id<m_innerSize)
0266       {
0267         Scalar rhsVal = m_rhsEval.coeff(IsRowMajor?m_lhsIter.outer():m_id,
0268                                         IsRowMajor?m_id:m_lhsIter.outer());
0269         if(m_lhsIter && m_lhsIter.index()==m_id)
0270         {
0271           m_value = m_functor(m_lhsIter.value(), rhsVal);
0272           ++m_lhsIter;
0273         }
0274         else
0275           m_value = m_functor(Scalar(0),rhsVal);
0276       }
0277
0278       return *this;
0279     }
0280
0281     EIGEN_STRONG_INLINE Scalar value() const { eigen_internal_assert(m_id<m_innerSize); return m_value; }
0282
0283     EIGEN_STRONG_INLINE StorageIndex index() const { return m_id; }
0284     EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); }
0285     EIGEN_STRONG_INLINE Index row() const { return IsRowMajor ? m_lhsIter.outer() : m_id; }
0286     EIGEN_STRONG_INLINE Index col() const { return IsRowMajor ? m_id : m_lhsIter.outer(); }
0287
0288     EIGEN_STRONG_INLINE operator bool() const { return m_id<m_innerSize; }
0289
0290   protected:
0291     LhsIterator m_lhsIter;
0292     const evaluator<Rhs> &m_rhsEval;
0293     const BinaryOp& m_functor;
0294     Scalar m_value;
0295     StorageIndex m_id;
0296     StorageIndex m_innerSize;
0297   };
0298
0299
0300   enum {
0301     CoeffReadCost = int(evaluator<Lhs>::CoeffReadCost) + int(evaluator<Rhs>::CoeffReadCost) + int(functor_traits<BinaryOp>::Cost),
0302     Flags = XprType::Flags
0303   };
0304
0305   explicit binary_evaluator(const XprType& xpr)
0306     : m_functor(xpr.functor()),
0307       m_lhsImpl(xpr.lhs()),
0308       m_rhsImpl(xpr.rhs()),
0309       m_expr(xpr)
0310   {
0311     EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits<BinaryOp>::Cost);
0312     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0313   }
0314
0315   inline Index nonZerosEstimate() const {
0316     return m_expr.size();
0317   }
0318
0319 protected:
0320   const BinaryOp m_functor;
0321   evaluator<Lhs> m_lhsImpl;
0322   evaluator<Rhs> m_rhsImpl;
0323   const XprType &m_expr;
0324 };
0325
0326 template<typename T,
0327          typename LhsKind   = typename evaluator_traits<typename T::Lhs>::Kind,
0328          typename RhsKind   = typename evaluator_traits<typename T::Rhs>::Kind,
0329          typename LhsScalar = typename traits<typename T::Lhs>::Scalar,
0330          typename RhsScalar = typename traits<typename T::Rhs>::Scalar> struct sparse_conjunction_evaluator;
0331
0332 // "sparse .* sparse"
0333 template<typename T1, typename T2, typename Lhs, typename Rhs>
0334 struct binary_evaluator<CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs>, IteratorBased, IteratorBased>
0335   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs> >
0336 {
0337   typedef CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs> XprType;
0338   typedef sparse_conjunction_evaluator<XprType> Base;
0339   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0340 };
0341 // "dense .* sparse"
0342 template<typename T1, typename T2, typename Lhs, typename Rhs>
0343 struct binary_evaluator<CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs>, IndexBased, IteratorBased>
0344   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs> >
0345 {
0346   typedef CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs> XprType;
0347   typedef sparse_conjunction_evaluator<XprType> Base;
0348   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0349 };
0350 // "sparse .* dense"
0351 template<typename T1, typename T2, typename Lhs, typename Rhs>
0352 struct binary_evaluator<CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs>, IteratorBased, IndexBased>
0353   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs> >
0354 {
0355   typedef CwiseBinaryOp<scalar_product_op<T1,T2>, Lhs, Rhs> XprType;
0356   typedef sparse_conjunction_evaluator<XprType> Base;
0357   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0358 };
0359
0360 // "sparse ./ dense"
0361 template<typename T1, typename T2, typename Lhs, typename Rhs>
0362 struct binary_evaluator<CwiseBinaryOp<scalar_quotient_op<T1,T2>, Lhs, Rhs>, IteratorBased, IndexBased>
0363   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_quotient_op<T1,T2>, Lhs, Rhs> >
0364 {
0365   typedef CwiseBinaryOp<scalar_quotient_op<T1,T2>, Lhs, Rhs> XprType;
0366   typedef sparse_conjunction_evaluator<XprType> Base;
0367   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0368 };
0369
0370 // "sparse && sparse"
0371 template<typename Lhs, typename Rhs>
0372 struct binary_evaluator<CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs>, IteratorBased, IteratorBased>
0373   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs> >
0374 {
0375   typedef CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs> XprType;
0376   typedef sparse_conjunction_evaluator<XprType> Base;
0377   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0378 };
0379 // "dense && sparse"
0380 template<typename Lhs, typename Rhs>
0381 struct binary_evaluator<CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs>, IndexBased, IteratorBased>
0382   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs> >
0383 {
0384   typedef CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs> XprType;
0385   typedef sparse_conjunction_evaluator<XprType> Base;
0386   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0387 };
0388 // "sparse && dense"
0389 template<typename Lhs, typename Rhs>
0390 struct binary_evaluator<CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs>, IteratorBased, IndexBased>
0391   : sparse_conjunction_evaluator<CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs> >
0392 {
0393   typedef CwiseBinaryOp<scalar_boolean_and_op, Lhs, Rhs> XprType;
0394   typedef sparse_conjunction_evaluator<XprType> Base;
0395   explicit binary_evaluator(const XprType& xpr) : Base(xpr) {}
0396 };
0397
0398 // "sparse ^ sparse"
0399 template<typename XprType>
0400 struct sparse_conjunction_evaluator<XprType, IteratorBased, IteratorBased>
0401   : evaluator_base<XprType>
0402 {
0403 protected:
0404   typedef typename XprType::Functor BinaryOp;
0405   typedef typename XprType::Lhs LhsArg;
0406   typedef typename XprType::Rhs RhsArg;
0407   typedef typename evaluator<LhsArg>::InnerIterator  LhsIterator;
0408   typedef typename evaluator<RhsArg>::InnerIterator  RhsIterator;
0409   typedef typename XprType::StorageIndex StorageIndex;
0410   typedef typename traits<XprType>::Scalar Scalar;
0411 public:
0412
0413   class InnerIterator
0414   {
0415   public:
0416
0417     EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer)
0418       : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor)
0419     {
0420       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
0421       {
0422         if (m_lhsIter.index() < m_rhsIter.index())
0423           ++m_lhsIter;
0424         else
0425           ++m_rhsIter;
0426       }
0427     }
0428
0429     EIGEN_STRONG_INLINE InnerIterator& operator++()
0430     {
0431       ++m_lhsIter;
0432       ++m_rhsIter;
0433       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
0434       {
0435         if (m_lhsIter.index() < m_rhsIter.index())
0436           ++m_lhsIter;
0437         else
0438           ++m_rhsIter;
0439       }
0440       return *this;
0441     }
0442
0443     EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); }
0444
0445     EIGEN_STRONG_INLINE StorageIndex index() const { return m_lhsIter.index(); }
0446     EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); }
0447     EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); }
0448     EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); }
0449
0450     EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); }
0451
0452   protected:
0453     LhsIterator m_lhsIter;
0454     RhsIterator m_rhsIter;
0455     const BinaryOp& m_functor;
0456   };
0457
0458
0459   enum {
0460     CoeffReadCost = int(evaluator<LhsArg>::CoeffReadCost) + int(evaluator<RhsArg>::CoeffReadCost) + int(functor_traits<BinaryOp>::Cost),
0461     Flags = XprType::Flags
0462   };
0463
0464   explicit sparse_conjunction_evaluator(const XprType& xpr)
0465     : m_functor(xpr.functor()),
0466       m_lhsImpl(xpr.lhs()),
0467       m_rhsImpl(xpr.rhs())
0468   {
0469     EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits<BinaryOp>::Cost);
0470     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0471   }
0472
0473   inline Index nonZerosEstimate() const {
0474     return (std::min)(m_lhsImpl.nonZerosEstimate(), m_rhsImpl.nonZerosEstimate());
0475   }
0476
0477 protected:
0478   const BinaryOp m_functor;
0479   evaluator<LhsArg> m_lhsImpl;
0480   evaluator<RhsArg> m_rhsImpl;
0481 };
0482
0483 // "dense ^ sparse"
0484 template<typename XprType>
0485 struct sparse_conjunction_evaluator<XprType, IndexBased, IteratorBased>
0486   : evaluator_base<XprType>
0487 {
0488 protected:
0489   typedef typename XprType::Functor BinaryOp;
0490   typedef typename XprType::Lhs LhsArg;
0491   typedef typename XprType::Rhs RhsArg;
0492   typedef evaluator<LhsArg> LhsEvaluator;
0493   typedef typename evaluator<RhsArg>::InnerIterator  RhsIterator;
0494   typedef typename XprType::StorageIndex StorageIndex;
0495   typedef typename traits<XprType>::Scalar Scalar;
0496 public:
0497
0498   class InnerIterator
0499   {
0500     enum { IsRowMajor = (int(RhsArg::Flags)&RowMajorBit)==RowMajorBit };
0501
0502   public:
0503
0504     EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer)
0505       : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_outer(outer)
0506     {}
0507
0508     EIGEN_STRONG_INLINE InnerIterator& operator++()
0509     {
0510       ++m_rhsIter;
0511       return *this;
0512     }
0513
0514     EIGEN_STRONG_INLINE Scalar value() const
0515     { return m_functor(m_lhsEval.coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); }
0516
0517     EIGEN_STRONG_INLINE StorageIndex index() const { return m_rhsIter.index(); }
0518     EIGEN_STRONG_INLINE Index outer() const { return m_rhsIter.outer(); }
0519     EIGEN_STRONG_INLINE Index row() const { return m_rhsIter.row(); }
0520     EIGEN_STRONG_INLINE Index col() const { return m_rhsIter.col(); }
0521
0522     EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; }
0523
0524   protected:
0525     const LhsEvaluator &m_lhsEval;
0526     RhsIterator m_rhsIter;
0527     const BinaryOp& m_functor;
0528     const Index m_outer;
0529   };
0530
0531
0532   enum {
0533     CoeffReadCost = int(evaluator<LhsArg>::CoeffReadCost) + int(evaluator<RhsArg>::CoeffReadCost) + int(functor_traits<BinaryOp>::Cost),
0534     Flags = XprType::Flags
0535   };
0536
0537   explicit sparse_conjunction_evaluator(const XprType& xpr)
0538     : m_functor(xpr.functor()),
0539       m_lhsImpl(xpr.lhs()),
0540       m_rhsImpl(xpr.rhs())
0541   {
0542     EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits<BinaryOp>::Cost);
0543     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0544   }
0545
0546   inline Index nonZerosEstimate() const {
0547     return m_rhsImpl.nonZerosEstimate();
0548   }
0549
0550 protected:
0551   const BinaryOp m_functor;
0552   evaluator<LhsArg> m_lhsImpl;
0553   evaluator<RhsArg> m_rhsImpl;
0554 };
0555
0556 // "sparse ^ dense"
0557 template<typename XprType>
0558 struct sparse_conjunction_evaluator<XprType, IteratorBased, IndexBased>
0559   : evaluator_base<XprType>
0560 {
0561 protected:
0562   typedef typename XprType::Functor BinaryOp;
0563   typedef typename XprType::Lhs LhsArg;
0564   typedef typename XprType::Rhs RhsArg;
0565   typedef typename evaluator<LhsArg>::InnerIterator LhsIterator;
0566   typedef evaluator<RhsArg> RhsEvaluator;
0567   typedef typename XprType::StorageIndex StorageIndex;
0568   typedef typename traits<XprType>::Scalar Scalar;
0569 public:
0570
0571   class InnerIterator
0572   {
0573     enum { IsRowMajor = (int(LhsArg::Flags)&RowMajorBit)==RowMajorBit };
0574
0575   public:
0576
0577     EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer)
0578       : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_outer(outer)
0579     {}
0580
0581     EIGEN_STRONG_INLINE InnerIterator& operator++()
0582     {
0583       ++m_lhsIter;
0584       return *this;
0585     }
0586
0587     EIGEN_STRONG_INLINE Scalar value() const
0588     { return m_functor(m_lhsIter.value(),
0589                        m_rhsEval.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
0590
0591     EIGEN_STRONG_INLINE StorageIndex index() const { return m_lhsIter.index(); }
0592     EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); }
0593     EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); }
0594     EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); }
0595
0596     EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }
0597
0598   protected:
0599     LhsIterator m_lhsIter;
0600     const evaluator<RhsArg> &m_rhsEval;
0601     const BinaryOp& m_functor;
0602     const Index m_outer;
0603   };
0604
0605
0606   enum {
0607     CoeffReadCost = int(evaluator<LhsArg>::CoeffReadCost) + int(evaluator<RhsArg>::CoeffReadCost) + int(functor_traits<BinaryOp>::Cost),
0608     Flags = XprType::Flags
0609   };
0610
0611   explicit sparse_conjunction_evaluator(const XprType& xpr)
0612     : m_functor(xpr.functor()),
0613       m_lhsImpl(xpr.lhs()),
0614       m_rhsImpl(xpr.rhs())
0615   {
0616     EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits<BinaryOp>::Cost);
0617     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0618   }
0619
0620   inline Index nonZerosEstimate() const {
0621     return m_lhsImpl.nonZerosEstimate();
0622   }
0623
0624 protected:
0625   const BinaryOp m_functor;
0626   evaluator<LhsArg> m_lhsImpl;
0627   evaluator<RhsArg> m_rhsImpl;
0628 };
0629
0630 }
0631
0632 /***************************************************************************
0633 * Implementation of SparseMatrixBase and SparseCwise functions/operators
0634 ***************************************************************************/
0635
0636 template<typename Derived>
0637 template<typename OtherDerived>
0638 Derived& SparseMatrixBase<Derived>::operator+=(const EigenBase<OtherDerived> &other)
0639 {
0640   call_assignment(derived(), other.derived(), internal::add_assign_op<Scalar,typename OtherDerived::Scalar>());
0641   return derived();
0642 }
0643
0644 template<typename Derived>
0645 template<typename OtherDerived>
0646 Derived& SparseMatrixBase<Derived>::operator-=(const EigenBase<OtherDerived> &other)
0647 {
0648   call_assignment(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
0649   return derived();
0650 }
0651
0652 template<typename Derived>
0653 template<typename OtherDerived>
0654 EIGEN_STRONG_INLINE Derived &
0655 SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other)
0656 {
0657   return derived() = derived() - other.derived();
0658 }
0659
0660 template<typename Derived>
0661 template<typename OtherDerived>
0662 EIGEN_STRONG_INLINE Derived &
0663 SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other)
0664 {
0665   return derived() = derived() + other.derived();
0666 }
0667
0668 template<typename Derived>
0669 template<typename OtherDerived>
0670 Derived& SparseMatrixBase<Derived>::operator+=(const DiagonalBase<OtherDerived>& other)
0671 {
0672   call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename OtherDerived::Scalar>());
0673   return derived();
0674 }
0675
0676 template<typename Derived>
0677 template<typename OtherDerived>
0678 Derived& SparseMatrixBase<Derived>::operator-=(const DiagonalBase<OtherDerived>& other)
0679 {
0680   call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename OtherDerived::Scalar>());
0681   return derived();
0682 }
0683
0684 template<typename Derived>
0685 template<typename OtherDerived>
0686 EIGEN_STRONG_INLINE const typename SparseMatrixBase<Derived>::template CwiseProductDenseReturnType<OtherDerived>::Type
0687 SparseMatrixBase<Derived>::cwiseProduct(const MatrixBase<OtherDerived> &other) const
0688 {
0689   return typename CwiseProductDenseReturnType<OtherDerived>::Type(derived(), other.derived());
0690 }
0691
0692 template<typename DenseDerived, typename SparseDerived>
0693 EIGEN_STRONG_INLINE const CwiseBinaryOp<internal::scalar_sum_op<typename DenseDerived::Scalar,typename SparseDerived::Scalar>, const DenseDerived, const SparseDerived>
0694 operator+(const MatrixBase<DenseDerived> &a, const SparseMatrixBase<SparseDerived> &b)
0695 {
0696   return CwiseBinaryOp<internal::scalar_sum_op<typename DenseDerived::Scalar,typename SparseDerived::Scalar>, const DenseDerived, const SparseDerived>(a.derived(), b.derived());
0697 }
0698
0699 template<typename SparseDerived, typename DenseDerived>
0700 EIGEN_STRONG_INLINE const CwiseBinaryOp<internal::scalar_sum_op<typename SparseDerived::Scalar,typename DenseDerived::Scalar>, const SparseDerived, const DenseDerived>
0701 operator+(const SparseMatrixBase<SparseDerived> &a, const MatrixBase<DenseDerived> &b)
0702 {
0703   return CwiseBinaryOp<internal::scalar_sum_op<typename SparseDerived::Scalar,typename DenseDerived::Scalar>, const SparseDerived, const DenseDerived>(a.derived(), b.derived());
0704 }
0705
0706 template<typename DenseDerived, typename SparseDerived>
0707 EIGEN_STRONG_INLINE const CwiseBinaryOp<internal::scalar_difference_op<typename DenseDerived::Scalar,typename SparseDerived::Scalar>, const DenseDerived, const SparseDerived>
0708 operator-(const MatrixBase<DenseDerived> &a, const SparseMatrixBase<SparseDerived> &b)
0709 {
0710   return CwiseBinaryOp<internal::scalar_difference_op<typename DenseDerived::Scalar,typename SparseDerived::Scalar>, const DenseDerived, const SparseDerived>(a.derived(), b.derived());
0711 }
0712
0713 template<typename SparseDerived, typename DenseDerived>
0714 EIGEN_STRONG_INLINE const CwiseBinaryOp<internal::scalar_difference_op<typename SparseDerived::Scalar,typename DenseDerived::Scalar>, const SparseDerived, const DenseDerived>
0715 operator-(const SparseMatrixBase<SparseDerived> &a, const MatrixBase<DenseDerived> &b)
0716 {
0717   return CwiseBinaryOp<internal::scalar_difference_op<typename SparseDerived::Scalar,typename DenseDerived::Scalar>, const SparseDerived, const DenseDerived>(a.derived(), b.derived());
0718 }
0719
0720 } // end namespace Eigen
0721
0722 #endif // EIGEN_SPARSE_CWISE_BINARY_OP_H