SelfadjointMatrixMatrix_MKL.h
Go to the documentation of this file.
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 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 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 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 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
#define EIGEN_MKL_HEMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
#define EIGEN_MKL_SYMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX)
#define EIGEN_MKL_SYMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX)
#define EIGEN_MKL_HEMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:04