Blas.hpp
Go to the documentation of this file.
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 #ifndef SPECTRE_DEBUG
13 #include <libxsmm.h>
14 #endif // ifndef SPECTRE_DEBUG
15 
16 #include "ErrorHandling/Assert.hpp"
17 #include "Utilities/Gsl.hpp"
18 
19 namespace blas_detail {
20 extern "C" {
21 double ddot_(const int& N, const double* X, const int& INCX, const double* Y,
22  const int& INCY);
23 
24 void dgemm_(const char& TRANSA, const char& TRANSB, const int& M, const int& N,
25  const int& K, const double& ALPHA, const double* A, const int& LDA,
26  const double* B, const int& LDB, const double& BETA,
27  const double* C, const int& LDC);
28 
29 void dgemv_(const char& TRANS, const int& M, const int& N, const double& ALPHA,
30  const double* A, const int& LDA, const double* X, const int& INCX,
31  const double& BETA, double* Y, const int& INCY);
32 } // extern "C"
33 } // namespace blas_detail
34 
35 /*!
36  * \ingroup UtilitiesGroup
37  * The dot product of two vectors.
38  *
39  * \param N the length of the vectors.
40  * \param X a pointer to the first element of the first vector.
41  * \param INCX the stride for the elements of the first vector.
42  * \param Y a pointer to the first element of the second vector.
43  * \param INCY the stride for the elements of the second vector.
44  * \return the dot product of the given vectors.
45  */
46 inline double ddot_(const size_t& N, const double* X, const size_t& INCX,
47  const double* Y, const size_t& INCY) {
48  // INCX and INCY are allowed to be negative by BLAS, but we never
49  // use them that way. If needed, they can be changed here, but then
50  // code providing values will also have to be changed to int to
51  // avoid warnings.
52  return blas_detail::ddot_(gsl::narrow_cast<int>(N), X,
53  gsl::narrow_cast<int>(INCX), Y,
54  gsl::narrow_cast<int>(INCY));
55 }
56 
57 // @{
58 /*!
59  * \ingroup UtilitiesGroup
60  * \brief Perform a matrix-matrix multiplication
61  *
62  * Perform the matrix-matrix multiplication
63  * \f[
64  * C = \alpha \mathrm{op}(A) \mathrm{op}(B) + \beta \mathrm{op}(C)
65  * \f]
66  *
67  * where \f$\mathrm{op}(A)\f$ represents either \f$A\f$ or \f$A^{T}\f$
68  * (transpose of \f$A\f$).
69  *
70  * LIBXSMM, which is much faster than BLAS for small matrices, can be called
71  * instead of BLAS by passing the template parameter `true`.
72  *
73  * \param TRANSA either 'N', 'T' or 'C', transposition of matrix A
74  * \param TRANSB either 'N', 'T' or 'C', transposition of matrix B
75  * \param M Number of rows in \f$\mathrm{op}(A)\f$
76  * \param N Number of columns in \f$\mathrm{op}(B)\f$ and \f$\mathrm{op}(C)\f$
77  * \param K Number of columns in \f$\mathrm{op}(A)\f$
78  * \param ALPHA specifies \f$\alpha\f$
79  * \param A Matrix \f$A\f$
80  * \param LDA Specifies first dimension of \f$\mathrm{op}(A)\f$
81  * \param B Matrix \f$B\f$
82  * \param LDB Specifies first dimension of \f$\mathrm{op}(B)\f$
83  * \param BETA specifies \f$\beta\f$
84  * \param C Matrix \f$C\f$
85  * \param LDC Specifies first dimension of \f$\mathrm{op}(C)\f$
86  * \tparam UseLibXsmm if `true` then use LIBXSMM
87  */
88 template <bool UseLibXsmm = false>
89 inline void dgemm_(const char& TRANSA, const char& TRANSB, const size_t& M,
90  const size_t& N, const size_t& K, const double& ALPHA,
91  const double* A, const size_t& LDA, const double* B,
92  const size_t& LDB, const double& BETA, double* C,
93  const size_t& LDC) {
94  ASSERT('N' == TRANSA or 'n' == TRANSA or 'T' == TRANSA or 't' == TRANSA or
95  'C' == TRANSA or 'c' == TRANSA,
96  "TRANSA must be upper or lower case N, T, or C. See the BLAS "
97  "documentation for help.");
98  ASSERT('N' == TRANSB or 'n' == TRANSB or 'T' == TRANSB or 't' == TRANSB or
99  'C' == TRANSB or 'c' == TRANSB,
100  "TRANSB must be upper or lower case N, T, or C. See the BLAS "
101  "documentation for help.");
103  TRANSA, TRANSB, gsl::narrow_cast<int>(M), gsl::narrow_cast<int>(N),
104  gsl::narrow_cast<int>(K), ALPHA, A, gsl::narrow_cast<int>(LDA), B,
105  gsl::narrow_cast<int>(LDB), BETA, C, gsl::narrow_cast<int>(LDC));
106 }
107 
108 // libxsmm is disabled in DEBUG builds because backtraces (from, for
109 // example, FPEs) do not work when the error occurs in libxsmm code.
110 #ifndef SPECTRE_DEBUG
111 template <>
112 inline void dgemm_<true>(const char& TRANSA, const char& TRANSB,
113  const size_t& M, const size_t& N, const size_t& K,
114  const double& ALPHA, const double* A,
115  const size_t& LDA, const double* B, const size_t& LDB,
116  const double& BETA, double* C, const size_t& LDC) {
117  ASSERT('N' == TRANSA or 'n' == TRANSA or 'T' == TRANSA or 't' == TRANSA or
118  'C' == TRANSA or 'c' == TRANSA,
119  "TRANSA must be upper or lower case N, T, or C. See the BLAS "
120  "documentation for help.");
121  ASSERT('N' == TRANSB or 'n' == TRANSB or 'T' == TRANSB or 't' == TRANSB or
122  'C' == TRANSB or 'c' == TRANSB,
123  "TRANSB must be upper or lower case N, T, or C. See the BLAS "
124  "documentation for help.");
125  const auto m = gsl::narrow_cast<int>(M), n = gsl::narrow_cast<int>(N),
126  k = gsl::narrow_cast<int>(K), lda = gsl::narrow_cast<int>(LDA),
127  ldb = gsl::narrow_cast<int>(LDB), ldc = gsl::narrow_cast<int>(LDC);
128  libxsmm_dgemm(&TRANSA, &TRANSB, &m, &n, &k, &ALPHA, A, &lda, B, &ldb, &BETA,
129  C, &ldc);
130 }
131 #endif // ifndef SPECTRE_DEBUG
132 // @}
133 
134 // @{
135 /*!
136  * \ingroup UtilitiesGroup
137  * \brief Perform a matrix-vector multiplication
138  *
139  * \f[
140  * y = \alpha \mathrm{op}(A) x + \beta y
141  * \f]
142  *
143  * where \f$\mathrm{op}(A)\f$ represents either \f$A\f$ or \f$A^{T}\f$
144  * (transpose of \f$A\f$).
145  *
146  * \param TRANS either 'N', 'T' or 'C', transposition of matrix A
147  * \param M Number of rows in \f$\mathrm{op}(A)\f$
148  * \param N Number of columns in \f$\mathrm{op}(A)\f$
149  * \param ALPHA specifies \f$\alpha\f$
150  * \param A Matrix \f$A\f$
151  * \param LDA Specifies first dimension of \f$\mathrm{op}(A)\f$
152  * \param X Vector \f$x\f$
153  * \param INCX Specifies the increment for the elements of \f$x\f$
154  * \param BETA Specifies \f$\beta\f$
155  * \param Y Vector \f$y\f$
156  * \param INCY Specifies the increment for the elements of \f$y\f$
157  */
158 inline void dgemv_(const char& TRANS, const size_t& M, const size_t& N,
159  const double& ALPHA, const double* A, const size_t& LDA,
160  const double* X, const size_t& INCX, const double& BETA,
161  double* Y, const size_t& INCY) {
162  ASSERT('N' == TRANS or 'n' == TRANS or 'T' == TRANS or 't' == TRANS or
163  'C' == TRANS or 'c' == TRANS,
164  "TRANS must be upper or lower case N, T, or C. See the BLAS "
165  "documentation for help.");
166  // INCX and INCY are allowed to be negative by BLAS, but we never
167  // use them that way. If needed, they can be changed here, but then
168  // code providing values will also have to be changed to int to
169  // avoid warnings.
170  blas_detail::dgemv_(TRANS, gsl::narrow_cast<int>(M), gsl::narrow_cast<int>(N),
171  ALPHA, A, gsl::narrow_cast<int>(LDA), X,
172  gsl::narrow_cast<int>(INCX), BETA, Y,
173  gsl::narrow_cast<int>(INCY));
174 }
175 //@}
double ddot_(const size_t &N, const double *X, const size_t &INCX, const double *Y, const size_t &INCY)
Definition: Blas.hpp:46
void dgemm_< true >(const char &TRANSA, const char &TRANSB, const size_t &M, const size_t &N, const size_t &K, const double &ALPHA, const double *A, const size_t &LDA, const double *B, const size_t &LDB, const double &BETA, double *C, const size_t &LDC)
Perform a matrix-matrix multiplication.
Definition: Blas.hpp:112
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
constexpr T narrow_cast(U &&u) noexcept
Cast u to a type T where the cast may result in narrowing.
Definition: Gsl.hpp:86
void dgemm_(const char &TRANSA, const char &TRANSB, const size_t &M, const size_t &N, const size_t &K, const double &ALPHA, const double *A, const size_t &LDA, const double *B, const size_t &LDB, const double &BETA, double *C, const size_t &LDC)
Perform a matrix-matrix multiplication.
Definition: Blas.hpp:89
Defines macro ASSERT.
Defines functions and classes from the GSL.
void dgemv_(const char &TRANS, const size_t &M, const size_t &N, const double &ALPHA, const double *A, const size_t &LDA, const double *X, const size_t &INCX, const double &BETA, double *Y, const size_t &INCY)
Perform a matrix-vector multiplication.
Definition: Blas.hpp:158
Definition: Blas.hpp:19