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