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

control_box_rst
Author(s): Christoph Rösmann
autogenerated on Wed Mar 2 2022 00:06:10