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 (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H 11 #define EIGEN_SPARSE_CWISE_BINARY_OP_H 12 13 namespace Eigen { 14 15 // Here we have to handle 3 cases: 16 // 1 - sparse op dense 17 // 2 - dense op sparse 18 // 3 - sparse op sparse 19 // We also need to implement a 4th iterator for: 20 // 4 - dense op dense 21 // Finally, we also need to distinguish between the product and other operations : 22 // configuration returned mode 23 // 1 - sparse op dense product sparse 24 // generic dense 25 // 2 - dense op sparse product sparse 26 // generic dense 27 // 3 - sparse op sparse product sparse 28 // generic sparse 29 // 4 - dense op dense product dense 30 // generic dense 31 32 namespace internal { 33 34 template<> struct promote_storage_type<Dense,Sparse> 35 { typedef Sparse ret; }; 36 37 template<> struct promote_storage_type<Sparse,Dense> 38 { typedef Sparse ret; }; 39 40 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived, 41 typename _LhsStorageMode = typename traits<Lhs>::StorageKind, 42 typename _RhsStorageMode = typename traits<Rhs>::StorageKind> 43 class sparse_cwise_binary_op_inner_iterator_selector; 44 45 } // end namespace internal 46 47 template<typename BinaryOp, typename Lhs, typename Rhs> 48 class CwiseBinaryOpImpl<BinaryOp, Lhs, Rhs, Sparse> 49 : public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> > 50 { 51 public: 52 class InnerIterator; 53 class ReverseInnerIterator; 54 typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> Derived; 55 EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) 56 CwiseBinaryOpImpl() 57 { 58 typedef typename internal::traits<Lhs>::StorageKind LhsStorageKind; 59 typedef typename internal::traits<Rhs>::StorageKind RhsStorageKind; 60 EIGEN_STATIC_ASSERT(( 61 (!internal::is_same<LhsStorageKind,RhsStorageKind>::value) 62 || ((Lhs::Flags&RowMajorBit) == (Rhs::Flags&RowMajorBit))), 63 THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH); 64 } 65 }; 66 67 template<typename BinaryOp, typename Lhs, typename Rhs> 68 class CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator 69 : public internal::sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs,typename CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator> 70 { 71 public: 72 typedef typename Lhs::Index Index; 73 typedef internal::sparse_cwise_binary_op_inner_iterator_selector< 74 BinaryOp,Lhs,Rhs, InnerIterator> Base; 75 76 // 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 91 namespace internal { 92 93 // sparse - sparse (generic) 94 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived> 95 class 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) 159 template<typename T, typename Lhs, typename Rhs, typename Derived> 160 class 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) 213 template<typename T, typename Lhs, typename Rhs, typename Derived> 214 class 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) 254 template<typename T, typename Lhs, typename Rhs, typename Derived> 255 class 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 299 template<typename Derived> 300 template<typename OtherDerived> 301 EIGEN_STRONG_INLINE Derived & 302 SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other) 303 { 304 return derived() = derived() - other.derived(); 305 } 306 307 template<typename Derived> 308 template<typename OtherDerived> 309 EIGEN_STRONG_INLINE Derived & 310 SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other) 311 { 312 return derived() = derived() + other.derived(); 313 } 314 315 template<typename Derived> 316 template<typename OtherDerived> 317 EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE 318 SparseMatrixBase<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
