Eigen  3.2.10
BlasUtil.h
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_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, int mr, int nr, bool ConjugateLhs=false, bool ConjugateRhs=false>
22 struct gebp_kernel;
23 
24 template<typename Scalar, typename Index, int nr, int StorageOrder, bool Conjugate = false, bool PanelMode=false>
25 struct gemm_pack_rhs;
26 
27 template<typename Scalar, typename Index, 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, typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version=Specialized>
38 struct general_matrix_vector_product;
39 
40 
41 template<bool Conjugate> struct conj_if;
42 
43 template<> struct conj_if<true> {
44  template<typename T>
45  inline T operator()(const T& x) const { return numext::conj(x); }
46  template<typename T>
47  inline T pconj(const T& x) const { return internal::pconj(x); }
48 };
49 
50 template<> struct conj_if<false> {
51  template<typename T>
52  inline const T& operator()(const T& x) const { return x; }
53  template<typename T>
54  inline const T& pconj(const T& x) const { return x; }
55 };
56 
57 // Generic implementation for custom complex types.
58 template<typename LhsScalar, typename RhsScalar, bool ConjLhs, bool ConjRhs>
59 struct conj_helper
60 {
61  typedef typename scalar_product_traits<LhsScalar,RhsScalar>::ReturnType Scalar;
62 
63  EIGEN_STRONG_INLINE Scalar pmadd(const LhsScalar& x, const RhsScalar& y, const Scalar& c) const
64  { return padd(c, pmul(x,y)); }
65 
66  EIGEN_STRONG_INLINE Scalar pmul(const LhsScalar& x, const RhsScalar& y) const
67  { return conj_if<ConjLhs>()(x) * conj_if<ConjRhs>()(y); }
68 };
69 
70 template<typename Scalar> struct conj_helper<Scalar,Scalar,false,false>
71 {
72  EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const { return internal::pmadd(x,y,c); }
73  EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const { return internal::pmul(x,y); }
74 };
75 
76 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, false,true>
77 {
78  typedef std::complex<RealScalar> Scalar;
79  EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
80  { return c + pmul(x,y); }
81 
82  EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
83  { 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)); }
84 };
85 
86 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,false>
87 {
88  typedef std::complex<RealScalar> Scalar;
89  EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
90  { return c + pmul(x,y); }
91 
92  EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
93  { 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)); }
94 };
95 
96 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,true>
97 {
98  typedef std::complex<RealScalar> Scalar;
99  EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
100  { return c + pmul(x,y); }
101 
102  EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
103  { 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)); }
104 };
105 
106 template<typename RealScalar,bool Conj> struct conj_helper<std::complex<RealScalar>, RealScalar, Conj,false>
107 {
108  typedef std::complex<RealScalar> Scalar;
109  EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const
110  { return padd(c, pmul(x,y)); }
111  EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const
112  { return conj_if<Conj>()(x)*y; }
113 };
114 
115 template<typename RealScalar,bool Conj> struct conj_helper<RealScalar, std::complex<RealScalar>, false,Conj>
116 {
117  typedef std::complex<RealScalar> Scalar;
118  EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const
119  { return padd(c, pmul(x,y)); }
120  EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const
121  { return x*conj_if<Conj>()(y); }
122 };
123 
124 template<typename From,typename To> struct get_factor {
125  static EIGEN_STRONG_INLINE To run(const From& x) { return x; }
126 };
127 
128 template<typename Scalar> struct get_factor<Scalar,typename NumTraits<Scalar>::Real> {
129  static EIGEN_STRONG_INLINE typename NumTraits<Scalar>::Real run(const Scalar& x) { return numext::real(x); }
130 };
131 
132 // Lightweight helper class to access matrix coefficients.
133 // Yes, this is somehow redundant with Map<>, but this version is much much lighter,
134 // and so I hope better compilation performance (time and code quality).
135 template<typename Scalar, typename Index, int StorageOrder>
136 class blas_data_mapper
137 {
138  public:
139  blas_data_mapper(Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
140  EIGEN_STRONG_INLINE Scalar& operator()(Index i, Index j)
141  { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; }
142  protected:
143  Scalar* EIGEN_RESTRICT m_data;
144  Index m_stride;
145 };
146 
147 // lightweight helper class to access matrix coefficients (const version)
148 template<typename Scalar, typename Index, int StorageOrder>
149 class const_blas_data_mapper
150 {
151  public:
152  const_blas_data_mapper(const Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
153  EIGEN_STRONG_INLINE const Scalar& operator()(Index i, Index j) const
154  { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; }
155  protected:
156  const Scalar* EIGEN_RESTRICT m_data;
157  Index m_stride;
158 };
159 
160 
161 /* Helper class to analyze the factors of a Product expression.
162  * In particular it allows to pop out operator-, scalar multiples,
163  * and conjugate */
164 template<typename XprType> struct blas_traits
165 {
166  typedef typename traits<XprType>::Scalar Scalar;
167  typedef const XprType& ExtractType;
168  typedef XprType _ExtractType;
169  enum {
170  IsComplex = NumTraits<Scalar>::IsComplex,
171  IsTransposed = false,
172  NeedToConjugate = false,
173  HasUsableDirectAccess = ( (int(XprType::Flags)&DirectAccessBit)
174  && ( bool(XprType::IsVectorAtCompileTime)
175  || int(inner_stride_at_compile_time<XprType>::ret) == 1)
176  ) ? 1 : 0
177  };
178  typedef typename conditional<bool(HasUsableDirectAccess),
179  ExtractType,
180  typename _ExtractType::PlainObject
181  >::type DirectLinearAccessType;
182  static inline ExtractType extract(const XprType& x) { return x; }
183  static inline const Scalar extractScalarFactor(const XprType&) { return Scalar(1); }
184 };
185 
186 // pop conjugate
187 template<typename Scalar, typename Xpr>
188 struct blas_traits<CwiseUnaryOp<scalar_conjugate_op<Scalar>, Xpr> >
189  : blas_traits<typename internal::remove_all<typename Xpr::Nested>::type>
190 {
191  typedef typename internal::remove_all<typename Xpr::Nested>::type NestedXpr;
192  typedef blas_traits<NestedXpr> Base;
193  typedef CwiseUnaryOp<scalar_conjugate_op<Scalar>, Xpr> XprType;
194  typedef typename Base::ExtractType ExtractType;
195 
196  enum {
197  IsComplex = NumTraits<Scalar>::IsComplex,
198  NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex
199  };
200  static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
201  static inline Scalar extractScalarFactor(const XprType& x) { return conj(Base::extractScalarFactor(x.nestedExpression())); }
202 };
203 
204 // pop scalar multiple
205 template<typename Scalar, typename Xpr>
206 struct blas_traits<CwiseUnaryOp<scalar_multiple_op<Scalar>, Xpr> >
207  : blas_traits<typename internal::remove_all<typename Xpr::Nested>::type>
208 {
209  typedef typename internal::remove_all<typename Xpr::Nested>::type NestedXpr;
210  typedef blas_traits<NestedXpr> Base;
211  typedef CwiseUnaryOp<scalar_multiple_op<Scalar>, Xpr> XprType;
212  typedef typename Base::ExtractType ExtractType;
213  static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
214  static inline Scalar extractScalarFactor(const XprType& x)
215  { return x.functor().m_other * Base::extractScalarFactor(x.nestedExpression()); }
216 };
217 
218 // pop opposite
219 template<typename Scalar, typename Xpr>
220 struct blas_traits<CwiseUnaryOp<scalar_opposite_op<Scalar>, Xpr> >
221  : blas_traits<typename internal::remove_all<typename Xpr::Nested>::type>
222 {
223  typedef typename internal::remove_all<typename Xpr::Nested>::type NestedXpr;
224  typedef blas_traits<NestedXpr> Base;
225  typedef CwiseUnaryOp<scalar_opposite_op<Scalar>, Xpr> XprType;
226  typedef typename Base::ExtractType ExtractType;
227  static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
228  static inline Scalar extractScalarFactor(const XprType& x)
229  { return - Base::extractScalarFactor(x.nestedExpression()); }
230 };
231 
232 // pop/push transpose
233 template<typename Xpr>
234 struct blas_traits<Transpose<Xpr> >
235  : blas_traits<typename internal::remove_all<typename Xpr::Nested>::type>
236 {
237  typedef typename internal::remove_all<typename Xpr::Nested>::type NestedXpr;
238  typedef typename NestedXpr::Scalar Scalar;
239  typedef blas_traits<NestedXpr> Base;
240  typedef Transpose<Xpr> XprType;
241  typedef Transpose<const typename Base::_ExtractType> ExtractType; // const to get rid of a compile error; anyway blas traits are only used on the RHS
242  typedef Transpose<const typename Base::_ExtractType> _ExtractType;
243  typedef typename conditional<bool(Base::HasUsableDirectAccess),
244  ExtractType,
245  typename ExtractType::PlainObject
246  >::type DirectLinearAccessType;
247  enum {
248  IsTransposed = Base::IsTransposed ? 0 : 1
249  };
250  static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
251  static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); }
252 };
253 
254 template<typename T>
255 struct blas_traits<const T>
256  : blas_traits<T>
257 {};
258 
259 template<typename T, bool HasUsableDirectAccess=blas_traits<T>::HasUsableDirectAccess>
260 struct extract_data_selector {
261  static const typename T::Scalar* run(const T& m)
262  {
263  return blas_traits<T>::extract(m).data();
264  }
265 };
266 
267 template<typename T>
268 struct extract_data_selector<T,false> {
269  static typename T::Scalar* run(const T&) { return 0; }
270 };
271 
272 template<typename T> const typename T::Scalar* extract_data(const T& m)
273 {
274  return extract_data_selector<T>::run(m);
275 }
276 
277 } // end namespace internal
278 
279 } // end namespace Eigen
280 
281 #endif // EIGEN_BLASUTIL_H
Definition: LDLT.h:16
Definition: StdDeque.h:50
Definition: Constants.h:266
Definition: Eigen_Colamd.h:50
const unsigned int DirectAccessBit
Definition: Constants.h:142