← Back to branch summary

~ubuntu-branches/ubuntu/saucy/blender/saucy-proposed

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/Core/products/SelfadjointMatrixMatrix_MKL.h

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/Eigen3/Eigen/src/Core/products/SelfadjointMatrixMatrix_MKL.h
 
1
/*
 
2
 Copyright (c) 2011, Intel Corporation. All rights reserved.
 
3
 
 
4
 Redistribution and use in source and binary forms, with or without modification,
 
5
 are permitted provided that the following conditions are met:
 
6
 
 
7
 * Redistributions of source code must retain the above copyright notice, this
 
8
   list of conditions and the following disclaimer.
 
9
 * Redistributions in binary form must reproduce the above copyright notice,
 
10
   this list of conditions and the following disclaimer in the documentation
 
11
   and/or other materials provided with the distribution.
 
12
 * Neither the name of Intel Corporation nor the names of its contributors may
 
13
   be used to endorse or promote products derived from this software without
 
14
   specific prior written permission.
 
15
 
 
16
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 
17
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 
18
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 
19
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 
20
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 
21
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 
22
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 
23
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 
24
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 
25
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
26
 
 
27
 ********************************************************************************
 
28
 *   Content : Eigen bindings to Intel(R) MKL
 
29
 *   Self adjoint matrix * matrix product functionality based on ?SYMM/?HEMM.
 
30
 ********************************************************************************
 
31
*/
 
32
 
 
33
#ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H
 
34
#define EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H
 
