Algol 68 Genie

single-python.c

     
   1  //! @file single-python.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 vector and matrix routines, Python look-a-likes.
  25  
  26  #include "a68g.h"
  27  
  28  #if defined (HAVE_GSL)
  29  
  30  #include "a68g-double.h"
  31  #include "a68g-genie.h"
  32  #include "a68g-prelude-gsl.h"
  33  #include "a68g-prelude.h"
  34  #include "a68g-torrix.h"
  35  
  36  //! Print REAL vector and matrix using GSL.
  37  
  38  #define FMT "%12.4g"
  39  
  40  void print_vector (gsl_vector *A, unt w)
  41  {
  42    unt N = SIZE (A);
  43    fprintf (stdout, "[%d]\n", N);
  44    if (w <= 0 || N <= 2 * w) {
  45      for (int i = 0; i < N; i++) {
  46        fprintf (stdout, FMT, gsl_vector_get (A, i));
  47      }
  48    } else {
  49      for (int i = 0; i < w; i++) {
  50        fprintf (stdout, FMT, gsl_vector_get (A, i));
  51      }
  52      fprintf (stdout, " ... ");
  53      for (int i = N - w; i < N; i++) {
  54        fprintf (stdout, FMT, gsl_vector_get (A, i));
  55      }
  56    }
  57    fprintf (stdout, "\n");
  58  }
  59  
  60  void print_row (gsl_matrix *A, unt m, unt w)
  61  {
  62    unt N = SIZE2 (A);
  63    if (w <= 0 || N <= 2 * w) {
  64      for (int i = 0; i < N; i++) {
  65        fprintf (stdout, FMT, gsl_matrix_get (A, m, i));
  66      }
  67    } else {
  68      for (int i = 0; i < w; i++) {
  69        fprintf (stdout, FMT, gsl_matrix_get (A, m, i));
  70      }
  71      fprintf (stdout, " ... ");
  72      for (int i = N - w; i < N; i++) {
  73        fprintf (stdout, FMT, gsl_matrix_get (A, m, i));
  74      }
  75    }
  76    fprintf (stdout, "\n");
  77  }
  78  
  79  void print_matrix (gsl_matrix *A, unt w)
  80  {
  81    size_t M = SIZE1 (A), N = SIZE2 (A);
  82    fprintf (stdout, "[" A68G_LU " " A68G_LU "]\n", (UNSIGNED_T) M, (UNSIGNED_T) N);
  83    if (w <= 0 || M <= 2 * w) {
  84      for (int i = 0; i < M; i++) {
  85        print_row (A, i, w);
  86      }
  87    } else {
  88      for (int i = 0; i < w; i++) {
  89        print_row (A, i, w);
  90      }
  91      fprintf (stdout, " ...\n");
  92      for (int i = M - w; i < M; i++) {
  93        print_row (A, i, w);
  94      }
  95    }
  96    fprintf (stdout, "\n");
  97  }
  98  
  99  //! @brief PROC print vector = ([] REAL v, INT width) VOID
 100  
 101  void genie_print_vector (NODE_T *p)
 102  {
 103    A68G_INT width;
 104    POP_OBJECT (p, &width, A68G_INT);
 105    gsl_vector *V = pop_vector (p, A68G_TRUE);
 106    print_vector (V, VALUE (&width));
 107    gsl_vector_free (V);
 108  }
 109  
 110  //! @brief PROC print matrix = ([, ] REAL v, INT width) VOID
 111  
 112  void genie_print_matrix (NODE_T *p)
 113  {
 114    A68G_INT width;
 115    POP_OBJECT (p, &width, A68G_INT);
 116    gsl_matrix *M = pop_matrix (p, A68G_TRUE);
 117    print_matrix (M, VALUE (&width));
 118    gsl_matrix_free (M);
 119  }
 120  
 121  //! @brief Convert VECTOR to [] REAL.
 122  
 123  A68G_ROW vector_to_row (NODE_T * p, gsl_vector * v)
 124  {
 125    unt len = (unt) (SIZE (v));
 126    A68G_ROW desc, row; A68G_ARRAY arr; A68G_TUPLE tup;
 127    NEW_ROW_1D (desc, row, arr, tup, M_ROW_REAL, M_REAL, len);
 128    BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
 129    int idx = VECTOR_OFFSET (&arr, &tup);
 130    unt inc = SPAN (&tup) * SLICE_SIZE (&arr);
 131    for (int k = 0; k < len; k++, idx += inc) {
 132      A68G_REAL *x = (A68G_REAL *) (base + idx);
 133      STATUS (x) = INIT_MASK;
 134      VALUE (x) = gsl_vector_get (v, (size_t) k);
 135      CHECK_REAL (p, VALUE (x));
 136    }
 137    return desc;
 138  }
 139  
 140  //! @brief Convert MATRIX to [, ] REAL.
 141  
 142  A68G_ROW matrix_to_row (NODE_T * p, gsl_matrix * a)
 143  {
 144    size_t len1 = SIZE1 (a), len2 = SIZE2 (a);
 145    A68G_REF desc; A68G_ROW row; A68G_ARRAY arr; A68G_TUPLE tup1, tup2;
 146    desc = heap_generator (p, M_ROW_ROW_REAL, DESCRIPTOR_SIZE (2));
 147    row = heap_generator_3 (p, M_ROW_ROW_REAL, len1, len2, SIZE (M_REAL));
 148    DIM (&arr) = 2;
 149    SLICE (&arr) = M_REAL;
 150    SLICE_SIZE (&arr) = SIZE (M_REAL);
 151    SLICE_OFFSET (&arr) = FIELD_OFFSET (&arr) = 0;
 152    ARRAY (&arr) = row;
 153    LWB (&tup1) = 1; UPB (&tup1) = len1; SPAN (&tup1) = 1;
 154    SHIFT (&tup1) = LWB (&tup1); K (&tup1) = 0;
 155    LWB (&tup2) = 1; UPB (&tup2) = len2; SPAN (&tup2) = ROW_SIZE (&tup1);
 156    SHIFT (&tup2) = LWB (&tup2) * SPAN (&tup2); K (&tup2) = 0;
 157    PUT_DESCRIPTOR2 (arr, tup1, tup2, &desc);
 158    BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
 159    int idx1 = MATRIX_OFFSET (&arr, &tup1, &tup2);
 160    int inc1 = SPAN (&tup1) * SLICE_SIZE (&arr), inc2 = SPAN (&tup2) * SLICE_SIZE (&arr);
 161    for (int k1 = 0; k1 < len1; k1++, idx1 += inc1) {
 162      for (int k2 = 0, idx2 = idx1; k2 < len2; k2++, idx2 += inc2) {
 163        A68G_REAL *x = (A68G_REAL *) (base + idx2);
 164        STATUS (x) = INIT_MASK;
 165        VALUE (x) = gsl_matrix_get (a, (size_t) k1, (size_t) k2);
 166        CHECK_REAL (p, VALUE (x));
 167      }
 168    }
 169    return desc;
 170  }
 171  
 172  //! @brief OP NORM = ([, ] REAL) REAL
 173  
 174  REAL_T matrix_norm (gsl_matrix *a)
 175  {
 176    #if (A68G_LEVEL >= 3)
 177      DOUBLE_T sum = 0.0;
 178      size_t M = SIZE1 (a), N = SIZE2 (a);
 179      for (int i = 0; i < M; i++) {
 180        for (int j = 0; j < N; j++) {
 181          DOUBLE_T z = gsl_matrix_get (a, i, j);
 182          sum += z * z;
 183        }
 184      }
 185      return ((REAL_T) sqrt_double (sum));
 186    #else
 187      REAL_T sum = 0.0;
 188      size_t M = SIZE1 (a), N = SIZE2 (a);
 189      for (int i = 0; i < M; i++) {
 190        for (int j = 0; j < N; j++) {
 191          REAL_T z = gsl_matrix_get (a, i, j);
 192          sum += z * z;
 193        }
 194      }
 195      return (sqrt (sum));
 196    #endif
 197  }
 198  
 199  void genie_matrix_norm (NODE_T * p)
 200  {
 201    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 202    gsl_matrix *a = pop_matrix (p, A68G_TRUE);
 203    PUSH_VALUE (p, matrix_norm (a), A68G_REAL);
 204    gsl_matrix_free (a);
 205    (void) gsl_set_error_handler (save_handler);
 206  }
 207  
 208  //! @brief OP HCAT = ([, ] REAL, [, ] REAL) [, ] REAL
 209  
 210  gsl_matrix *matrix_hcat (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
 211  {
 212    size_t Mv = SIZE1 (v), Nv = SIZE2 (v), Mu, Nu;
 213    if (u == NO_REAL_MATRIX) {
 214      Mu = Nu = 0;
 215    } else {
 216      Mu = SIZE1 (u), Nu = SIZE2 (u);
 217    }
 218    if (Mu == 0 && Nu == 0) {
 219      Mu = Mv;
 220    } else {
 221      MATH_RTE (p, Mu != Mv, M_INT, "number of rows does not match");
 222    }
 223    gsl_matrix *w = gsl_matrix_calloc (Mu, Nu + Nv);
 224    for (int i = 0; i < Mu; i++) {
 225      unt k = 0;
 226      for (int j = 0; j < Nu; j++) {
 227        gsl_matrix_set (w, i, k++, gsl_matrix_get (u, i, j));
 228      }
 229      for (int j = 0; j < Nv; j++) {
 230        gsl_matrix_set (w, i, k++, gsl_matrix_get (v, i, j));
 231      }
 232    }
 233    return (w);
 234  }
 235  
 236  void genie_matrix_hcat (NODE_T *p)
 237  {
 238  // Yield [u v].
 239    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 240    gsl_matrix *v = pop_matrix (p, A68G_TRUE);
 241    gsl_matrix *u = pop_matrix (p, A68G_TRUE);
 242    gsl_matrix *w = matrix_hcat (p, u, v);
 243    push_matrix (p, w);
 244    gsl_matrix_free (u);
 245    gsl_matrix_free (v);
 246    gsl_matrix_free (w);
 247    (void) gsl_set_error_handler (save_handler);
 248  }
 249  
 250  //! @brief OP VCAT = ([, ] REAL, [, ] REAL) [, ] REAL
 251  
 252  gsl_matrix *matrix_vcat (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
 253  {
 254    size_t Mv = SIZE1 (v), Nv = SIZE2 (v), Mu, Nu;
 255    if (u == NO_REAL_MATRIX) {
 256      Mu = Nu = 0;
 257    } else {
 258      Mu = SIZE1 (u), Nu = SIZE2 (u);
 259    }
 260    if (Mu == 0 && Nu == 0) {
 261      Nu = Nv;
 262    } else {  
 263      MATH_RTE (p, Nu != Nv, M_INT, "number of columns does not match");
 264    }
 265    gsl_matrix *w = gsl_matrix_calloc (Mu + Mv, Nu);
 266    for (int j = 0; j < Nu; j++) {
 267      unt k = 0;
 268      for (int i = 0; i < Mu; i++) {
 269        gsl_matrix_set (w, k++, j, gsl_matrix_get (u, i, j));
 270      }
 271      for (int i = 0; i < Mv; i++) {
 272        gsl_matrix_set (w, k++, j, gsl_matrix_get (v, i, j));
 273      }
 274    }
 275    return (w);
 276  }
 277  
 278  void genie_matrix_vcat (NODE_T *p)
 279  {
 280  // Yield [u; v].
 281    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 282    gsl_matrix *v = pop_matrix (p, A68G_TRUE);
 283    gsl_matrix *u = pop_matrix (p, A68G_TRUE);
 284    gsl_matrix *w = matrix_vcat (p, u, v);
 285    push_matrix (p, w);
 286    gsl_matrix_free (u);
 287    gsl_matrix_free (v);
 288    gsl_matrix_free (w);
 289    (void) gsl_set_error_handler (save_handler);
 290  }
 291  
 292  gsl_matrix *mat_before_ab (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
 293  {
 294  // Form A := A BEFORE B.
 295    gsl_matrix *w = matrix_hcat (p, u, v);
 296    if (u != NO_REAL_MATRIX) {
 297      a68g_matrix_free (u); // Deallocate A, otherwise caller leaks memory.
 298    }
 299    return w;
 300  }
 301  
 302  gsl_matrix *mat_over_ab (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
 303  {
 304  // Form A := A OVER B.
 305    gsl_matrix *w = matrix_vcat (p, u, v);
 306    if (u != NO_REAL_MATRIX) {
 307      a68g_matrix_free (u); // Deallocate A, otherwise caller leaks memory.
 308    }
 309    return w;
 310  }
 311  
 312  #endif