Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Declares the interfaces for the BLAS used.
6 : ///
7 : /// Wrappers are defined to perform casts from different integer types
8 : /// when the natural type in C++ differs from the BLAS argument.
9 :
10 : #pragma once
11 :
12 : #include <complex>
13 :
14 : #ifndef SPECTRE_DEBUG
15 : #include <libxsmm.h>
16 : #endif // ifndef SPECTRE_DEBUG
17 : #include <gsl/gsl_cblas.h>
18 :
19 : #include "Utilities/ErrorHandling/Assert.hpp"
20 : #include "Utilities/Gsl.hpp"
21 :
22 : namespace blas_detail {
23 : extern "C" {
24 : double ddot_(const int& N, const double* X, const int& INCX, const double* Y,
25 : const int& INCY);
26 :
27 : // The final two arguments are the "hidden" lengths of the first two.
28 : // https://gcc.gnu.org/onlinedocs/gfortran/Argument-passing-conventions.html
29 : void dgemm_(const char& TRANSA, const char& TRANSB, const int& M, const int& N,
30 : const int& K, const double& ALPHA, const double* A, const int& LDA,
31 : const double* B, const int& LDB, const double& BETA,
32 : const double* C, const int& LDC, size_t, size_t);
33 : void zgemm_(const char& TRANSA, const char& TRANSB, const int& M, const int& N,
34 : const int& K, const std::complex<double>& ALPHA,
35 : const std::complex<double>* A, const int& LDA,
36 : const std::complex<double>* B, const int& LDB,
37 : const std::complex<double>& BETA, const std::complex<double>* C,
38 : const int& LDC, size_t, size_t);
39 :
40 : // The final argument is the "hidden" length of the first one.
41 : // https://gcc.gnu.org/onlinedocs/gfortran/Argument-passing-conventions.html
42 : void dgemv_(const char& TRANS, const int& M, const int& N, const double& ALPHA,
43 : const double* A, const int& LDA, const double* X, const int& INCX,
44 : const double& BETA, double* Y, const int& INCY, size_t);
45 : } // extern "C"
46 : } // namespace blas_detail
47 :
48 : /*!
49 : * \brief Disable OpenBLAS multithreading since it conflicts with Charm++
50 : * parallelism
51 : *
52 : * Add this function to the `charm_init_node_funcs` of any executable that uses
53 : * BLAS routines.
54 : *
55 : * Details: https://github.com/xianyi/OpenBLAS/wiki/Faq#multi-threaded
56 : */
57 1 : void disable_openblas_multithreading();
58 :
59 : /// @{
60 : /*!
61 : * \ingroup UtilitiesGroup
62 : * The dot product of two vectors.
63 : *
64 : * \param N the length of the vectors.
65 : * \param X a pointer to the first element of the first vector.
66 : * \param INCX the stride for the elements of the first vector.
67 : * \param Y a pointer to the first element of the second vector.
68 : * \param INCY the stride for the elements of the second vector.
69 : * \return the dot product of the given vectors.
70 : */
71 1 : inline double ddot_(const size_t& N, const double* X, const size_t& INCX,
72 : const double* Y, const size_t& INCY) {
73 : // INCX and INCY are allowed to be negative by BLAS, but we never
74 : // use them that way. If needed, they can be changed here, but then
75 : // code providing values will also have to be changed to int to
76 : // avoid warnings.
77 : return blas_detail::ddot_(gsl::narrow_cast<int>(N), X,
78 : gsl::narrow_cast<int>(INCX), Y,
79 : gsl::narrow_cast<int>(INCY));
80 : }
81 : /// The unconjugated complex dot product $x \cdot y$. See `zdotc_` for the
82 : /// conjugated complex dot product, which is the standard dot product on the
83 : /// vector space of complex numbers.
84 1 : inline std::complex<double> zdotu_(const size_t& N,
85 : const std::complex<double>* X,
86 : const size_t& INCX,
87 : const std::complex<double>* Y,
88 : const size_t& INCY) {
89 : // The complex result of the BLAS zdot* functions is sometimes returned by
90 : // value and sometimes returned by reference, depending on the Fortran
91 : // compiler settings. By using the cblas interface we ensure a consistent
92 : // behavior.
93 : std::complex<double> result;
94 : cblas_zdotu_sub(gsl::narrow_cast<int>(N), X, gsl::narrow_cast<int>(INCX), Y,
95 : gsl::narrow_cast<int>(INCY), &result);
96 : return result;
97 : }
98 : /// The conjugated complex dot product $\bar{x} \cdot y$. This is the standard
99 : /// dot product on the vector space of complex numbers.
100 1 : inline std::complex<double> zdotc_(const size_t& N,
101 : const std::complex<double>* X,
102 : const size_t& INCX,
103 : const std::complex<double>* Y,
104 : const size_t& INCY) {
105 : std::complex<double> result;
106 : cblas_zdotc_sub(gsl::narrow_cast<int>(N), X, gsl::narrow_cast<int>(INCX), Y,
107 : gsl::narrow_cast<int>(INCY), &result);
108 : return result;
109 : }
110 : /// @}
111 :
112 : /// @{
113 : /*!
114 : * \ingroup UtilitiesGroup
115 : * \brief Perform a matrix-matrix multiplication
116 : *
117 : * Perform the matrix-matrix multiplication
118 : * \f[
119 : * C = \alpha \mathrm{op}(A) \mathrm{op}(B) + \beta \mathrm{op}(C)
120 : * \f]
121 : *
122 : * where \f$\mathrm{op}(A)\f$ represents either \f$A\f$ or \f$A^{T}\f$
123 : * (transpose of \f$A\f$).
124 : *
125 : * LIBXSMM, which is much faster than BLAS for small matrices, can be called
126 : * instead of BLAS by passing the template parameter `true`.
127 : *
128 : * \param TRANSA either 'N', 'T' or 'C', transposition of matrix A
129 : * \param TRANSB either 'N', 'T' or 'C', transposition of matrix B
130 : * \param M Number of rows in \f$\mathrm{op}(A)\f$
131 : * \param N Number of columns in \f$\mathrm{op}(B)\f$ and \f$\mathrm{op}(C)\f$
132 : * \param K Number of columns in \f$\mathrm{op}(A)\f$
133 : * \param ALPHA specifies \f$\alpha\f$
134 : * \param A Matrix \f$A\f$
135 : * \param LDA Specifies first dimension of \f$\mathrm{op}(A)\f$
136 : * \param B Matrix \f$B\f$
137 : * \param LDB Specifies first dimension of \f$\mathrm{op}(B)\f$
138 : * \param BETA specifies \f$\beta\f$
139 : * \param C Matrix \f$C\f$
140 : * \param LDC Specifies first dimension of \f$\mathrm{op}(C)\f$
141 : * \tparam UseLibXsmm if `true` then use LIBXSMM
142 : */
143 : template <bool UseLibXsmm = false>
144 1 : inline void dgemm_(const char& TRANSA, const char& TRANSB, const size_t& M,
145 : const size_t& N, const size_t& K, const double& ALPHA,
146 : const double* A, const size_t& LDA, const double* B,
147 : const size_t& LDB, const double& BETA, double* C,
148 : const size_t& LDC) {
149 : ASSERT('N' == TRANSA or 'n' == TRANSA or 'T' == TRANSA or 't' == TRANSA or
150 : 'C' == TRANSA or 'c' == TRANSA,
151 : "TRANSA must be upper or lower case N, T, or C. See the BLAS "
152 : "documentation for help.");
153 : ASSERT('N' == TRANSB or 'n' == TRANSB or 'T' == TRANSB or 't' == TRANSB or
154 : 'C' == TRANSB or 'c' == TRANSB,
155 : "TRANSB must be upper or lower case N, T, or C. See the BLAS "
156 : "documentation for help.");
157 :
158 : // On some BLAS implementations (e.g. Accelerate.framework on macOS) a call
159 : // with a zero-sized dimension aborts instead of acting as a no-op. Treat
160 : // these cases explicitly so we behave consistently across platforms.
161 : if (M == 0 or N == 0 or K == 0) {
162 : return;
163 : }
164 :
165 : blas_detail::dgemm_(
166 : TRANSA, TRANSB, gsl::narrow_cast<int>(M), gsl::narrow_cast<int>(N),
167 : gsl::narrow_cast<int>(K), ALPHA, A, gsl::narrow_cast<int>(LDA), B,
168 : gsl::narrow_cast<int>(LDB), BETA, C, gsl::narrow_cast<int>(LDC), 1, 1);
169 : }
170 : template <bool UseLibXsmm = false>
171 1 : inline void zgemm_(const char& TRANSA, const char& TRANSB, const size_t& M,
172 : const size_t& N, const size_t& K,
173 : const std::complex<double>& ALPHA,
174 : const std::complex<double>* A, const size_t& LDA,
175 : const std::complex<double>* B, const size_t& LDB,
176 : const std::complex<double>& BETA, std::complex<double>* C,
177 : const size_t& LDC) {
178 : ASSERT('N' == TRANSA or 'n' == TRANSA or 'T' == TRANSA or 't' == TRANSA or
179 : 'C' == TRANSA or 'c' == TRANSA,
180 : "TRANSA must be upper or lower case N, T, or C. See the BLAS "
181 : "documentation for help.");
182 : ASSERT('N' == TRANSB or 'n' == TRANSB or 'T' == TRANSB or 't' == TRANSB or
183 : 'C' == TRANSB or 'c' == TRANSB,
184 : "TRANSB must be upper or lower case N, T, or C. See the BLAS "
185 : "documentation for help.");
186 :
187 : // On some BLAS implementations (e.g. Accelerate.framework on macOS) a call
188 : // with a zero-sized dimension aborts instead of acting as a no-op. Treat
189 : // these cases explicitly so we behave consistently across platforms.
190 : if (M == 0 or N == 0 or K == 0) {
191 : return;
192 : }
193 :
194 : blas_detail::zgemm_(
195 : TRANSA, TRANSB, gsl::narrow_cast<int>(M), gsl::narrow_cast<int>(N),
196 : gsl::narrow_cast<int>(K), ALPHA, A, gsl::narrow_cast<int>(LDA), B,
197 : gsl::narrow_cast<int>(LDB), BETA, C, gsl::narrow_cast<int>(LDC), 1, 1);
198 : }
199 :
200 : // libxsmm is disabled in DEBUG builds because backtraces (from, for
201 : // example, FPEs) do not work when the error occurs in libxsmm code.
202 : #ifndef SPECTRE_DEBUG
203 : template <>
204 1 : inline void dgemm_<true>(const char& TRANSA, const char& TRANSB,
205 : const size_t& M, const size_t& N, const size_t& K,
206 : const double& ALPHA, const double* A,
207 : const size_t& LDA, const double* B, const size_t& LDB,
208 : const double& BETA, double* C, const size_t& LDC) {
209 : ASSERT('N' == TRANSA or 'n' == TRANSA or 'T' == TRANSA or 't' == TRANSA or
210 : 'C' == TRANSA or 'c' == TRANSA,
211 : "TRANSA must be upper or lower case N, T, or C. See the BLAS "
212 : "documentation for help.");
213 : ASSERT('N' == TRANSB or 'n' == TRANSB or 'T' == TRANSB or 't' == TRANSB or
214 : 'C' == TRANSB or 'c' == TRANSB,
215 : "TRANSB must be upper or lower case N, T, or C. See the BLAS "
216 : "documentation for help.");
217 :
218 : // On some BLAS implementations (e.g. Accelerate.framework on macOS) a call
219 : // with a zero-sized dimension aborts instead of acting as a no-op. Treat
220 : // these cases explicitly so we behave consistently across platforms.
221 : if (M == 0 or N == 0 or K == 0) {
222 : return;
223 : }
224 :
225 : // LIBXSMM's GEMM kernels only implement ALPHA == 1 and BETA in {0, 1}. Up to
226 : // v1.16.1 it forwarded any other values to the BLAS library, but since
227 : // v2.0.0 it silently computes as if ALPHA == 1 and BETA == 1 instead, which
228 : // gives wrong results (e.g. Evolution/DgSubcell/Matrices.cpp passes
229 : // ALPHA == 2). Take the BLAS path for the values LIBXSMM doesn't handle.
230 : if (ALPHA != 1.0 or (BETA != 0.0 and BETA != 1.0)) {
231 : blas_detail::dgemm_(
232 : TRANSA, TRANSB, gsl::narrow_cast<int>(M), gsl::narrow_cast<int>(N),
233 : gsl::narrow_cast<int>(K), ALPHA, A, gsl::narrow_cast<int>(LDA), B,
234 : gsl::narrow_cast<int>(LDB), BETA, C, gsl::narrow_cast<int>(LDC), 1, 1);
235 : return;
236 : }
237 :
238 : const auto m = gsl::narrow_cast<int>(M);
239 : const auto n = gsl::narrow_cast<int>(N);
240 : const auto k = gsl::narrow_cast<int>(K);
241 : const auto lda = gsl::narrow_cast<int>(LDA);
242 : const auto ldb = gsl::narrow_cast<int>(LDB);
243 : const auto ldc = gsl::narrow_cast<int>(LDC);
244 : libxsmm_dgemm(&TRANSA, &TRANSB, &m, &n, &k, &ALPHA, A, &lda, B, &ldb, &BETA,
245 : C, &ldc);
246 : }
247 : #endif // ifndef SPECTRE_DEBUG
248 : /// @}
249 :
250 : /// @{
251 : /*!
252 : * \ingroup UtilitiesGroup
253 : * \brief Perform a matrix-vector multiplication
254 : *
255 : * \f[
256 : * y = \alpha \mathrm{op}(A) x + \beta y
257 : * \f]
258 : *
259 : * where \f$\mathrm{op}(A)\f$ represents either \f$A\f$ or \f$A^{T}\f$
260 : * (transpose of \f$A\f$).
261 : *
262 : * \param TRANS either 'N', 'T' or 'C', transposition of matrix A
263 : * \param M Number of rows in \f$\mathrm{op}(A)\f$
264 : * \param N Number of columns in \f$\mathrm{op}(A)\f$
265 : * \param ALPHA specifies \f$\alpha\f$
266 : * \param A Matrix \f$A\f$
267 : * \param LDA Specifies first dimension of \f$\mathrm{op}(A)\f$
268 : * \param X Vector \f$x\f$
269 : * \param INCX Specifies the increment for the elements of \f$x\f$
270 : * \param BETA Specifies \f$\beta\f$
271 : * \param Y Vector \f$y\f$
272 : * \param INCY Specifies the increment for the elements of \f$y\f$
273 : */
274 1 : inline void dgemv_(const char& TRANS, const size_t& M, const size_t& N,
275 : const double& ALPHA, const double* A, const size_t& LDA,
276 : const double* X, const size_t& INCX, const double& BETA,
277 : double* Y, const size_t& INCY) {
278 : ASSERT('N' == TRANS or 'n' == TRANS or 'T' == TRANS or 't' == TRANS or
279 : 'C' == TRANS or 'c' == TRANS,
280 : "TRANS must be upper or lower case N, T, or C. See the BLAS "
281 : "documentation for help.");
282 :
283 : // On some BLAS implementations (e.g. Accelerate.framework on macOS) a call
284 : // with a zero-sized dimension aborts instead of acting as a no-op. Treat
285 : // these cases explicitly so we behave consistently across platforms.
286 : if (M == 0 or N == 0) {
287 : return;
288 : }
289 :
290 : // INCX and INCY are allowed to be negative by BLAS, but we never
291 : // use them that way. If needed, they can be changed here, but then
292 : // code providing values will also have to be changed to int to
293 : // avoid warnings.
294 : blas_detail::dgemv_(TRANS, gsl::narrow_cast<int>(M), gsl::narrow_cast<int>(N),
295 : ALPHA, A, gsl::narrow_cast<int>(LDA), X,
296 : gsl::narrow_cast<int>(INCX), BETA, Y,
297 : gsl::narrow_cast<int>(INCY), 1);
298 : }
299 : /// @}
|
