Eigen  5.0.1
 
Loading...
Searching...
No Matches
GeneralMatrixMatrix_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 * General matrix-matrix product functionality based on ?GEMM.
30 ********************************************************************************
31*/
32// SPDX-License-Identifier: BSD-3-Clause
33
34#ifndef EIGEN_GENERAL_MATRIX_MATRIX_BLAS_H
35#define EIGEN_GENERAL_MATRIX_MATRIX_BLAS_H
36
37// IWYU pragma: private
38#include "../InternalHeaderCheck.h"
39
40namespace Eigen {
41
42namespace internal {
43
44/**********************************************************************
45 * This file implements general matrix-matrix multiplication using BLAS
46 * gemm function via partial specialization of
47 * general_matrix_matrix_product::run(..) method for float, double,
48 * std::complex<float> and std::complex<double> types
49 **********************************************************************/
50
51// gemm specialization
52
53#define EIGEN_BLAS_GEMM_SPECIALIZATION(EIGTYPE, EIGPREFIX, BLASTYPE, BLASFUNC) \
54 template <typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> \
55 struct general_matrix_matrix_product<Index, EIGTYPE, LhsStorageOrder, ConjugateLhs, EIGTYPE, RhsStorageOrder, \
56 ConjugateRhs, ColMajor, 1> { \
57 typedef gebp_traits<EIGTYPE, EIGTYPE> Traits; \
58 \
59 static void run(Index rows, Index cols, Index depth, const EIGTYPE* lhs_, Index lhsStride, const EIGTYPE* rhs_, \
60 Index rhsStride, EIGTYPE* res, Index resIncr, Index resStride, EIGTYPE alpha, \
61 level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/, GemmParallelInfo<Index>* /*info = 0*/) { \
62 using std::conj; \
63 if (rows == 0 || cols == 0 || depth == 0) return; \
64 EIGEN_ONLY_USED_FOR_DEBUG(resIncr); \
65 eigen_assert(resIncr == 1); \
66 char transa, transb; \
67 BlasIndex m, n, k, lda, ldb, ldc; \
68 const EIGTYPE *a, *b; \
69 EIGTYPE beta(1); \
70 MatrixX##EIGPREFIX a_tmp, b_tmp; \
71 \
72 /* Set transpose options */ \
73 transa = (LhsStorageOrder == RowMajor) ? ((ConjugateLhs) ? 'C' : 'T') : 'N'; \
74 transb = (RhsStorageOrder == RowMajor) ? ((ConjugateRhs) ? 'C' : 'T') : 'N'; \
75 \
76 /* Set m, n, k */ \
77 m = convert_index<BlasIndex>(rows); \
78 n = convert_index<BlasIndex>(cols); \
79 k = convert_index<BlasIndex>(depth); \
80 \
81 /* Set lda, ldb, ldc */ \
82 lda = convert_index<BlasIndex>(lhsStride); \
83 ldb = convert_index<BlasIndex>(rhsStride); \
84 ldc = convert_index<BlasIndex>(resStride); \
85 \
86 /* Set a, b, c */ \
87 EIGEN_IF_CONSTEXPR ((LhsStorageOrder == ColMajor) && (ConjugateLhs)) { \
88 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(lhs_, m, k, OuterStride<>(lhsStride)); \
89 a_tmp = lhs.conjugate(); \
90 a = a_tmp.data(); \
91 lda = convert_index<BlasIndex>(a_tmp.outerStride()); \
92 } else \
93 a = lhs_; \
94 \
95 EIGEN_IF_CONSTEXPR ((RhsStorageOrder == ColMajor) && (ConjugateRhs)) { \
96 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(rhs_, k, n, OuterStride<>(rhsStride)); \
97 b_tmp = rhs.conjugate(); \
98 b = b_tmp.data(); \
99 ldb = convert_index<BlasIndex>(b_tmp.outerStride()); \
100 } else \
101 b = rhs_; \
102 \
103 BLASFUNC(&transa, &transb, &m, &n, &k, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, &lda, \
104 (const BLASTYPE*)b, &ldb, (const BLASTYPE*)&numext::real_ref(beta), (BLASTYPE*)res, &ldc); \
105 } \
106 };
107
108#ifdef EIGEN_USE_MKL
109EIGEN_BLAS_GEMM_SPECIALIZATION(double, d, double, dgemm)
110EIGEN_BLAS_GEMM_SPECIALIZATION(float, f, float, sgemm)
111EIGEN_BLAS_GEMM_SPECIALIZATION(dcomplex, cd, MKL_Complex16, zgemm)
112EIGEN_BLAS_GEMM_SPECIALIZATION(scomplex, cf, MKL_Complex8, cgemm)
113#else
114EIGEN_BLAS_GEMM_SPECIALIZATION(double, d, double, EIGEN_BLAS_SYM(dgemm))
115EIGEN_BLAS_GEMM_SPECIALIZATION(float, f, float, EIGEN_BLAS_SYM(sgemm))
116EIGEN_BLAS_GEMM_SPECIALIZATION(dcomplex, cd, double, EIGEN_BLAS_SYM(zgemm))
117EIGEN_BLAS_GEMM_SPECIALIZATION(scomplex, cf, float, EIGEN_BLAS_SYM(cgemm))
118#endif
119
120// If OpenBLAS with BUILD_BFLOAT16=1 support is available,
121// use sbgemm for bfloat16.
122#if EIGEN_USE_OPENBLAS_BFLOAT16
123
124extern "C" {
125// OpenBLAS prototype.
126void EIGEN_BLAS_SYM(sbgemm)(const char* trans_a, const char* trans_b, const EIGEN_BLAS_INT* M, const EIGEN_BLAS_INT* N,
127 const EIGEN_BLAS_INT* K, const float* alpha, const Eigen::bfloat16* A,
128 const EIGEN_BLAS_INT* lda, const Eigen::bfloat16* B, const EIGEN_BLAS_INT* ldb,
129 const float* beta, float* C, const EIGEN_BLAS_INT* ldc);
130} // extern "C"
131
132template <typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs>
133struct general_matrix_matrix_product<Index, Eigen::bfloat16, LhsStorageOrder, ConjugateLhs, Eigen::bfloat16,
134 RhsStorageOrder, ConjugateRhs, ColMajor, 1> {
135 typedef gebp_traits<Eigen::bfloat16, Eigen::bfloat16> Traits;
136
137 static void run(Index rows, Index cols, Index depth, const Eigen::bfloat16* lhs_, Index lhsStride,
138 const Eigen::bfloat16* rhs_, Index rhsStride, Eigen::bfloat16* res, Index resIncr, Index resStride,
139 Eigen::bfloat16 alpha, level3_blocking<Eigen::bfloat16, Eigen::bfloat16>& /*blocking*/,
140 GemmParallelInfo<Index>* /*info = 0*/) {
141 using std::conj;
142 if (rows == 0 || cols == 0 || depth == 0) return;
143 EIGEN_ONLY_USED_FOR_DEBUG(resIncr);
144 eigen_assert(resIncr == 1);
145 char transa, transb;
146 BlasIndex m, n, k, lda, ldb, ldc;
147 const Eigen::bfloat16 *a, *b;
148
149 float falpha = static_cast<float>(alpha);
150 float fbeta = float(1.0);
151
152 using MatrixXbf16 = Matrix<Eigen::bfloat16, Dynamic, Dynamic>;
153 MatrixXbf16 a_tmp, b_tmp;
154 MatrixXf r_tmp;
155
156 /* Set transpose options */
157 transa = (LhsStorageOrder == RowMajor) ? ((ConjugateLhs) ? 'C' : 'T') : 'N';
158 transb = (RhsStorageOrder == RowMajor) ? ((ConjugateRhs) ? 'C' : 'T') : 'N';
159
160 /* Set m, n, k */
161 m = convert_index<BlasIndex>(rows);
162 n = convert_index<BlasIndex>(cols);
163 k = convert_index<BlasIndex>(depth);
164
165 /* Set lda, ldb, ldc */
166 lda = convert_index<BlasIndex>(lhsStride);
167 ldb = convert_index<BlasIndex>(rhsStride);
168 ldc = convert_index<BlasIndex>(m);
169
170 /* Set a, b, c */
171 EIGEN_IF_CONSTEXPR ((LhsStorageOrder == ColMajor) && (ConjugateLhs)) {
172 Map<const MatrixXbf16, 0, OuterStride<> > lhs(lhs_, m, k, OuterStride<>(lhsStride));
173 a_tmp = lhs.conjugate();
174 a = a_tmp.data();
175 lda = convert_index<BlasIndex>(a_tmp.outerStride());
176 } else {
177 a = lhs_;
178 }
179
180 EIGEN_IF_CONSTEXPR ((RhsStorageOrder == ColMajor) && (ConjugateRhs)) {
181 Map<const MatrixXbf16, 0, OuterStride<> > rhs(rhs_, k, n, OuterStride<>(rhsStride));
182 b_tmp = rhs.conjugate();
183 b = b_tmp.data();
184 ldb = convert_index<BlasIndex>(b_tmp.outerStride());
185 } else {
186 b = rhs_;
187 }
188
189 // Evaluate to a temporary intermediate array.
190 r_tmp.resize(m, n);
191
192 EIGEN_BLAS_SYM(sbgemm)
193 (&transa, &transb, &m, &n, &k, (const float*)&numext::real_ref(falpha), a, &lda, b, &ldb,
194 (const float*)&numext::real_ref(fbeta), r_tmp.data(), &ldc);
195
196 // Cast to the output.
197 Map<MatrixXbf16, 0, OuterStride<> > result(res, m, n, OuterStride<>(resStride));
198 result = r_tmp.cast<Eigen::bfloat16>();
199 }
200};
201
202#endif // EIGEN_USE_OPENBLAS_BFLOAT16
203
204#undef EIGEN_BLAS_GEMM_SPECIALIZATION
205} // namespace internal
206
207} // end namespace Eigen
208
209#endif // EIGEN_GENERAL_MATRIX_MATRIX_BLAS_H
Eigen::ColMajor
@ ColMajor
Definition Constants.h:319
Eigen::RowMajor
@ RowMajor
Definition Constants.h:321
Eigen::MatrixXf
Matrix< float, Dynamic, Dynamic > MatrixXf
Dynamic×Dynamic matrix of type float.
Definition Matrix.h:488