35
 
 
36
namespace Eigen { 
 
37
 
 
38
namespace internal {
 
39
 
 
40
 
 
41
/* Optimized selfadjoint matrix * matrix (?SYMM/?HEMM) product */
 
42
 
 
43
#define EIGEN_MKL_SYMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
 
44
template <typename Index, \
 
45
          int LhsStorageOrder, bool ConjugateLhs, \
 
46
          int RhsStorageOrder, bool ConjugateRhs> \
 
47
struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,true,ConjugateLhs,RhsStorageOrder,false,ConjugateRhs,ColMajor> \
 
48
{\
 
49
\
 
50
  static EIGEN_DONT_INLINE void run( \
 
51
    Index rows, Index cols, \
 
52
    const EIGTYPE* _lhs, Index lhsStride, \
 
53
    const EIGTYPE* _rhs, Index rhsStride, \
 
54
    EIGTYPE* res,        Index resStride, \
 
55
    EIGTYPE alpha) \
 
56
  { \
 
57
    char side='L', uplo='L'; \
 
58
    MKL_INT m, n, lda, ldb, ldc; \
 
59
    const EIGTYPE *a, *b; \
 
60
    MKLTYPE alpha_, beta_; \
 
61
    MatrixX##EIGPREFIX b_tmp; \
 
62
    EIGTYPE myone(1);\
 
63
\
 
64
/* Set transpose options */ \
 
65
/* Set m, n, k */ \
 
66
    m = (MKL_INT)rows;  \
 
67
    n = (MKL_INT)cols;  \
 
68
\
 
69
/* Set alpha_ & beta_ */ \
 
70
    assign_scalar_eig2mkl(alpha_, alpha); \
 
71
    assign_scalar_eig2mkl(beta_, myone); \
 
72
\
 
73
/* Set lda, ldb, ldc */ \
 
74
    lda = (MKL_INT)lhsStride; \
 
75
    ldb = (MKL_INT)rhsStride; \
 
76
    ldc = (MKL_INT)resStride; \
 
77
\
 
78
/* Set a, b, c */ \
 
79
    if (LhsStorageOrder==RowMajor) uplo='U'; \
 
80
    a = _lhs; \
 
81
\
 
82
    if (RhsStorageOrder==RowMajor) { \
 
83
      Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \
 
84
      b_tmp = rhs.adjoint(); \
 
85
      b = b_tmp.data(); \
 
86
      ldb = b_tmp.outerStride(); \
 
87
    } else b = _rhs; \
 
88
\
 
89
    MKLPREFIX##symm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
 
90
\
 
91
  } \
 
92
};
 
93
 
 
94
 
 
95
#define EIGEN_MKL_HEMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
 
96
template <typename Index, \
 
97
          int LhsStorageOrder, bool ConjugateLhs, \
 
98
          int RhsStorageOrder, bool ConjugateRhs> \
 
99
struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,true,ConjugateLhs,RhsStorageOrder,false,ConjugateRhs,ColMajor> \
 
100
{\
 
101
  static EIGEN_DONT_INLINE void run( \
 
102
    Index rows, Index cols, \
 
103
    const EIGTYPE* _lhs, Index lhsStride, \
 
104
    const EIGTYPE* _rhs, Index rhsStride, \
 
105
    EIGTYPE* res,        Index resStride, \
 
106
    EIGTYPE alpha) \
 
107
  { \
 
108
    char side='L', uplo='L'; \
 
109
    MKL_INT m, n, lda, ldb, ldc; \
 
110
    const EIGTYPE *a, *b; \
 
111
    MKLTYPE alpha_, beta_; \
 
112
    MatrixX##EIGPREFIX b_tmp; \
 
113
    Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder> a_tmp; \
 
114
    EIGTYPE myone(1); \
 
115
\
 
116
/* Set transpose options */ \
 
117
/* Set m, n, k */ \
 
118
    m = (MKL_INT)rows; \
 
119
    n = (MKL_INT)cols; \
 
120
\
 
121
/* Set alpha_ & beta_ */ \
 
122
    assign_scalar_eig2mkl(alpha_, alpha); \
 
123
    assign_scalar_eig2mkl(beta_, myone); \
 
124
\
 
125
/* Set lda, ldb, ldc */ \
 
126
    lda = (MKL_INT)lhsStride; \
 
127
    ldb = (MKL_INT)rhsStride; \
 
128
    ldc = (MKL_INT)resStride; \
 
129
\
 
130
/* Set a, b, c */ \
 
131
    if (((LhsStorageOrder==ColMajor) && ConjugateLhs) || ((LhsStorageOrder==RowMajor) && (!ConjugateLhs))) { \
 
132
      Map<const Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder>, 0, OuterStride<> > lhs(_lhs,m,m,OuterStride<>(lhsStride)); \
 
133
      a_tmp = lhs.conjugate(); \
 
134
      a = a_tmp.data(); \
 
135
      lda = a_tmp.outerStride(); \
 
136
    } else a = _lhs; \
 
137
    if (LhsStorageOrder==RowMajor) uplo='U'; \
 
138
\
 
139
    if (RhsStorageOrder==ColMajor && (!ConjugateRhs)) { \
 
140
       b = _rhs; } \
 
141
    else { \
 
142
      if (RhsStorageOrder==ColMajor && ConjugateRhs) { \
 
143
        Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,m,n,OuterStride<>(rhsStride)); \
 
144
        b_tmp = rhs.conjugate(); \
 
145
      } else \
 
146
      if (ConjugateRhs) { \
 
147
        Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \
 
148
        b_tmp = rhs.adjoint(); \
 
149
      } else { \
 
150
        Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \
 
151
        b_tmp = rhs.transpose(); \
 
152
      } \
 
153
      b = b_tmp.data(); \
 
154
      ldb = b_tmp.outerStride(); \
 
155
    } \
 
156
\
 
157
    MKLPREFIX##hemm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
 
158
\
 
159
  } \
 
160
};
 
161
 
 
162
EIGEN_MKL_SYMM_L(double, double, d, d)
 
163
EIGEN_MKL_SYMM_L(float, float, f, s)
 
164
EIGEN_MKL_HEMM_L(dcomplex, MKL_Complex16, cd, z)
 
165
EIGEN_MKL_HEMM_L(scomplex, MKL_Complex8, cf, c)
 
166
 
 
167
 
 
168
/* Optimized matrix * selfadjoint matrix (?SYMM/?HEMM) product */
 
169
 
 
170
#define EIGEN_MKL_SYMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
 
171
template <typename Index, \
 
172
          int LhsStorageOrder, bool ConjugateLhs, \
 
173
          int RhsStorageOrder, bool ConjugateRhs> \
 
174
struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,false,ConjugateLhs,RhsStorageOrder,true,ConjugateRhs,ColMajor> \
 
175
{\
 
176
\
 
177
  static EIGEN_DONT_INLINE void run( \
 
178
    Index rows, Index cols, \
 
179
    const EIGTYPE* _lhs, Index lhsStride, \
 
180
    const EIGTYPE* _rhs, Index rhsStride, \
 
181
    EIGTYPE* res,        Index resStride, \
 
182
    EIGTYPE alpha) \
 
183
  { \
 
184
    char side='R', uplo='L'; \
 
185
    MKL_INT m, n, lda, ldb, ldc; \
 
186
    const EIGTYPE *a, *b; \
 
187
    MKLTYPE alpha_, beta_; \
 
188
    MatrixX##EIGPREFIX b_tmp; \
 
189
    EIGTYPE myone(1);\
 
190
\
 
191
/* Set m, n, k */ \
 
192
    m = (MKL_INT)rows;  \
 
193
    n = (MKL_INT)cols;  \
 
194
\
 
195
/* Set alpha_ & beta_ */ \
 
196
    assign_scalar_eig2mkl(alpha_, alpha); \
 
197
    assign_scalar_eig2mkl(beta_, myone); \
 
198
\
 
199
/* Set lda, ldb, ldc */ \
 
200
    lda = (MKL_INT)rhsStride; \
 
201
    ldb = (MKL_INT)lhsStride; \
 
202
    ldc = (MKL_INT)resStride; \
 
203
\
 
204
/* Set a, b, c */ \
 
205
    if (RhsStorageOrder==RowMajor) uplo='U'; \
 
206
    a = _rhs; \
 
207
\
 
208
    if (LhsStorageOrder==RowMajor) { \
 
209
      Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(rhsStride)); \
 
210
      b_tmp = lhs.adjoint(); \
 
211
      b = b_tmp.data(); \
 
212
      ldb = b_tmp.outerStride(); \
 
213
    } else b = _lhs; \
 
214
\
 
215
    MKLPREFIX##symm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
 
216
\
 
217
  } \
 
218
};
 
