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-2025 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 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) * ELEM_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 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 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
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
