1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_SELFCWISEBINARYOP_H
11#define EIGEN_SELFCWISEBINARYOP_H
12
13namespace Eigen {
14
15/** \class SelfCwiseBinaryOp
16  * \ingroup Core_Module
17  *
18  * \internal
19  *
20  * \brief Internal helper class for optimizing operators like +=, -=
21  *
22  * This is a pseudo expression class re-implementing the copyCoeff/copyPacket
23  * method to directly performs a +=/-= operations in an optimal way. In particular,
24  * this allows to make sure that the input/output data are loaded only once using
25  * aligned packet loads.
26  *
27  * \sa class SwapWrapper for a similar trick.
28  */
29
30namespace internal {
31template<typename BinaryOp, typename Lhs, typename Rhs>
32struct traits<SelfCwiseBinaryOp<BinaryOp,Lhs,Rhs> >
33  : traits<CwiseBinaryOp<BinaryOp,Lhs,Rhs> >
34{
35  enum {
36    // Note that it is still a good idea to preserve the DirectAccessBit
37    // so that assign can correctly align the data.
38    Flags = traits<CwiseBinaryOp<BinaryOp,Lhs,Rhs> >::Flags | (Lhs::Flags&DirectAccessBit) | (Lhs::Flags&LvalueBit),
39    OuterStrideAtCompileTime = Lhs::OuterStrideAtCompileTime,
40    InnerStrideAtCompileTime = Lhs::InnerStrideAtCompileTime
41  };
42};
43}
44
45template<typename BinaryOp, typename Lhs, typename Rhs> class SelfCwiseBinaryOp
46  : public internal::dense_xpr_base< SelfCwiseBinaryOp<BinaryOp, Lhs, Rhs> >::type
47{
48  public:
49
50    typedef typename internal::dense_xpr_base<SelfCwiseBinaryOp>::type Base;
51    EIGEN_DENSE_PUBLIC_INTERFACE(SelfCwiseBinaryOp)
52
53    typedef typename internal::packet_traits<Scalar>::type Packet;
54
55    inline SelfCwiseBinaryOp(Lhs& xpr, const BinaryOp& func = BinaryOp()) : m_matrix(xpr), m_functor(func) {}
56
57    inline Index rows() const { return m_matrix.rows(); }
58    inline Index cols() const { return m_matrix.cols(); }
59    inline Index outerStride() const { return m_matrix.outerStride(); }
60    inline Index innerStride() const { return m_matrix.innerStride(); }
61    inline const Scalar* data() const { return m_matrix.data(); }
62
63    // note that this function is needed by assign to correctly align loads/stores
64    // TODO make Assign use .data()
65    inline Scalar& coeffRef(Index row, Index col)
66    {
67      EIGEN_STATIC_ASSERT_LVALUE(Lhs)
68      return m_matrix.const_cast_derived().coeffRef(row, col);
69    }
70    inline const Scalar& coeffRef(Index row, Index col) const
71    {
72      return m_matrix.coeffRef(row, col);
73    }
74
75    // note that this function is needed by assign to correctly align loads/stores
76    // TODO make Assign use .data()
77    inline Scalar& coeffRef(Index index)
78    {
79      EIGEN_STATIC_ASSERT_LVALUE(Lhs)
80      return m_matrix.const_cast_derived().coeffRef(index);
81    }
82    inline const Scalar& coeffRef(Index index) const
83    {
84      return m_matrix.const_cast_derived().coeffRef(index);
85    }
86
87    template<typename OtherDerived>
88    void copyCoeff(Index row, Index col, const DenseBase<OtherDerived>& other)
89    {
90      OtherDerived& _other = other.const_cast_derived();
91      eigen_internal_assert(row >= 0 && row < rows()
92                         && col >= 0 && col < cols());
93      Scalar& tmp = m_matrix.coeffRef(row,col);
94      tmp = m_functor(tmp, _other.coeff(row,col));
95    }
96
97    template<typename OtherDerived>
98    void copyCoeff(Index index, const DenseBase<OtherDerived>& other)
99    {
100      OtherDerived& _other = other.const_cast_derived();
101      eigen_internal_assert(index >= 0 && index < m_matrix.size());
102      Scalar& tmp = m_matrix.coeffRef(index);
103      tmp = m_functor(tmp, _other.coeff(index));
104    }
105
106    template<typename OtherDerived, int StoreMode, int LoadMode>
107    void copyPacket(Index row, Index col, const DenseBase<OtherDerived>& other)
108    {
109      OtherDerived& _other = other.const_cast_derived();
110      eigen_internal_assert(row >= 0 && row < rows()
111                        && col >= 0 && col < cols());
112      m_matrix.template writePacket<StoreMode>(row, col,
113        m_functor.packetOp(m_matrix.template packet<StoreMode>(row, col),_other.template packet<LoadMode>(row, col)) );
114    }
115
116    template<typename OtherDerived, int StoreMode, int LoadMode>
117    void copyPacket(Index index, const DenseBase<OtherDerived>& other)
118    {
119      OtherDerived& _other = other.const_cast_derived();
120      eigen_internal_assert(index >= 0 && index < m_matrix.size());
121      m_matrix.template writePacket<StoreMode>(index,
122        m_functor.packetOp(m_matrix.template packet<StoreMode>(index),_other.template packet<LoadMode>(index)) );
123    }
124
125    // reimplement lazyAssign to handle complex *= real
126    // see CwiseBinaryOp ctor for details
127    template<typename RhsDerived>
128    EIGEN_STRONG_INLINE SelfCwiseBinaryOp& lazyAssign(const DenseBase<RhsDerived>& rhs)
129    {
130      EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs,RhsDerived)
131      EIGEN_CHECK_BINARY_COMPATIBILIY(BinaryOp,typename Lhs::Scalar,typename RhsDerived::Scalar);
132
133    #ifdef EIGEN_DEBUG_ASSIGN
134      internal::assign_traits<SelfCwiseBinaryOp, RhsDerived>::debug();
135    #endif
136      eigen_assert(rows() == rhs.rows() && cols() == rhs.cols());
137      internal::assign_impl<SelfCwiseBinaryOp, RhsDerived>::run(*this,rhs.derived());
138    #ifndef EIGEN_NO_DEBUG
139      this->checkTransposeAliasing(rhs.derived());
140    #endif
141      return *this;
142    }
143
144    // overloaded to honor evaluation of special matrices
145    // maybe another solution would be to not use SelfCwiseBinaryOp
146    // at first...
147    SelfCwiseBinaryOp& operator=(const Rhs& _rhs)
148    {
149      typename internal::nested<Rhs>::type rhs(_rhs);
150      return Base::operator=(rhs);
151    }
152
153    Lhs& expression() const
154    {
155      return m_matrix;
156    }
157
158    const BinaryOp& functor() const
159    {
160      return m_functor;
161    }
162
163  protected:
164    Lhs& m_matrix;
165    const BinaryOp& m_functor;
166
167  private:
168    SelfCwiseBinaryOp& operator=(const SelfCwiseBinaryOp&);
169};
170
171template<typename Derived>
172inline Derived& DenseBase<Derived>::operator*=(const Scalar& other)
173{
174  typedef typename Derived::PlainObject PlainObject;
175  SelfCwiseBinaryOp<internal::scalar_product_op<Scalar>, Derived, typename PlainObject::ConstantReturnType> tmp(derived());
176  tmp = PlainObject::Constant(rows(),cols(),other);
177  return derived();
178}
179
180template<typename Derived>
181inline Derived& DenseBase<Derived>::operator/=(const Scalar& other)
182{
183  typedef typename internal::conditional<NumTraits<Scalar>::IsInteger,
184                                        internal::scalar_quotient_op<Scalar>,
185                                        internal::scalar_product_op<Scalar> >::type BinOp;
186  typedef typename Derived::PlainObject PlainObject;
187  SelfCwiseBinaryOp<BinOp, Derived, typename PlainObject::ConstantReturnType> tmp(derived());
188  Scalar actual_other;
189  if(NumTraits<Scalar>::IsInteger)  actual_other = other;
190  else                              actual_other = Scalar(1)/other;
191  tmp = PlainObject::Constant(rows(),cols(), actual_other);
192  return derived();
193}
194
195} // end namespace Eigen
196
197#endif // EIGEN_SELFCWISEBINARYOP_H
198