Eigen  5.0.1-dev+60122df6
 
Loading...
Searching...
No Matches
SelfadjointMatrixMatrix_BLAS.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 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
39namespace Eigen {
40
41namespace 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
149EIGEN_BLAS_SYMM_L(double, double, d, dsymm)
150EIGEN_BLAS_SYMM_L(float, float, f, ssymm)
151EIGEN_BLAS_HEMM_L(dcomplex, MKL_Complex16, cd, zhemm)
152EIGEN_BLAS_HEMM_L(scomplex, MKL_Complex8, cf, chemm)
153#else
154EIGEN_BLAS_SYMM_L(double, double, d, dsymm_)
155EIGEN_BLAS_SYMM_L(float, float, f, ssymm_)
156EIGEN_BLAS_HEMM_L(dcomplex, double, cd, zhemm_)
157EIGEN_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
263EIGEN_BLAS_SYMM_R(double, double, d, dsymm)
264EIGEN_BLAS_SYMM_R(float, float, f, ssymm)
265EIGEN_BLAS_HEMM_R(dcomplex, MKL_Complex16, cd, zhemm)
266EIGEN_BLAS_HEMM_R(scomplex, MKL_Complex8, cf, chemm)
267#else
268EIGEN_BLAS_SYMM_R(double, double, d, dsymm_)
269EIGEN_BLAS_SYMM_R(float, float, f, ssymm_)
270EIGEN_BLAS_HEMM_R(dcomplex, double, cd, zhemm_)
271EIGEN_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
Namespace containing all symbols from the Eigen library.
Definition B01_Experimental.dox:1