Algol 68 Genie

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-2026 J. Marcel van der Veer [algol68g@algol68genie.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 a68g_vector_free (gsl_vector *V)
  33  {
  34    if (V != NO_REAL_VECTOR) {
  35      gsl_vector_free (V);
  36    }
  37  }
  38  
  39  void a68g_matrix_free (gsl_matrix *M)
  40  {
  41    if (M != NO_REAL_MATRIX) {
  42      gsl_matrix_free (M);
  43    }
  44  }
  45  
  46  void a68g_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