Home | History | Annotate | Download | only in util 
      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 (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_BLASUTIL_H
     11 #define EIGEN_BLASUTIL_H
     12 
     13 // This file contains many lightweight helper classes used to
     14 // implement and control fast level 2 and level 3 BLAS-like routines.
     15 
     16 namespace Eigen {
     17 
     18 namespace internal {
     19 
     20 // forward declarations
     21 template<typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs=false, bool ConjugateRhs=false>
     22 struct gebp_kernel;
     23 
     24 template<typename Scalar, typename Index, typename DataMapper, int nr, int StorageOrder, bool Conjugate = false, bool PanelMode=false>
     25 struct gemm_pack_rhs;
     26 
     27 template<typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, int StorageOrder, bool Conjugate = false, bool PanelMode = false>
     28 struct gemm_pack_lhs;
     29 
     30 template<
     31   typename Index,
     32   typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
     33   typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs,
     34   int ResStorageOrder>
     35 struct general_matrix_matrix_product;
     36 
     37 template<typename Index,
     38          typename LhsScalar, typename LhsMapper, int LhsStorageOrder, bool ConjugateLhs,
     39          typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version=Specialized>
     40 struct general_matrix_vector_product;
     41 
     42 
     43 template<bool Conjugate> struct conj_if;
     44 
     45 template<> struct conj_if<true> {
     46   template<typename T>
     47   inline T operator()(const T& x) const { return numext::conj(x); }
     48   template<typename T>
     49   inline T pconj(const T& x) const { return internal::pconj(x); }
     50 };
     51 
     52 template<> struct conj_if<false> {
     53   template<typename T>
     54   inline const T& operator()(const T& x) const { return x; }
     55   template<typename T>
     56   inline const T& pconj(const T& x) const { return x; }
     57 };
     58 
     59 // Generic implementation for custom complex types.
     60 template<typename LhsScalar, typename RhsScalar, bool ConjLhs, bool ConjRhs>
     61 struct conj_helper
     62 {
     63   typedef typename ScalarBinaryOpTraits<LhsScalar,RhsScalar>::ReturnType Scalar;
     64 
     65   EIGEN_STRONG_INLINE Scalar pmadd(const LhsScalar& x, const RhsScalar& y, const Scalar& c) const
     66   { return padd(c, pmul(x,y)); }
     67 
     68   EIGEN_STRONG_INLINE Scalar pmul(const LhsScalar& x, const RhsScalar& y) const
     69   { return conj_if<ConjLhs>()(x) *  conj_if<ConjRhs>()(y); }
     70 };
     71 
     72 template<typename Scalar> struct conj_helper<Scalar,Scalar,false,false>
     73 {
     74   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const { return internal::pmadd(x,y,c); }
     75   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const { return internal::pmul(x,y); }
     76 };
     77 
     78 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, false,true>
     79 {
     80   typedef std::complex<RealScalar> Scalar;
     81   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
     82   { return c + pmul(x,y); }
     83 
     84   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
     85   { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::imag(x)*numext::real(y) - numext::real(x)*numext::imag(y)); }
     86 };
     87 
     88 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,false>
     89 {
     90   typedef std::complex<RealScalar> Scalar;
     91   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
     92   { return c + pmul(x,y); }
     93 
     94   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
     95   { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); }
     96 };
     97 
     98 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,true>
     99 {
    100   typedef std::complex<RealScalar> Scalar;
    101   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
    102   { return c + pmul(x,y); }
    103 
    104   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
    105   { return Scalar(numext::real(x)*numext::real(y) - numext::imag(x)*numext::imag(y), - numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); }
    106 };
    107 
    108 template<typename RealScalar,bool Conj> struct conj_helper<std::complex<RealScalar>, RealScalar, Conj,false>
    109 {
    110   typedef std::complex<RealScalar> Scalar;
    111   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const
    112   { return padd(c, pmul(x,y)); }
    113   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const
    114   { return conj_if<Conj>()(x)*y; }
    115 };
    116 
    117 template<typename RealScalar,bool Conj> struct conj_helper<RealScalar, std::complex<RealScalar>, false,Conj>
    118 {
    119   typedef std::complex<RealScalar> Scalar;
    120   EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const
    121   { return padd(c, pmul(x,y)); }
    122   EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const
    123   { return x*conj_if<Conj>()(y); }
    124 };
    125 
    126 template<typename From,typename To> struct get_factor {
    127   EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE To run(const From& x) { return To(x); }
    128 };
    129 
    130 template<typename Scalar> struct get_factor<Scalar,typename NumTraits<Scalar>::Real> {
    131   EIGEN_DEVICE_FUNC
    132   static EIGEN_STRONG_INLINE typename NumTraits<Scalar>::Real run(const Scalar& x) { return numext::real(x); }
    133 };
    134 
    135 
    136 template<typename Scalar, typename Index>
    137 class BlasVectorMapper {
    138   public:
    139   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE BlasVectorMapper(Scalar *data) : m_data(data) {}
    140 
    141   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar operator()(Index i) const {
    142     return m_data[i];
    143   }
    144   template <typename Packet, int AlignmentType>
    145   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet load(Index i) const {
    146     return ploadt<Packet, AlignmentType>(m_data + i);
    147   }
    148 
    149   template <typename Packet>
    150   EIGEN_DEVICE_FUNC bool aligned(Index i) const {
    151     return (UIntPtr(m_data+i)%sizeof(Packet))==0;
    152   }
    153 
    154   protected:
    155   Scalar* m_data;
    156 };
    157 
    158 template<typename Scalar, typename Index, int AlignmentType>
    159 class BlasLinearMapper {
    160   public:
    161   typedef typename packet_traits<Scalar>::type Packet;
    162   typedef typename packet_traits<Scalar>::half HalfPacket;
    163 
    164   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE BlasLinearMapper(Scalar *data) : m_data(data) {}
    165 
    166   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void prefetch(int i) const {
    167     internal::prefetch(&operator()(i));
    168   }
    169 
    170   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar& operator()(Index i) const {
    171     return m_data[i];
    172   }
    173 
    174   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet loadPacket(Index i) const {
    175     return ploadt<Packet, AlignmentType>(m_data + i);
    176   }
    177 
    178   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE HalfPacket loadHalfPacket(Index i) const {
    179     return ploadt<HalfPacket, AlignmentType>(m_data + i);
    180   }
    181 
    182   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void storePacket(Index i, const Packet &p) const {
    183     pstoret<Scalar, Packet, AlignmentType>(m_data + i, p);
    184   }
    185 
    186   protected:
    187   Scalar *m_data;
    188 };
    189 
    190 // Lightweight helper class to access matrix coefficients.
    191 template<typename Scalar, typename Index, int StorageOrder, int AlignmentType = Unaligned>
    192 class blas_data_mapper {
    193   public:
    194   typedef typename packet_traits<Scalar>::type Packet;
    195   typedef typename packet_traits<Scalar>::half HalfPacket;
    196 
    197   typedef BlasLinearMapper<Scalar, Index, AlignmentType> LinearMapper;
    198   typedef BlasVectorMapper<Scalar, Index> VectorMapper;
    199 
    200   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE blas_data_mapper(Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
    201 
    202   EIGEN_DEVICE_FUNC  EIGEN_ALWAYS_INLINE blas_data_mapper<Scalar, Index, StorageOrder, AlignmentType>
    203   getSubMapper(Index i, Index j) const {
    204     return blas_data_mapper<Scalar, Index, StorageOrder, AlignmentType>(&operator()(i, j), m_stride);
    205   }
    206 
    207   EIGEN_DEVICE_FUNC  EIGEN_ALWAYS_INLINE LinearMapper getLinearMapper(Index i, Index j) const {
    208     return LinearMapper(&operator()(i, j));
    209   }
    210 
    211   EIGEN_DEVICE_FUNC  EIGEN_ALWAYS_INLINE VectorMapper getVectorMapper(Index i, Index j) const {
    212     return VectorMapper(&operator()(i, j));
    213   }
    214 
    215 
    216   EIGEN_DEVICE_FUNC
    217   EIGEN_ALWAYS_INLINE Scalar& operator()(Index i, Index j) const {
    218     return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride];
    219   }
    220 
    221   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet loadPacket(Index i, Index j) const {
    222     return ploadt<Packet, AlignmentType>(&operator()(i, j));
    223   }
    224 
    225   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE HalfPacket loadHalfPacket(Index i, Index j) const {
    226     return ploadt<HalfPacket, AlignmentType>(&operator()(i, j));
    227   }
    228 
    229   template<typename SubPacket>
    230   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void scatterPacket(Index i, Index j, const SubPacket &p) const {
    231     pscatter<Scalar, SubPacket>(&operator()(i, j), p, m_stride);
    232   }
    233 
    234   template<typename SubPacket>
    235   EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE SubPacket gatherPacket(Index i, Index j) const {
    236     return pgather<Scalar, SubPacket>(&operator()(i, j), m_stride);
    237   }
    238 
    239   EIGEN_DEVICE_FUNC const Index stride() const { return m_stride; }
    240   EIGEN_DEVICE_FUNC const Scalar* data() const { return m_data; }
    241 
    242   EIGEN_DEVICE_FUNC Index firstAligned(Index size) const {
    243     if (UIntPtr(m_data)%sizeof(Scalar)) {
    244       return -1;
    245     }
    246     return internal::first_default_aligned(m_data, size);
    247   }
    248 
    249   protected:
    250   Scalar* EIGEN_RESTRICT m_data;
    251   const Index m_stride;
    252 };
    253 
    254 // lightweight helper class to access matrix coefficients (const version)
    255 template<typename Scalar, typename Index, int StorageOrder>
    256 class const_blas_data_mapper : public blas_data_mapper<const Scalar, Index, StorageOrder> {
    257   public:
    258   EIGEN_ALWAYS_INLINE const_blas_data_mapper(const Scalar *data, Index stride) : blas_data_mapper<const Scalar, Index, StorageOrder>(data, stride) {}
    259 
    260   EIGEN_ALWAYS_INLINE const_blas_data_mapper<Scalar, Index, StorageOrder> getSubMapper(Index i, Index j) const {
    261     return const_blas_data_mapper<Scalar, Index, StorageOrder>(&(this->operator()(i, j)), this->m_stride);
    262   }
    263 };
    264 
    265 
    266 /* Helper class to analyze the factors of a Product expression.
    267  * In particular it allows to pop out operator-, scalar multiples,
    268  * and conjugate */
    269 template<typename XprType> struct blas_traits
    270 {
    271   typedef typename traits<XprType>::Scalar Scalar;
    272   typedef const XprType& ExtractType;
    273   typedef XprType _ExtractType;
    274   enum {
    275     IsComplex = NumTraits<Scalar>::IsComplex,
    276     IsTransposed = false,
    277     NeedToConjugate = false,
    278     HasUsableDirectAccess = (    (int(XprType::Flags)&DirectAccessBit)
    279                               && (   bool(XprType::IsVectorAtCompileTime)
    280                                   || int(inner_stride_at_compile_time<XprType>::ret) == 1)
    281                              ) ?  1 : 0
    282   };
    283   typedef typename conditional<bool(HasUsableDirectAccess),
    284     ExtractType,
    285     typename _ExtractType::PlainObject
    286     >::type DirectLinearAccessType;
    287   static inline ExtractType extract(const XprType& x) { return x; }
    288   static inline const Scalar extractScalarFactor(const XprType&) { return Scalar(1); }
    289 };
    290 
    291 // pop conjugate
    292 template<typename Scalar, typename NestedXpr>
    293 struct blas_traits<CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> >
    294  : blas_traits<NestedXpr>
    295 {
    296   typedef blas_traits<NestedXpr> Base;
    297   typedef CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> XprType;
    298   typedef typename Base::ExtractType ExtractType;
    299 
    300   enum {
    301     IsComplex = NumTraits<Scalar>::IsComplex,
    302     NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex
    303   };
    304   static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
    305   static inline Scalar extractScalarFactor(const XprType& x) { return conj(Base::extractScalarFactor(x.nestedExpression())); }
    306 };
    307 
    308 // pop scalar multiple
    309 template<typename Scalar, typename NestedXpr, typename Plain>
    310 struct blas_traits<CwiseBinaryOp<scalar_product_op<Scalar>, const CwiseNullaryOp<scalar_constant_op<Scalar>,Plain>, NestedXpr> >
    311  : blas_traits<NestedXpr>
    312 {
    313   typedef blas_traits<NestedXpr> Base;
    314   typedef CwiseBinaryOp<scalar_product_op<Scalar>, const CwiseNullaryOp<scalar_constant_op<Scalar>,Plain>, NestedXpr> XprType;
    315   typedef typename Base::ExtractType ExtractType;
    316   static inline ExtractType extract(const XprType& x) { return Base::extract(x.rhs()); }
    317   static inline Scalar extractScalarFactor(const XprType& x)
    318   { return x.lhs().functor().m_other * Base::extractScalarFactor(x.rhs()); }
    319 };
    320 template<typename Scalar, typename NestedXpr, typename Plain>
    321 struct blas_traits<CwiseBinaryOp<scalar_product_op<Scalar>, NestedXpr, const CwiseNullaryOp<scalar_constant_op<Scalar>,Plain> > >
    322  : blas_traits<NestedXpr>
    323 {
    324   typedef blas_traits<NestedXpr> Base;
    325   typedef CwiseBinaryOp<scalar_product_op<Scalar>, NestedXpr, const CwiseNullaryOp<scalar_constant_op<Scalar>,Plain> > XprType;
    326   typedef typename Base::ExtractType ExtractType;
    327   static inline ExtractType extract(const XprType& x) { return Base::extract(x.lhs()); }
    328   static inline Scalar extractScalarFactor(const XprType& x)
    329   { return Base::extractScalarFactor(x.lhs()) * x.rhs().functor().m_other; }
    330 };
    331 template<typename Scalar, typename Plain1, typename Plain2>
    332 struct blas_traits<CwiseBinaryOp<scalar_product_op<Scalar>, const CwiseNullaryOp<scalar_constant_op<Scalar>,Plain1>,
    333                                                             const CwiseNullaryOp<scalar_constant_op<Scalar>,Plain2> > >
    334  : blas_traits<CwiseNullaryOp<scalar_constant_op<Scalar>,Plain1> >
    335 {};
    336 
    337 // pop opposite
    338 template<typename Scalar, typename NestedXpr>
    339 struct blas_traits<CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> >
    340  : blas_traits<NestedXpr>
    341 {
    342   typedef blas_traits<NestedXpr> Base;
    343   typedef CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> XprType;
    344   typedef typename Base::ExtractType ExtractType;
    345   static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
    346   static inline Scalar extractScalarFactor(const XprType& x)
    347   { return - Base::extractScalarFactor(x.nestedExpression()); }
    348 };
    349 
    350 // pop/push transpose
    351 template<typename NestedXpr>
    352 struct blas_traits<Transpose<NestedXpr> >
    353  : blas_traits<NestedXpr>
    354 {
    355   typedef typename NestedXpr::Scalar Scalar;
    356   typedef blas_traits<NestedXpr> Base;
    357   typedef Transpose<NestedXpr> XprType;
    358   typedef Transpose<const typename Base::_ExtractType>  ExtractType; // const to get rid of a compile error; anyway blas traits are only used on the RHS
    359   typedef Transpose<const typename Base::_ExtractType> _ExtractType;
    360   typedef typename conditional<bool(Base::HasUsableDirectAccess),
    361     ExtractType,
    362     typename ExtractType::PlainObject
    363     >::type DirectLinearAccessType;
    364   enum {
    365     IsTransposed = Base::IsTransposed ? 0 : 1
    366   };
    367   static inline ExtractType extract(const XprType& x) { return ExtractType(Base::extract(x.nestedExpression())); }
    368   static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); }
    369 };
    370 
    371 template<typename T>
    372 struct blas_traits<const T>
    373      : blas_traits<T>
    374 {};
    375 
    376 template<typename T, bool HasUsableDirectAccess=blas_traits<T>::HasUsableDirectAccess>
    377 struct extract_data_selector {
    378   static const typename T::Scalar* run(const T& m)
    379   {
    380     return blas_traits<T>::extract(m).data();
    381   }
    382 };
    383 
    384 template<typename T>
    385 struct extract_data_selector<T,false> {
    386   static typename T::Scalar* run(const T&) { return 0; }
    387 };
    388 
    389 template<typename T> const typename T::Scalar* extract_data(const T& m)
    390 {
    391   return extract_data_selector<T>::run(m);
    392 }
    393 
    394 } // end namespace internal
    395 
    396 } // end namespace Eigen
    397 
    398 #endif // EIGEN_BLASUTIL_H
    399