← Back to branch summary

~choreonoid/choreonoid/debian

« back to all changes in this revision

Viewing changes to thirdparty/eigen3/Eigen/src/Core/util/BlasUtil.h

  • Committer: Thomas Moulard
  • Date: 2012-10-23 12:43:24 UTC
  • Revision ID: git-v1:351cf736ad49bc7a9a7b9767dee760a013517a5d
Tags: upstream/1.1.0
Imported Upstream version 1.1.0

Show diffs side-by-side

added added

removed removed

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