## single-python.c

```
1  //! @file single-python.c
2  //! @author J. Marcel van der Veer
3
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
10  //!
11  //! This program is free software; you can redistribute it and/or modify it
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
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    unt M = SIZE1 (A), N = SIZE2 (A);
82    fprintf (stdout, "[%d, %d]\n", M, 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    A68_INT width;
104    POP_OBJECT (p, &width, A68_INT);
105    gsl_vector *V = pop_vector (p, A68_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    A68_INT width;
115    POP_OBJECT (p, &width, A68_INT);
116    gsl_matrix *M = pop_matrix (p, A68_TRUE);
117    print_matrix (M, VALUE (&width));
118    gsl_matrix_free (M);
119  }
120
121  //! @brief Convert VECTOR to [] REAL.
122
123  A68_ROW vector_to_row (NODE_T * p, gsl_vector * v)
124  {
125    unt len = (unt) (SIZE (v));
126    A68_ROW desc, row; A68_ARRAY arr; A68_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) * ELEM_SIZE (&arr);
131    for (int k = 0; k < len; k++, idx += inc) {
132      A68_REAL *x = (A68_REAL *) (base + idx);
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  A68_ROW matrix_to_row (NODE_T * p, gsl_matrix * a)
143  {
144    int len1 = (int) (SIZE1 (a)), len2 = (int) (SIZE2 (a));
145    A68_REF desc; A68_ROW row; A68_ARRAY arr; A68_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    MOID (&arr) = M_REAL;
150    ELEM_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) * ELEM_SIZE (&arr), inc2 = SPAN (&tup2) * ELEM_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        A68_REAL *x = (A68_REAL *) (base + idx2);
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 (A68_LEVEL >= 3)
177    DOUBLE_T sum = 0.0;
178    unt 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    unt 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, A68_TRUE);
203    PUSH_VALUE (p, matrix_norm (a), A68_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    unt Mv = SIZE1 (v), Nv = SIZE2 (v);
213    unt Mu, Nu;
214    if (u == NO_REAL_MATRIX) {
215      Mu = Nu = 0;
216    } else {
217      Mu = SIZE1 (u), Nu = SIZE2 (u);
218    }
219    if (Mu == 0 && Nu == 0) {
220      Mu = Mv;
221    } else {
222      MATH_RTE (p, Mu != Mv, M_INT, "number of rows does not match");
223    }
224    gsl_matrix *w = gsl_matrix_calloc (Mu, Nu + Nv);
225    for (int i = 0; i < Mu; i++) {
226      unt k = 0;
227      for (int j = 0; j < Nu; j++) {
228        gsl_matrix_set (w, i, k++, gsl_matrix_get (u, i, j));
229      }
230      for (int j = 0; j < Nv; j++) {
231        gsl_matrix_set (w, i, k++, gsl_matrix_get (v, i, j));
232      }
233    }
234    return (w);
235  }
236
237  void genie_matrix_hcat (NODE_T *p)
238  {
239  // Yield [u v].
240    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
241    gsl_matrix *v = pop_matrix (p, A68_TRUE);
242    gsl_matrix *u = pop_matrix (p, A68_TRUE);
243    gsl_matrix *w = matrix_hcat (p, u, v);
244    push_matrix (p, w);
245    gsl_matrix_free (u);
246    gsl_matrix_free (v);
247    gsl_matrix_free (w);
248    (void) gsl_set_error_handler (save_handler);
249  }
250
251  //! @brief OP VCAT = ([, ] REAL, [, ] REAL) [, ] REAL
252
253  gsl_matrix *matrix_vcat (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
254  {
255    unt Mv = SIZE1 (v), Nv = SIZE2 (v);
256    unt Mu, Nu;
257    if (u == NO_REAL_MATRIX) {
258      Mu = Nu = 0;
259    } else {
260      Mu = SIZE1 (u), Nu = SIZE2 (u);
261    }
262    if (Mu == 0 && Nu == 0) {
263      Nu = Nv;
264    } else {
265      MATH_RTE (p, Nu != Nv, M_INT, "number of columns does not match");
266    }
267    gsl_matrix *w = gsl_matrix_calloc (Mu + Mv, Nu);
268    for (int j = 0; j < Nu; j++) {
269      unt k = 0;
270      for (int i = 0; i < Mu; i++) {
271        gsl_matrix_set (w, k++, j, gsl_matrix_get (u, i, j));
272      }
273      for (int i = 0; i < Mv; i++) {
274        gsl_matrix_set (w, k++, j, gsl_matrix_get (v, i, j));
275      }
276    }
277    return (w);
278  }
279
280  void genie_matrix_vcat (NODE_T *p)
281  {
282  // Yield [u; v].
283    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
284    gsl_matrix *v = pop_matrix (p, A68_TRUE);
285    gsl_matrix *u = pop_matrix (p, A68_TRUE);
286    gsl_matrix *w = matrix_vcat (p, u, v);
287    push_matrix (p, w);
288    gsl_matrix_free (u);
289    gsl_matrix_free (v);
290    gsl_matrix_free (w);
291    (void) gsl_set_error_handler (save_handler);
292  }
293
294  gsl_matrix *mat_before_ab (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
295  {
296  // Form A := A BEFORE B.
297    gsl_matrix *w = matrix_hcat (p, u, v);
298    if (u != NO_REAL_MATRIX) {
299      a68_matrix_free (u); // Deallocate A, otherwise caller leaks memory.
300    }
301    return w;
302  }
303
304  gsl_matrix *mat_over_ab (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
305  {
306  // Form A := A OVER B.
307    gsl_matrix *w = matrix_vcat (p, u, v);
308    if (u != NO_REAL_MATRIX) {
309      a68_matrix_free (u); // Deallocate A, otherwise caller leaks memory.
310    }
311    return w;
312  }
313
314  #endif

```

 © 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)