Eigen/src/Sparse/SparseCwiseBinaryOp.h - mirror - Git at Google

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra. Eigen itself is part of the KDE project.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
 // version 3 of the License, or (at your option) any later version.
 //
 // Alternatively, you can redistribute it and/or
 // modify it under the terms of the GNU General Public License as
 // published by the Free Software Foundation; either version 2 of
 // the License, or (at your option) any later version.
 //
 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
 // GNU General Public License for more details.
 //
 // You should have received a copy of the GNU Lesser General Public
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.

 #ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H
 #define EIGEN_SPARSE_CWISE_BINARY_OP_H

 // Here we have to handle 3 cases:
 //  1 - sparse op dense
 //  2 - dense op sparse
 //  3 - sparse op sparse
 // We also need to implement a 4th iterator for:
 //  4 - dense op dense
 // Finally, we also need to distinguish between the product and other operations :
 //                configuration      returned mode
 //  1 - sparse op dense    product      sparse
 //                         generic      dense
 //  2 - dense op sparse    product      sparse
 //                         generic      dense
 //  3 - sparse op sparse   product      sparse
 //                         generic      sparse
 //  4 - dense op dense     product      dense
 //                         generic      dense

 template<typename BinaryOp, typename Lhs, typename Rhs>
 struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
 {
   typedef typename ei_result_of<
                      BinaryOp(
                        typename Lhs::Scalar,
                        typename Rhs::Scalar
                      )
                    >::type Scalar;
   typedef typename Lhs::Nested LhsNested;
   typedef typename Rhs::Nested RhsNested;
   typedef typename ei_unref<LhsNested>::type _LhsNested;
   typedef typename ei_unref<RhsNested>::type _RhsNested;
   enum {
     LhsCoeffReadCost = _LhsNested::CoeffReadCost,
     RhsCoeffReadCost = _RhsNested::CoeffReadCost,
     LhsFlags = _LhsNested::Flags,
     RhsFlags = _RhsNested::Flags,
     RowsAtCompileTime = Lhs::RowsAtCompileTime,
     ColsAtCompileTime = Lhs::ColsAtCompileTime,
     MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
     Flags = (int(LhsFlags) | int(RhsFlags)) & HereditaryBits,
     CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost
   };
 };

 template<typename BinaryOp, typename Lhs, typename Rhs>
 class SparseCwiseBinaryOp : ei_no_assignment_operator,
   public SparseMatrixBase<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
 {
   public:

     class InnerIterator;

     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseBinaryOp)
     typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested;
     typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested;
     typedef typename ei_unref<LhsNested>::type _LhsNested;
     typedef typename ei_unref<RhsNested>::type _RhsNested;

     EIGEN_STRONG_INLINE SparseCwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
     {
       EIGEN_STATIC_ASSERT((_LhsNested::Flags&RowMajorBit)==(_RhsNested::Flags&RowMajorBit),
         BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER)
       EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret
                            ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret)
                            : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)),
         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
       // require the sizes to match
       EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
       ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
     }

     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
     EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }

     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
     EIGEN_STRONG_INLINE const BinaryOp& functor() const { return m_functor; }

   protected:
     const LhsNested m_lhs;
     const RhsNested m_rhs;
     const BinaryOp m_functor;
 };

 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived,
   int _LhsStorageMode = int(Lhs::Flags) & SparseBit,
   int _RhsStorageMode = int(Rhs::Flags) & SparseBit>
 class ei_sparse_cwise_binary_op_inner_iterator_selector;

 template<typename BinaryOp, typename Lhs, typename Rhs>
 class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
   : public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator>
 {
   public:
     typedef ei_sparse_cwise_binary_op_inner_iterator_selector<
       BinaryOp,Lhs,Rhs, InnerIterator> Base;

     EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer)
       : Base(binOp,outer)
     {}
 };

 /***************************************************************************
 * Implementation of inner-iterators
 ***************************************************************************/

 // template<typename T> struct ei_func_is_conjunction { enum { ret = false }; };
 // template<typename T> struct ei_func_is_conjunction<ei_scalar_product_op<T> > { enum { ret = true }; };

 // TODO generalize the ei_scalar_product_op specialization to all conjunctions if any !

 // sparse - sparse  (generic)
 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
 class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived, IsSparse, IsSparse>
 {
     typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
     typedef typename ei_traits<CwiseBinaryXpr>::Scalar Scalar;
     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
     typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
     typedef typename _LhsNested::InnerIterator LhsIterator;
     typedef typename _RhsNested::InnerIterator RhsIterator;
   public:

     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
       : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
     {
       this->operator++();
     }

     EIGEN_STRONG_INLINE Derived& operator++()
     {
       if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
       {
         m_id = m_lhsIter.index();
         m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
         ++m_lhsIter;
         ++m_rhsIter;
       }
       else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index())))
       {
         m_id = m_lhsIter.index();
         m_value = m_functor(m_lhsIter.value(), Scalar(0));
         ++m_lhsIter;
       }
       else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index())))
       {
         m_id = m_rhsIter.index();
         m_value = m_functor(Scalar(0), m_rhsIter.value());
         ++m_rhsIter;
       }
       else
       {
         m_id = -1;
       }
       return *static_cast<Derived*>(this);
     }

     EIGEN_STRONG_INLINE Scalar value() const { return m_value; }

     EIGEN_STRONG_INLINE int index() const { return m_id; }
     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
     EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }

     EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }

   protected:
     LhsIterator m_lhsIter;
     RhsIterator m_rhsIter;
     const BinaryOp& m_functor;
     Scalar m_value;
     int m_id;
 };

 // sparse - sparse  (product)
 template<typename T, typename Lhs, typename Rhs, typename Derived>
 class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsSparse>
 {
     typedef ei_scalar_product_op<T> BinaryFunc;
     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
     typedef typename CwiseBinaryXpr::Scalar Scalar;
     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
     typedef typename _LhsNested::InnerIterator LhsIterator;
     typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
     typedef typename _RhsNested::InnerIterator RhsIterator;
   public:

     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
       : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
     {
       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
       {
         if (m_lhsIter.index() < m_rhsIter.index())
           ++m_lhsIter;
         else
           ++m_rhsIter;
       }
     }

     EIGEN_STRONG_INLINE Derived& operator++()
     {
       ++m_lhsIter;
       ++m_rhsIter;
       while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
       {
         if (m_lhsIter.index() < m_rhsIter.index())
           ++m_lhsIter;
         else
           ++m_rhsIter;
       }
       return *static_cast<Derived*>(this);
     }

     EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); }

     EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
     EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }

     EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); }

   protected:
     LhsIterator m_lhsIter;
     RhsIterator m_rhsIter;
     const BinaryFunc& m_functor;
 };

 // sparse - dense  (product)
 template<typename T, typename Lhs, typename Rhs, typename Derived>
 class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsDense>
 {
     typedef ei_scalar_product_op<T> BinaryFunc;
     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
     typedef typename CwiseBinaryXpr::Scalar Scalar;
     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
     typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested;
     typedef typename _LhsNested::InnerIterator LhsIterator;
     enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
   public:

     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
       : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
     {}

     EIGEN_STRONG_INLINE Derived& operator++()
     {
       ++m_lhsIter;
       return *static_cast<Derived*>(this);
     }

     EIGEN_STRONG_INLINE Scalar value() const
     { return m_functor(m_lhsIter.value(),
                        m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }

     EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
     EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }

     EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }

   protected:
     const RhsNested m_rhs;
     LhsIterator m_lhsIter;
     const BinaryFunc m_functor;
     const int m_outer;
 };

 // sparse - dense  (product)
 template<typename T, typename Lhs, typename Rhs, typename Derived>
 class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsDense, IsSparse>
 {
     typedef ei_scalar_product_op<T> BinaryFunc;
     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
     typedef typename CwiseBinaryXpr::Scalar Scalar;
     typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
     typedef typename _RhsNested::InnerIterator RhsIterator;
     enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit };
   public:

     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
       : m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer)
     {}

     EIGEN_STRONG_INLINE Derived& operator++()
     {
       ++m_rhsIter;
       return *static_cast<Derived*>(this);
     }

     EIGEN_STRONG_INLINE Scalar value() const
     { return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); }

     EIGEN_STRONG_INLINE int index() const { return m_rhsIter.index(); }
     EIGEN_STRONG_INLINE int row() const { return m_rhsIter.row(); }
     EIGEN_STRONG_INLINE int col() const { return m_rhsIter.col(); }

     EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; }

   protected:
     const CwiseBinaryXpr& m_xpr;
     RhsIterator m_rhsIter;
     const BinaryFunc& m_functor;
     const int m_outer;
 };


 /***************************************************************************
 * Implementation of SparseMatrixBase and SparseCwise functions/operators
 ***************************************************************************/

 template<typename Derived>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>,
                                  Derived, OtherDerived>
 SparseMatrixBase<Derived>::operator-(const SparseMatrixBase<OtherDerived> &other) const
 {
   return SparseCwiseBinaryOp<ei_scalar_difference_op<Scalar>,
                        Derived, OtherDerived>(derived(), other.derived());
 }

 template<typename Derived>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE Derived &
 SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other)
 {
   return *this = derived() - other.derived();
 }

 template<typename Derived>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
 SparseMatrixBase<Derived>::operator+(const SparseMatrixBase<OtherDerived> &other) const
 {
   return SparseCwiseBinaryOp<ei_scalar_sum_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
 }

 template<typename Derived>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE Derived &
 SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other)
 {
   return *this = derived() + other.derived();
 }

 template<typename ExpressionType>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
 SparseCwise<ExpressionType>::operator*(const SparseMatrixBase<OtherDerived> &other) const
 {
   return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
 }

 template<typename ExpressionType>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
 SparseCwise<ExpressionType>::operator*(const MatrixBase<OtherDerived> &other) const
 {
   return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
 }

 // template<typename ExpressionType>
 // template<typename OtherDerived>
 // EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
 // SparseCwise<ExpressionType>::operator/(const SparseMatrixBase<OtherDerived> &other) const
 // {
 //   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
 // }
 //
 // template<typename ExpressionType>
 // template<typename OtherDerived>
 // EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
 // SparseCwise<ExpressionType>::operator/(const MatrixBase<OtherDerived> &other) const
 // {
 //   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
 // }

 template<typename ExpressionType>
 template<typename OtherDerived>
 inline ExpressionType& SparseCwise<ExpressionType>::operator*=(const SparseMatrixBase<OtherDerived> &other)
 {
   return m_matrix.const_cast_derived() = _expression() * other.derived();
 }

 // template<typename ExpressionType>
 // template<typename OtherDerived>
 // inline ExpressionType& SparseCwise<ExpressionType>::operator/=(const SparseMatrixBase<OtherDerived> &other)
 // {
 //   return m_matrix.const_cast_derived() = *this / other;
 // }

 template<typename ExpressionType>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)
 SparseCwise<ExpressionType>::min(const SparseMatrixBase<OtherDerived> &other) const
 {
   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)(_expression(), other.derived());
 }

 template<typename ExpressionType>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)
 SparseCwise<ExpressionType>::max(const SparseMatrixBase<OtherDerived> &other) const
 {
   return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)(_expression(), other.derived());
 }

 // template<typename Derived>
 // template<typename CustomBinaryOp, typename OtherDerived>
 // EIGEN_STRONG_INLINE const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
 // SparseMatrixBase<Derived>::binaryExpr(const SparseMatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const
 // {
 //   return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(derived(), other.derived(), func);
 // }

 #endif // EIGEN_SPARSE_CWISE_BINARY_OP_H
	// This file is part of Eigen, a lightweight C++ template library
	// for linear algebra. Eigen itself is part of the KDE project.
	//
	// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
	//
	// Eigen is free software; you can redistribute it and/or
	// modify it under the terms of the GNU Lesser General Public
	// License as published by the Free Software Foundation; either
	// version 3 of the License, or (at your option) any later version.
	//
	// Alternatively, you can redistribute it and/or
	// modify it under the terms of the GNU General Public License as
	// published by the Free Software Foundation; either version 2 of
	// the License, or (at your option) any later version.
	//
	// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
	// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
	// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
	// GNU General Public License for more details.
	//
	// You should have received a copy of the GNU Lesser General Public
	// License and a copy of the GNU General Public License along with
	// Eigen. If not, see <http://www.gnu.org/licenses/>.

	#ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H
	#define EIGEN_SPARSE_CWISE_BINARY_OP_H

	// Here we have to handle 3 cases:
	// 1 - sparse op dense
	// 2 - dense op sparse
	// 3 - sparse op sparse
	// We also need to implement a 4th iterator for:
	// 4 - dense op dense
	// Finally, we also need to distinguish between the product and other operations :
	// configuration returned mode
	// 1 - sparse op dense product sparse
	// generic dense
	// 2 - dense op sparse product sparse
	// generic dense
	// 3 - sparse op sparse product sparse
	// generic sparse
	// 4 - dense op dense product dense
	// generic dense

	template<typename BinaryOp, typename Lhs, typename Rhs>
	struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
	{
	typedef typename ei_result_of<
	BinaryOp(
	typename Lhs::Scalar,
	typename Rhs::Scalar
	)
	>::type Scalar;
	typedef typename Lhs::Nested LhsNested;
	typedef typename Rhs::Nested RhsNested;
	typedef typename ei_unref<LhsNested>::type _LhsNested;
	typedef typename ei_unref<RhsNested>::type _RhsNested;
	enum {
	LhsCoeffReadCost = _LhsNested::CoeffReadCost,
	RhsCoeffReadCost = _RhsNested::CoeffReadCost,
	LhsFlags = _LhsNested::Flags,
	RhsFlags = _RhsNested::Flags,
	RowsAtCompileTime = Lhs::RowsAtCompileTime,
	ColsAtCompileTime = Lhs::ColsAtCompileTime,
	MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
	MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
	Flags = (int(LhsFlags) \| int(RhsFlags)) & HereditaryBits,
	CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost
	};
	};

	template<typename BinaryOp, typename Lhs, typename Rhs>
	class SparseCwiseBinaryOp : ei_no_assignment_operator,
	public SparseMatrixBase<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
	{
	public:

	class InnerIterator;

	EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseBinaryOp)
	typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested;
	typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested;
	typedef typename ei_unref<LhsNested>::type _LhsNested;
	typedef typename ei_unref<RhsNested>::type _RhsNested;

	EIGEN_STRONG_INLINE SparseCwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
	: m_lhs(lhs), m_rhs(rhs), m_functor(func)
	{
	EIGEN_STATIC_ASSERT((_LhsNested::Flags&RowMajorBit)==(_RhsNested::Flags&RowMajorBit),
	BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER)
	EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret
	? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret)
	: int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)),
	YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
	// require the sizes to match
	EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
	ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
	}

	EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
	EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }

	EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
	EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
	EIGEN_STRONG_INLINE const BinaryOp& functor() const { return m_functor; }

	protected:
	const LhsNested m_lhs;
	const RhsNested m_rhs;
	const BinaryOp m_functor;
	};

	template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived,
	int _LhsStorageMode = int(Lhs::Flags) & SparseBit,
	int _RhsStorageMode = int(Rhs::Flags) & SparseBit>
	class ei_sparse_cwise_binary_op_inner_iterator_selector;

	template<typename BinaryOp, typename Lhs, typename Rhs>
	class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
	: public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator>
	{
	public:
	typedef ei_sparse_cwise_binary_op_inner_iterator_selector<
	BinaryOp,Lhs,Rhs, InnerIterator> Base;

	EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer)
	: Base(binOp,outer)
	{}
	};

	/***************************************************************************
	* Implementation of inner-iterators
	***************************************************************************/

	// template<typename T> struct ei_func_is_conjunction { enum { ret = false }; };
	// template<typename T> struct ei_func_is_conjunction<ei_scalar_product_op<T> > { enum { ret = true }; };

	// TODO generalize the ei_scalar_product_op specialization to all conjunctions if any !

	// sparse - sparse (generic)
	template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
	class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived, IsSparse, IsSparse>
	{
	typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
	typedef typename ei_traits<CwiseBinaryXpr>::Scalar Scalar;
	typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
	typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
	typedef typename _LhsNested::InnerIterator LhsIterator;
	typedef typename _RhsNested::InnerIterator RhsIterator;
	public:

	EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
	: m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
	{
	this->operator++();
	}

	EIGEN_STRONG_INLINE Derived& operator++()
	{
	if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
	{
	m_id = m_lhsIter.index();
	m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
	++m_lhsIter;
	++m_rhsIter;
	}
	else if (m_lhsIter && (!m_rhsIter \|\| (m_lhsIter.index() < m_rhsIter.index())))
	{
	m_id = m_lhsIter.index();
	m_value = m_functor(m_lhsIter.value(), Scalar(0));
	++m_lhsIter;
	}
	else if (m_rhsIter && (!m_lhsIter \|\| (m_lhsIter.index() > m_rhsIter.index())))
	{
	m_id = m_rhsIter.index();
	m_value = m_functor(Scalar(0), m_rhsIter.value());
	++m_rhsIter;
	}
	else
	{
	m_id = -1;
	}
	return static_cast<Derived>(this);
	}

	EIGEN_STRONG_INLINE Scalar value() const { return m_value; }

	EIGEN_STRONG_INLINE int index() const { return m_id; }
	EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
	EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }

	EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }

	protected:
	LhsIterator m_lhsIter;
	RhsIterator m_rhsIter;
	const BinaryOp& m_functor;
	Scalar m_value;
	int m_id;
	};

	// sparse - sparse (product)
	template<typename T, typename Lhs, typename Rhs, typename Derived>
	class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsSparse>
	{
	typedef ei_scalar_product_op<T> BinaryFunc;
	typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
	typedef typename CwiseBinaryXpr::Scalar Scalar;
	typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
	typedef typename _LhsNested::InnerIterator LhsIterator;
	typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
	typedef typename _RhsNested::InnerIterator RhsIterator;
	public:

	EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
	: m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor())
	{
	while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
	{
	if (m_lhsIter.index() < m_rhsIter.index())
	++m_lhsIter;
	else
	++m_rhsIter;
	}
	}

	EIGEN_STRONG_INLINE Derived& operator++()
	{
	++m_lhsIter;
	++m_rhsIter;
	while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index()))
	{
	if (m_lhsIter.index() < m_rhsIter.index())
	++m_lhsIter;
	else
	++m_rhsIter;
	}
	return static_cast<Derived>(this);
	}

	EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); }

	EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
	EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
	EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }

	EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); }

	protected:
	LhsIterator m_lhsIter;
	RhsIterator m_rhsIter;
	const BinaryFunc& m_functor;
	};

	// sparse - dense (product)
	template<typename T, typename Lhs, typename Rhs, typename Derived>
	class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsDense>
	{
	typedef ei_scalar_product_op<T> BinaryFunc;
	typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
	typedef typename CwiseBinaryXpr::Scalar Scalar;
	typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
	typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested;
	typedef typename _LhsNested::InnerIterator LhsIterator;
	enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
	public:

	EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
	: m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
	{}

	EIGEN_STRONG_INLINE Derived& operator++()
	{
	++m_lhsIter;
	return static_cast<Derived>(this);
	}

	EIGEN_STRONG_INLINE Scalar value() const
	{ return m_functor(m_lhsIter.value(),
	m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }

	EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
	EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
	EIGEN_STRONG_INLINE int col() const { return m_lhsIter.col(); }

	EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }

	protected:
	const RhsNested m_rhs;
	LhsIterator m_lhsIter;
	const BinaryFunc m_functor;
	const int m_outer;
	};

	// sparse - dense (product)
	template<typename T, typename Lhs, typename Rhs, typename Derived>
	class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsDense, IsSparse>
	{
	typedef ei_scalar_product_op<T> BinaryFunc;
	typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
	typedef typename CwiseBinaryXpr::Scalar Scalar;
	typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
	typedef typename _RhsNested::InnerIterator RhsIterator;
	enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit };
	public:

	EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
	: m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer)
	{}

	EIGEN_STRONG_INLINE Derived& operator++()
	{
	++m_rhsIter;
	return static_cast<Derived>(this);
	}

	EIGEN_STRONG_INLINE Scalar value() const
	{ return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); }

	EIGEN_STRONG_INLINE int index() const { return m_rhsIter.index(); }
	EIGEN_STRONG_INLINE int row() const { return m_rhsIter.row(); }
	EIGEN_STRONG_INLINE int col() const { return m_rhsIter.col(); }

	EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; }

	protected:
	const CwiseBinaryXpr& m_xpr;
	RhsIterator m_rhsIter;
	const BinaryFunc& m_functor;
	const int m_outer;
	};


	/***************************************************************************
	* Implementation of SparseMatrixBase and SparseCwise functions/operators
	***************************************************************************/

	template<typename Derived>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>,
	Derived, OtherDerived>
	SparseMatrixBase<Derived>::operator-(const SparseMatrixBase<OtherDerived> &other) const
	{
	return SparseCwiseBinaryOp<ei_scalar_difference_op<Scalar>,
	Derived, OtherDerived>(derived(), other.derived());
	}

	template<typename Derived>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE Derived &
	SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other)
	{
	return *this = derived() - other.derived();
	}

	template<typename Derived>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
	SparseMatrixBase<Derived>::operator+(const SparseMatrixBase<OtherDerived> &other) const
	{
	return SparseCwiseBinaryOp<ei_scalar_sum_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
	}

	template<typename Derived>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE Derived &
	SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other)
	{
	return *this = derived() + other.derived();
	}

	template<typename ExpressionType>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
	SparseCwise<ExpressionType>::operator*(const SparseMatrixBase<OtherDerived> &other) const
	{
	return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
	}

	template<typename ExpressionType>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
	SparseCwise<ExpressionType>::operator*(const MatrixBase<OtherDerived> &other) const
	{
	return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
	}

	// template<typename ExpressionType>
	// template<typename OtherDerived>
	// EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
	// SparseCwise<ExpressionType>::operator/(const SparseMatrixBase<OtherDerived> &other) const
	// {
	// return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
	// }
	//
	// template<typename ExpressionType>
	// template<typename OtherDerived>
	// EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
	// SparseCwise<ExpressionType>::operator/(const MatrixBase<OtherDerived> &other) const
	// {
	// return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
	// }

	template<typename ExpressionType>
	template<typename OtherDerived>
	inline ExpressionType& SparseCwise<ExpressionType>::operator*=(const SparseMatrixBase<OtherDerived> &other)
	{
	return m_matrix.const_cast_derived() = _expression() * other.derived();
	}

	// template<typename ExpressionType>
	// template<typename OtherDerived>
	// inline ExpressionType& SparseCwise<ExpressionType>::operator/=(const SparseMatrixBase<OtherDerived> &other)
	// {
	// return m_matrix.const_cast_derived() = *this / other;
	// }

	template<typename ExpressionType>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)
	SparseCwise<ExpressionType>::min(const SparseMatrixBase<OtherDerived> &other) const
	{
	return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)(_expression(), other.derived());
	}

	template<typename ExpressionType>
	template<typename OtherDerived>
	EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)
	SparseCwise<ExpressionType>::max(const SparseMatrixBase<OtherDerived> &other) const
	{
	return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)(_expression(), other.derived());
	}

	// template<typename Derived>
	// template<typename CustomBinaryOp, typename OtherDerived>
	// EIGEN_STRONG_INLINE const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
	// SparseMatrixBase<Derived>::binaryExpr(const SparseMatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const
	// {
	// return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(derived(), other.derived(), func);
	// }

	#endif // EIGEN_SPARSE_CWISE_BINARY_OP_H