Michael Ernst Peter / Eigen

src/Core/util/BlasUtil.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 5 //
ykuroda 0:13a5d365ba16 6 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 7 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 9
ykuroda 0:13a5d365ba16 10 #ifndef EIGEN_BLASUTIL_H
ykuroda 0:13a5d365ba16 11 #define EIGEN_BLASUTIL_H
ykuroda 0:13a5d365ba16 12
ykuroda 0:13a5d365ba16 13 // This file contains many lightweight helper classes used to
ykuroda 0:13a5d365ba16 14 // implement and control fast level 2 and level 3 BLAS-like routines.
ykuroda 0:13a5d365ba16 15
ykuroda 0:13a5d365ba16 16 namespace Eigen {
ykuroda 0:13a5d365ba16 17
ykuroda 0:13a5d365ba16 18 namespace internal {
ykuroda 0:13a5d365ba16 19
ykuroda 0:13a5d365ba16 20 // forward declarations
ykuroda 0:13a5d365ba16 21 template<typename LhsScalar, typename RhsScalar, typename Index, int mr, int nr, bool ConjugateLhs=false, bool ConjugateRhs=false>
ykuroda 0:13a5d365ba16 22 struct gebp_kernel;
ykuroda 0:13a5d365ba16 23
ykuroda 0:13a5d365ba16 24 template<typename Scalar, typename Index, int nr, int StorageOrder, bool Conjugate = false, bool PanelMode=false>
ykuroda 0:13a5d365ba16 25 struct gemm_pack_rhs;
ykuroda 0:13a5d365ba16 26
ykuroda 0:13a5d365ba16 27 template<typename Scalar, typename Index, int Pack1, int Pack2, int StorageOrder, bool Conjugate = false, bool PanelMode = false>
ykuroda 0:13a5d365ba16 28 struct gemm_pack_lhs;
ykuroda 0:13a5d365ba16 29
ykuroda 0:13a5d365ba16 30 template<
ykuroda 0:13a5d365ba16 31 typename Index,
ykuroda 0:13a5d365ba16 32 typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
ykuroda 0:13a5d365ba16 33 typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs,
ykuroda 0:13a5d365ba16 34 int ResStorageOrder>
ykuroda 0:13a5d365ba16 35 struct general_matrix_matrix_product;
ykuroda 0:13a5d365ba16 36
ykuroda 0:13a5d365ba16 37 template<typename Index, typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version=Specialized>
ykuroda 0:13a5d365ba16 38 struct general_matrix_vector_product;
ykuroda 0:13a5d365ba16 39
ykuroda 0:13a5d365ba16 40
ykuroda 0:13a5d365ba16 41 template<bool Conjugate> struct conj_if;
ykuroda 0:13a5d365ba16 42
ykuroda 0:13a5d365ba16 43 template<> struct conj_if<true> {
ykuroda 0:13a5d365ba16 44 template<typename T>
ykuroda 0:13a5d365ba16 45 inline T operator()(const T& x) { return numext::conj(x); }
ykuroda 0:13a5d365ba16 46 template<typename T>
ykuroda 0:13a5d365ba16 47 inline T pconj(const T& x) { return internal::pconj(x); }
ykuroda 0:13a5d365ba16 48 };
ykuroda 0:13a5d365ba16 49
ykuroda 0:13a5d365ba16 50 template<> struct conj_if<false> {
ykuroda 0:13a5d365ba16 51 template<typename T>
ykuroda 0:13a5d365ba16 52 inline const T& operator()(const T& x) { return x; }
ykuroda 0:13a5d365ba16 53 template<typename T>
ykuroda 0:13a5d365ba16 54 inline const T& pconj(const T& x) { return x; }
ykuroda 0:13a5d365ba16 55 };
ykuroda 0:13a5d365ba16 56
ykuroda 0:13a5d365ba16 57 template<typename Scalar> struct conj_helper<Scalar,Scalar,false,false>
ykuroda 0:13a5d365ba16 58 {
ykuroda 0:13a5d365ba16 59 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const { return internal::pmadd(x,y,c); }
ykuroda 0:13a5d365ba16 60 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const { return internal::pmul(x,y); }
ykuroda 0:13a5d365ba16 61 };
ykuroda 0:13a5d365ba16 62
ykuroda 0:13a5d365ba16 63 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, false,true>
ykuroda 0:13a5d365ba16 64 {
ykuroda 0:13a5d365ba16 65 typedef std::complex<RealScalar> Scalar;
ykuroda 0:13a5d365ba16 66 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
ykuroda 0:13a5d365ba16 67 { return c + pmul(x,y); }
ykuroda 0:13a5d365ba16 68
ykuroda 0:13a5d365ba16 69 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
ykuroda 0:13a5d365ba16 70 { 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)); }
ykuroda 0:13a5d365ba16 71 };
ykuroda 0:13a5d365ba16 72
ykuroda 0:13a5d365ba16 73 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,false>
ykuroda 0:13a5d365ba16 74 {
ykuroda 0:13a5d365ba16 75 typedef std::complex<RealScalar> Scalar;
ykuroda 0:13a5d365ba16 76 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
ykuroda 0:13a5d365ba16 77 { return c + pmul(x,y); }
ykuroda 0:13a5d365ba16 78
ykuroda 0:13a5d365ba16 79 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
ykuroda 0:13a5d365ba16 80 { 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)); }
ykuroda 0:13a5d365ba16 81 };
ykuroda 0:13a5d365ba16 82
ykuroda 0:13a5d365ba16 83 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,true>
ykuroda 0:13a5d365ba16 84 {
ykuroda 0:13a5d365ba16 85 typedef std::complex<RealScalar> Scalar;
ykuroda 0:13a5d365ba16 86 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
ykuroda 0:13a5d365ba16 87 { return c + pmul(x,y); }
ykuroda 0:13a5d365ba16 88
ykuroda 0:13a5d365ba16 89 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
ykuroda 0:13a5d365ba16 90 { 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)); }
ykuroda 0:13a5d365ba16 91 };
ykuroda 0:13a5d365ba16 92
ykuroda 0:13a5d365ba16 93 template<typename RealScalar,bool Conj> struct conj_helper<std::complex<RealScalar>, RealScalar, Conj,false>
ykuroda 0:13a5d365ba16 94 {
ykuroda 0:13a5d365ba16 95 typedef std::complex<RealScalar> Scalar;
ykuroda 0:13a5d365ba16 96 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const
ykuroda 0:13a5d365ba16 97 { return padd(c, pmul(x,y)); }
ykuroda 0:13a5d365ba16 98 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const
ykuroda 0:13a5d365ba16 99 { return conj_if<Conj>()(x)*y; }
ykuroda 0:13a5d365ba16 100 };
ykuroda 0:13a5d365ba16 101
ykuroda 0:13a5d365ba16 102 template<typename RealScalar,bool Conj> struct conj_helper<RealScalar, std::complex<RealScalar>, false,Conj>
ykuroda 0:13a5d365ba16 103 {
ykuroda 0:13a5d365ba16 104 typedef std::complex<RealScalar> Scalar;
ykuroda 0:13a5d365ba16 105 EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const
ykuroda 0:13a5d365ba16 106 { return padd(c, pmul(x,y)); }
ykuroda 0:13a5d365ba16 107 EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const
ykuroda 0:13a5d365ba16 108 { return x*conj_if<Conj>()(y); }
ykuroda 0:13a5d365ba16 109 };
ykuroda 0:13a5d365ba16 110
ykuroda 0:13a5d365ba16 111 template<typename From,typename To> struct get_factor {
ykuroda 0:13a5d365ba16 112 static EIGEN_STRONG_INLINE To run(const From& x) { return x; }
ykuroda 0:13a5d365ba16 113 };
ykuroda 0:13a5d365ba16 114
ykuroda 0:13a5d365ba16 115 template<typename Scalar> struct get_factor<Scalar,typename NumTraits<Scalar>::Real> {
ykuroda 0:13a5d365ba16 116 static EIGEN_STRONG_INLINE typename NumTraits<Scalar>::Real run(const Scalar& x) { return numext::real(x); }
ykuroda 0:13a5d365ba16 117 };
ykuroda 0:13a5d365ba16 118
ykuroda 0:13a5d365ba16 119 // Lightweight helper class to access matrix coefficients.
ykuroda 0:13a5d365ba16 120 // Yes, this is somehow redundant with Map<>, but this version is much much lighter,
ykuroda 0:13a5d365ba16 121 // and so I hope better compilation performance (time and code quality).
ykuroda 0:13a5d365ba16 122 template<typename Scalar, typename Index, int StorageOrder>
ykuroda 0:13a5d365ba16 123 class blas_data_mapper
ykuroda 0:13a5d365ba16 124 {
ykuroda 0:13a5d365ba16 125 public:
ykuroda 0:13a5d365ba16 126 blas_data_mapper(Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
ykuroda 0:13a5d365ba16 127 EIGEN_STRONG_INLINE Scalar& operator()(Index i, Index j)
ykuroda 0:13a5d365ba16 128 { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; }
ykuroda 0:13a5d365ba16 129 protected:
ykuroda 0:13a5d365ba16 130 Scalar* EIGEN_RESTRICT m_data;
ykuroda 0:13a5d365ba16 131 Index m_stride;
ykuroda 0:13a5d365ba16 132 };
ykuroda 0:13a5d365ba16 133
ykuroda 0:13a5d365ba16 134 // lightweight helper class to access matrix coefficients (const version)
ykuroda 0:13a5d365ba16 135 template<typename Scalar, typename Index, int StorageOrder>
ykuroda 0:13a5d365ba16 136 class const_blas_data_mapper
ykuroda 0:13a5d365ba16 137 {
ykuroda 0:13a5d365ba16 138 public:
ykuroda 0:13a5d365ba16 139 const_blas_data_mapper(const Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
ykuroda 0:13a5d365ba16 140 EIGEN_STRONG_INLINE const Scalar& operator()(Index i, Index j) const
ykuroda 0:13a5d365ba16 141 { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; }
ykuroda 0:13a5d365ba16 142 protected:
ykuroda 0:13a5d365ba16 143 const Scalar* EIGEN_RESTRICT m_data;
ykuroda 0:13a5d365ba16 144 Index m_stride;
ykuroda 0:13a5d365ba16 145 };
ykuroda 0:13a5d365ba16 146
ykuroda 0:13a5d365ba16 147
ykuroda 0:13a5d365ba16 148 /* Helper class to analyze the factors of a Product expression.
ykuroda 0:13a5d365ba16 149 * In particular it allows to pop out operator-, scalar multiples,
ykuroda 0:13a5d365ba16 150 * and conjugate */
ykuroda 0:13a5d365ba16 151 template<typename XprType> struct blas_traits
ykuroda 0:13a5d365ba16 152 {
ykuroda 0:13a5d365ba16 153 typedef typename traits<XprType>::Scalar Scalar;
ykuroda 0:13a5d365ba16 154 typedef const XprType& ExtractType;
ykuroda 0:13a5d365ba16 155 typedef XprType _ExtractType;
ykuroda 0:13a5d365ba16 156 enum {
ykuroda 0:13a5d365ba16 157 IsComplex = NumTraits<Scalar>::IsComplex,
ykuroda 0:13a5d365ba16 158 IsTransposed = false,
ykuroda 0:13a5d365ba16 159 NeedToConjugate = false,
ykuroda 0:13a5d365ba16 160 HasUsableDirectAccess = ( (int(XprType::Flags)&DirectAccessBit)
ykuroda 0:13a5d365ba16 161 && ( bool(XprType::IsVectorAtCompileTime)
ykuroda 0:13a5d365ba16 162 || int(inner_stride_at_compile_time<XprType>::ret) == 1)
ykuroda 0:13a5d365ba16 163 ) ? 1 : 0
ykuroda 0:13a5d365ba16 164 };
ykuroda 0:13a5d365ba16 165 typedef typename conditional<bool(HasUsableDirectAccess),
ykuroda 0:13a5d365ba16 166 ExtractType,
ykuroda 0:13a5d365ba16 167 typename _ExtractType::PlainObject
ykuroda 0:13a5d365ba16 168 >::type DirectLinearAccessType;
ykuroda 0:13a5d365ba16 169 static inline ExtractType extract(const XprType& x) { return x; }
ykuroda 0:13a5d365ba16 170 static inline const Scalar extractScalarFactor(const XprType&) { return Scalar(1); }
ykuroda 0:13a5d365ba16 171 };
ykuroda 0:13a5d365ba16 172
ykuroda 0:13a5d365ba16 173 // pop conjugate
ykuroda 0:13a5d365ba16 174 template<typename Scalar, typename NestedXpr>
ykuroda 0:13a5d365ba16 175 struct blas_traits<CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> >
ykuroda 0:13a5d365ba16 176 : blas_traits<NestedXpr>
ykuroda 0:13a5d365ba16 177 {
ykuroda 0:13a5d365ba16 178 typedef blas_traits<NestedXpr> Base;
ykuroda 0:13a5d365ba16 179 typedef CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> XprType;
ykuroda 0:13a5d365ba16 180 typedef typename Base::ExtractType ExtractType;
ykuroda 0:13a5d365ba16 181
ykuroda 0:13a5d365ba16 182 enum {
ykuroda 0:13a5d365ba16 183 IsComplex = NumTraits<Scalar>::IsComplex,
ykuroda 0:13a5d365ba16 184 NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex
ykuroda 0:13a5d365ba16 185 };
ykuroda 0:13a5d365ba16 186 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 187 static inline Scalar extractScalarFactor(const XprType& x) { return conj(Base::extractScalarFactor(x.nestedExpression())); }
ykuroda 0:13a5d365ba16 188 };
ykuroda 0:13a5d365ba16 189
ykuroda 0:13a5d365ba16 190 // pop scalar multiple
ykuroda 0:13a5d365ba16 191 template<typename Scalar, typename NestedXpr>
ykuroda 0:13a5d365ba16 192 struct blas_traits<CwiseUnaryOp<scalar_multiple_op<Scalar>, NestedXpr> >
ykuroda 0:13a5d365ba16 193 : blas_traits<NestedXpr>
ykuroda 0:13a5d365ba16 194 {
ykuroda 0:13a5d365ba16 195 typedef blas_traits<NestedXpr> Base;
ykuroda 0:13a5d365ba16 196 typedef CwiseUnaryOp<scalar_multiple_op<Scalar>, NestedXpr> XprType;
ykuroda 0:13a5d365ba16 197 typedef typename Base::ExtractType ExtractType;
ykuroda 0:13a5d365ba16 198 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 199 static inline Scalar extractScalarFactor(const XprType& x)
ykuroda 0:13a5d365ba16 200 { return x.functor().m_other * Base::extractScalarFactor(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 201 };
ykuroda 0:13a5d365ba16 202
ykuroda 0:13a5d365ba16 203 // pop opposite
ykuroda 0:13a5d365ba16 204 template<typename Scalar, typename NestedXpr>
ykuroda 0:13a5d365ba16 205 struct blas_traits<CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> >
ykuroda 0:13a5d365ba16 206 : blas_traits<NestedXpr>
ykuroda 0:13a5d365ba16 207 {
ykuroda 0:13a5d365ba16 208 typedef blas_traits<NestedXpr> Base;
ykuroda 0:13a5d365ba16 209 typedef CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> XprType;
ykuroda 0:13a5d365ba16 210 typedef typename Base::ExtractType ExtractType;
ykuroda 0:13a5d365ba16 211 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 212 static inline Scalar extractScalarFactor(const XprType& x)
ykuroda 0:13a5d365ba16 213 { return - Base::extractScalarFactor(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 214 };
ykuroda 0:13a5d365ba16 215
ykuroda 0:13a5d365ba16 216 // pop/push transpose
ykuroda 0:13a5d365ba16 217 template<typename NestedXpr>
ykuroda 0:13a5d365ba16 218 struct blas_traits<Transpose<NestedXpr> >
ykuroda 0:13a5d365ba16 219 : blas_traits<NestedXpr>
ykuroda 0:13a5d365ba16 220 {
ykuroda 0:13a5d365ba16 221 typedef typename NestedXpr::Scalar Scalar;
ykuroda 0:13a5d365ba16 222 typedef blas_traits<NestedXpr> Base;
ykuroda 0:13a5d365ba16 223 typedef Transpose<NestedXpr> XprType;
ykuroda 0:13a5d365ba16 224 typedef Transpose<const typename Base::_ExtractType> ExtractType; // const to get rid of a compile error; anyway blas traits are only used on the RHS
ykuroda 0:13a5d365ba16 225 typedef Transpose<const typename Base::_ExtractType> _ExtractType;
ykuroda 0:13a5d365ba16 226 typedef typename conditional<bool(Base::HasUsableDirectAccess),
ykuroda 0:13a5d365ba16 227 ExtractType,
ykuroda 0:13a5d365ba16 228 typename ExtractType::PlainObject
ykuroda 0:13a5d365ba16 229 >::type DirectLinearAccessType;
ykuroda 0:13a5d365ba16 230 enum {
ykuroda 0:13a5d365ba16 231 IsTransposed = Base::IsTransposed ? 0 : 1
ykuroda 0:13a5d365ba16 232 };
ykuroda 0:13a5d365ba16 233 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 234 static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); }
ykuroda 0:13a5d365ba16 235 };
ykuroda 0:13a5d365ba16 236
ykuroda 0:13a5d365ba16 237 template<typename T>
ykuroda 0:13a5d365ba16 238 struct blas_traits<const T>
ykuroda 0:13a5d365ba16 239 : blas_traits<T>
ykuroda 0:13a5d365ba16 240 {};
ykuroda 0:13a5d365ba16 241
ykuroda 0:13a5d365ba16 242 template<typename T, bool HasUsableDirectAccess=blas_traits<T>::HasUsableDirectAccess>
ykuroda 0:13a5d365ba16 243 struct extract_data_selector {
ykuroda 0:13a5d365ba16 244 static const typename T::Scalar* run(const T& m)
ykuroda 0:13a5d365ba16 245 {
ykuroda 0:13a5d365ba16 246 return blas_traits<T>::extract(m).data();
ykuroda 0:13a5d365ba16 247 }
ykuroda 0:13a5d365ba16 248 };
ykuroda 0:13a5d365ba16 249
ykuroda 0:13a5d365ba16 250 template<typename T>
ykuroda 0:13a5d365ba16 251 struct extract_data_selector<T,false> {
ykuroda 0:13a5d365ba16 252 static typename T::Scalar* run(const T&) { return 0; }
ykuroda 0:13a5d365ba16 253 };
ykuroda 0:13a5d365ba16 254
ykuroda 0:13a5d365ba16 255 template<typename T> const typename T::Scalar* extract_data(const T& m)
ykuroda 0:13a5d365ba16 256 {
ykuroda 0:13a5d365ba16 257 return extract_data_selector<T>::run(m);
ykuroda 0:13a5d365ba16 258 }
ykuroda 0:13a5d365ba16 259
ykuroda 0:13a5d365ba16 260 } // end namespace internal
ykuroda 0:13a5d365ba16 261
ykuroda 0:13a5d365ba16 262 } // end namespace Eigen
ykuroda 0:13a5d365ba16 263
ykuroda 0:13a5d365ba16 264 #endif // EIGEN_BLASUTIL_H