219
 
 
220
 
 
221
#define EIGEN_MKL_HEMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
 
222
template <typename Index, \
 
223
          int LhsStorageOrder, bool ConjugateLhs, \
 
224
          int RhsStorageOrder, bool ConjugateRhs> \
 
225
struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,false,ConjugateLhs,RhsStorageOrder,true,ConjugateRhs,ColMajor> \
 
226
{\
 
227
  static EIGEN_DONT_INLINE void run( \
 
228
    Index rows, Index cols, \
 
229
    const EIGTYPE* _lhs, Index lhsStride, \
 
230
    const EIGTYPE* _rhs, Index rhsStride, \
 
231
    EIGTYPE* res,        Index resStride, \
 
232
    EIGTYPE alpha) \
 
233
  { \
 
234
    char side='R', uplo='L'; \
 
235
    MKL_INT m, n, lda, ldb, ldc; \
 
236
    const EIGTYPE *a, *b; \
 
237
    MKLTYPE alpha_, beta_; \
 
238
    MatrixX##EIGPREFIX b_tmp; \
 
239
    Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder> a_tmp; \
 
240
    EIGTYPE myone(1); \
 
241
\
 
242
/* Set m, n, k */ \
 
243
    m = (MKL_INT)rows; \
 
244
    n = (MKL_INT)cols; \
 
245
\
 
246
/* Set alpha_ & beta_ */ \
 
247
    assign_scalar_eig2mkl(alpha_, alpha); \
 
248
    assign_scalar_eig2mkl(beta_, myone); \
 
249
\
 
250
/* Set lda, ldb, ldc */ \
 
251
    lda = (MKL_INT)rhsStride; \
 
252
    ldb = (MKL_INT)lhsStride; \
 
253
    ldc = (MKL_INT)resStride; \
 
254
\
 
255
/* Set a, b, c */ \
 
256
    if (((RhsStorageOrder==ColMajor) && ConjugateRhs) || ((RhsStorageOrder==RowMajor) && (!ConjugateRhs))) { \
 
257
      Map<const Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder>, 0, OuterStride<> > rhs(_rhs,n,n,OuterStride<>(rhsStride)); \
 
258
      a_tmp = rhs.conjugate(); \
 
259
      a = a_tmp.data(); \
 
260
      lda = a_tmp.outerStride(); \
 
261
    } else a = _rhs; \
 
262
    if (RhsStorageOrder==RowMajor) uplo='U'; \
 
263
\
 
264
    if (LhsStorageOrder==ColMajor && (!ConjugateLhs)) { \
 
265
       b = _lhs; } \
 
266
    else { \
 
267
      if (LhsStorageOrder==ColMajor && ConjugateLhs) { \
 
268
        Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,m,n,OuterStride<>(lhsStride)); \
 
269
        b_tmp = lhs.conjugate(); \
 
270
      } else \
 
271
      if (ConjugateLhs) { \
 
272
        Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \
 
273
        b_tmp = lhs.adjoint(); \
 
274
      } else { \
 
275
        Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \
 
276
        b_tmp = lhs.transpose(); \
 
277
      } \
 
278
      b = b_tmp.data(); \
 
279
      ldb = b_tmp.outerStride(); \
 
280
    } \
 
281
\
 
282
    MKLPREFIX##hemm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
 
283
  } \
 
284
};
 
285
 
 
286
EIGEN_MKL_SYMM_R(double, double, d, d)
 
287
EIGEN_MKL_SYMM_R(float, float, f, s)
 
288
EIGEN_MKL_HEMM_R(dcomplex, MKL_Complex16, cd, z)
 
289
EIGEN_MKL_HEMM_R(scomplex, MKL_Complex8, cf, c)
 
290
 
 
291
} // end namespace internal
 
292
 
 
293
} // end namespace Eigen
 
294
 
 
295
#endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H