single-blas.c
1 //! @file single-blas.c
2 //! @author J. Marcel van der Veer
3
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2024 J. Marcel van der Veer [algol68g@xs4all.nl].
8
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! REAL GSL BLAS support.
25
26 #include "a68g.h"
27
28 #if defined (HAVE_GSL)
29
30 #include "a68g-torrix.h"
31
32 void a68_vector_free (gsl_vector *V)
33 {
34 if (V != NO_REAL_VECTOR) {
35 gsl_vector_free (V);
36 }
37 }
38
39 void a68_matrix_free (gsl_matrix *M)
40 {
41 if (M != NO_REAL_MATRIX) {
42 gsl_matrix_free (M);
43 }
44 }
45
46 void a68_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
47 double alpha, gsl_matrix *A, gsl_matrix *B,
48 double beta, gsl_matrix **C)
49 {
50 // Wrapper for gsl_blas_dgemm, allocates result matrix C if needed.
51 // GEMM from BLAS computes C := alpha * TransA (A) * TransB (B) + beta * C
52 if ((*C) == NO_REAL_MATRIX) {
53 unt N = (TransA == SELF ? SIZE1 (A) : SIZE2 (A));
54 unt M = (TransB == SELF ? SIZE2 (B) : SIZE1 (B));
55 (*C) = gsl_matrix_calloc (N, M); // NxM * MxP gives NxP.
56 }
57 ASSERT_GSL (gsl_blas_dgemm (TransA, TransB, alpha, A, B, beta, (*C)))
58 }
59
60 #endif
