SelfadjointMatrixMatrix_BLAS.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 BLAS F77
29  * Self adjoint matrix * matrix product functionality based on ?SYMM/?HEMM.
30  ********************************************************************************
31 */
32 
33 #ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_BLAS_H
34 #define EIGEN_SELFADJOINT_MATRIX_MATRIX_BLAS_H
35 
36 // IWYU pragma: private
37 #include "../InternalHeaderCheck.h"
38 
39 namespace Eigen {
40 
41 namespace internal {
42 
43 /* Optimized selfadjoint matrix * matrix (?SYMM/?HEMM) product */
44 
45 #define EIGEN_BLAS_SYMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC) \
46  template <typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> \
47  struct product_selfadjoint_matrix<EIGTYPE, Index, LhsStorageOrder, true, ConjugateLhs, RhsStorageOrder, false, \
48  ConjugateRhs, ColMajor, 1> { \
49  static void run(Index rows, Index cols, const EIGTYPE* _lhs, Index lhsStride, const EIGTYPE* _rhs, \
50  Index rhsStride, EIGTYPE* res, Index resIncr, Index resStride, EIGTYPE alpha, \
51  level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/) { \
52  if (rows == 0 || cols == 0) return; \
53  EIGEN_ONLY_USED_FOR_DEBUG(resIncr); \
54  eigen_assert(resIncr == 1); \
55  char side = 'L', uplo = 'L'; \
56  BlasIndex m, n, lda, ldb, ldc; \
57  const EIGTYPE *a, *b; \
58  EIGTYPE beta(1); \
59  MatrixX##EIGPREFIX b_tmp; \
60  \
61  /* Set transpose options */ \
62  /* Set m, n, k */ \
63  m = convert_index<BlasIndex>(rows); \
64  n = convert_index<BlasIndex>(cols); \
65  \
66  /* Set lda, ldb, ldc */ \
67  lda = convert_index<BlasIndex>(lhsStride); \
68  ldb = convert_index<BlasIndex>(rhsStride); \
69  ldc = convert_index<BlasIndex>(resStride); \
70  \
71  /* Set a, b, c */ \
72  if (LhsStorageOrder == RowMajor) uplo = 'U'; \
73  a = _lhs; \
74  \
75  if (RhsStorageOrder == RowMajor) { \
76  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs, n, m, OuterStride<>(rhsStride)); \
77  b_tmp = rhs.adjoint(); \
78  b = b_tmp.data(); \
79  ldb = convert_index<BlasIndex>(b_tmp.outerStride()); \
80  } else \
81  b = _rhs; \
82  \
83  BLASFUNC(&side, &uplo, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, \
84  (const BLASTYPE*)b, &ldb, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)res, &ldc); \
85  } \
86  };
87 
88 #define EIGEN_BLAS_HEMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC) \
89  template <typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> \
90  struct product_selfadjoint_matrix<EIGTYPE, Index, LhsStorageOrder, true, ConjugateLhs, RhsStorageOrder, false, \
91  ConjugateRhs, ColMajor, 1> { \
92  static void run(Index rows, Index cols, const EIGTYPE* _lhs, Index lhsStride, const EIGTYPE* _rhs, \
93  Index rhsStride, EIGTYPE* res, Index resIncr, Index resStride, EIGTYPE alpha, \
94  level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/) { \
95  if (rows == 0 || cols == 0) return; \
96  EIGEN_ONLY_USED_FOR_DEBUG(resIncr); \
97  eigen_assert(resIncr == 1); \
98  char side = 'L', uplo = 'L'; \
99  BlasIndex m, n, lda, ldb, ldc; \
100  const EIGTYPE *a, *b; \
101  EIGTYPE beta(1); \
102  MatrixX##EIGPREFIX b_tmp; \
103  Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder> a_tmp; \
104  \
105  /* Set transpose options */ \
106  /* Set m, n, k */ \
107  m = convert_index<BlasIndex>(rows); \
108  n = convert_index<BlasIndex>(cols); \
109  \
110  /* Set lda, ldb, ldc */ \
111  lda = convert_index<BlasIndex>(lhsStride); \
112  ldb = convert_index<BlasIndex>(rhsStride); \
113  ldc = convert_index<BlasIndex>(resStride); \
114  \
115  /* Set a, b, c */ \
116  if (((LhsStorageOrder == ColMajor) && ConjugateLhs) || ((LhsStorageOrder == RowMajor) && (!ConjugateLhs))) { \
117  Map<const Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder>, 0, OuterStride<> > lhs( \
118  _lhs, m, m, OuterStride<>(lhsStride)); \
119  a_tmp = lhs.conjugate(); \
120  a = a_tmp.data(); \
121  lda = convert_index<BlasIndex>(a_tmp.outerStride()); \
122  } else \
123  a = _lhs; \
124  if (LhsStorageOrder == RowMajor) uplo = 'U'; \
125  \
126  if (RhsStorageOrder == ColMajor && (!ConjugateRhs)) { \
127  b = _rhs; \
128  } else { \
129  if (RhsStorageOrder == ColMajor && ConjugateRhs) { \
130  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs, m, n, OuterStride<>(rhsStride)); \
131  b_tmp = rhs.conjugate(); \
132  } else if (ConjugateRhs) { \
133  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs, n, m, OuterStride<>(rhsStride)); \
134  b_tmp = rhs.adjoint(); \
135  } else { \
136  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs, n, m, OuterStride<>(rhsStride)); \
137  b_tmp = rhs.transpose(); \
138  } \
139  b = b_tmp.data(); \
140  ldb = convert_index<BlasIndex>(b_tmp.outerStride()); \
141  } \
142  \
143  BLASFUNC(&side, &uplo, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, \
144  (const BLASTYPE*)b, &ldb, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)res, &ldc); \
145  } \
146  };
147 
148 #ifdef EIGEN_USE_MKL
149 EIGEN_BLAS_SYMM_L(double, double, d, dsymm)
150 EIGEN_BLAS_SYMM_L(float, float, f, ssymm)
151 EIGEN_BLAS_HEMM_L(dcomplex, MKL_Complex16, cd, zhemm)
152 EIGEN_BLAS_HEMM_L(scomplex, MKL_Complex8, cf, chemm)
153 #else
154 EIGEN_BLAS_SYMM_L(double, double, d, dsymm_)
155 EIGEN_BLAS_SYMM_L(float, float, f, ssymm_)
156 EIGEN_BLAS_HEMM_L(dcomplex, double, cd, zhemm_)
157 EIGEN_BLAS_HEMM_L(scomplex, float, cf, chemm_)
158 #endif
159 
160 /* Optimized matrix * selfadjoint matrix (?SYMM/?HEMM) product */
161 
162 #define EIGEN_BLAS_SYMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC) \
163  template <typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> \
164  struct product_selfadjoint_matrix<EIGTYPE, Index, LhsStorageOrder, false, ConjugateLhs, RhsStorageOrder, true, \
165  ConjugateRhs, ColMajor, 1> { \
166  static void run(Index rows, Index cols, const EIGTYPE* _lhs, Index lhsStride, const EIGTYPE* _rhs, \
167  Index rhsStride, EIGTYPE* res, Index resIncr, Index resStride, EIGTYPE alpha, \
168  level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/) { \
169  if (rows == 0 || cols == 0) return; \
170  EIGEN_ONLY_USED_FOR_DEBUG(resIncr); \
171  eigen_assert(resIncr == 1); \
172  char side = 'R', uplo = 'L'; \
173  BlasIndex m, n, lda, ldb, ldc; \
174  const EIGTYPE *a, *b; \
175  EIGTYPE beta(1); \
176  MatrixX##EIGPREFIX b_tmp; \
177  \
178  /* Set m, n, k */ \
179  m = convert_index<BlasIndex>(rows); \
180  n = convert_index<BlasIndex>(cols); \
181  \
182  /* Set lda, ldb, ldc */ \
183  lda = convert_index<BlasIndex>(rhsStride); \
184  ldb = convert_index<BlasIndex>(lhsStride); \
185  ldc = convert_index<BlasIndex>(resStride); \
186  \
187  /* Set a, b, c */ \
188  if (RhsStorageOrder == RowMajor) uplo = 'U'; \
189  a = _rhs; \
190  \
191  if (LhsStorageOrder == RowMajor) { \
192  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs, n, m, OuterStride<>(rhsStride)); \
193  b_tmp = lhs.adjoint(); \
194  b = b_tmp.data(); \
195  ldb = convert_index<BlasIndex>(b_tmp.outerStride()); \
196  } else \
197  b = _lhs; \
198  \
199  BLASFUNC(&side, &uplo, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, \
200  (const BLASTYPE*)b, &ldb, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)res, &ldc); \
201  } \
202  };
203 
204 #define EIGEN_BLAS_HEMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC) \
205  template <typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> \
206  struct product_selfadjoint_matrix<EIGTYPE, Index, LhsStorageOrder, false, ConjugateLhs, RhsStorageOrder, true, \
207  ConjugateRhs, ColMajor, 1> { \
208  static void run(Index rows, Index cols, const EIGTYPE* _lhs, Index lhsStride, const EIGTYPE* _rhs, \
209  Index rhsStride, EIGTYPE* res, Index resIncr, Index resStride, EIGTYPE alpha, \
210  level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/) { \
211  EIGEN_ONLY_USED_FOR_DEBUG(resIncr); \
212  eigen_assert(resIncr == 1); \
213  char side = 'R', uplo = 'L'; \
214  BlasIndex m, n, lda, ldb, ldc; \
215  const EIGTYPE *a, *b; \
216  EIGTYPE beta(1); \
217  MatrixX##EIGPREFIX b_tmp; \
218  Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder> a_tmp; \
219  \
220  /* Set m, n, k */ \
221  m = convert_index<BlasIndex>(rows); \
222  n = convert_index<BlasIndex>(cols); \
223  \
224  /* Set lda, ldb, ldc */ \
225  lda = convert_index<BlasIndex>(rhsStride); \
226  ldb = convert_index<BlasIndex>(lhsStride); \
227  ldc = convert_index<BlasIndex>(resStride); \
228  \
229  /* Set a, b, c */ \
230  if (((RhsStorageOrder == ColMajor) && ConjugateRhs) || ((RhsStorageOrder == RowMajor) && (!ConjugateRhs))) { \
231  Map<const Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder>, 0, OuterStride<> > rhs( \
232  _rhs, n, n, OuterStride<>(rhsStride)); \
233  a_tmp = rhs.conjugate(); \
234  a = a_tmp.data(); \
235  lda = convert_index<BlasIndex>(a_tmp.outerStride()); \
236  } else \
237  a = _rhs; \
238  if (RhsStorageOrder == RowMajor) uplo = 'U'; \
239  \
240  if (LhsStorageOrder == ColMajor && (!ConjugateLhs)) { \
241  b = _lhs; \
242  } else { \
243  if (LhsStorageOrder == ColMajor && ConjugateLhs) { \
244  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs, m, n, OuterStride<>(lhsStride)); \
245  b_tmp = lhs.conjugate(); \
246  } else if (ConjugateLhs) { \
247  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs, n, m, OuterStride<>(lhsStride)); \
248  b_tmp = lhs.adjoint(); \
249  } else { \
250  Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs, n, m, OuterStride<>(lhsStride)); \
251  b_tmp = lhs.transpose(); \
252  } \
253  b = b_tmp.data(); \
254  ldb = convert_index<BlasIndex>(b_tmp.outerStride()); \
255  } \
256  \
257  BLASFUNC(&side, &uplo, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, \
258  (const BLASTYPE*)b, &ldb, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)res, &ldc); \
259  } \
260  };
261 
262 #ifdef EIGEN_USE_MKL
263 EIGEN_BLAS_SYMM_R(double, double, d, dsymm)
264 EIGEN_BLAS_SYMM_R(float, float, f, ssymm)
265 EIGEN_BLAS_HEMM_R(dcomplex, MKL_Complex16, cd, zhemm)
266 EIGEN_BLAS_HEMM_R(scomplex, MKL_Complex8, cf, chemm)
267 #else
268 EIGEN_BLAS_SYMM_R(double, double, d, dsymm_)
269 EIGEN_BLAS_SYMM_R(float, float, f, ssymm_)
270 EIGEN_BLAS_HEMM_R(dcomplex, double, cd, zhemm_)
271 EIGEN_BLAS_HEMM_R(scomplex, float, cf, chemm_)
272 #endif
273 } // end namespace internal
274 
275 } // end namespace Eigen
276 
277 #endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_BLAS_H
EIGEN_BLAS_HEMM_L
#define EIGEN_BLAS_HEMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC)
Definition: SelfadjointMatrixMatrix_BLAS.h:88
EIGEN_BLAS_HEMM_R
#define EIGEN_BLAS_HEMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC)
Definition: SelfadjointMatrixMatrix_BLAS.h:204
EIGEN_BLAS_SYMM_R
#define EIGEN_BLAS_SYMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC)
Definition: SelfadjointMatrixMatrix_BLAS.h:162
EIGEN_BLAS_SYMM_L
#define EIGEN_BLAS_SYMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASFUNC)
Definition: SelfadjointMatrixMatrix_BLAS.h:45
ssymm
int BLASFUNC() ssymm(char *, char *, int *, int *, float *, float *, int *, float *, int *, float *, float *, int *)
dsymm
int BLASFUNC() dsymm(char *, char *, int *, int *, double *, double *, int *, double *, int *, double *, double *, int *)
chemm
int BLASFUNC() chemm(char *, char *, int *, int *, float *, float *, int *, float *, int *, float *, float *, int *)
zhemm
int BLASFUNC() zhemm(char *, char *, int *, int *, double *, double *, int *, double *, int *, double *, double *, int *)
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
Eigen
Namespace containing all symbols from the Eigen library.
Definition: bench_norm.cpp:70
Eigen::dcomplex
std::complex< double > dcomplex
Definition: MKL_support.h:128
Eigen::scomplex
std::complex< float > scomplex
Definition: MKL_support.h:129
internal
Definition: Eigen_Colamd.h:49