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-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 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);
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 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);
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 (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
