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

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:04:03