single-torrix-gsl.c
1 //! @file single-torrix-gsl.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 GSL vector and matrix routines.
25
26 #include "a68g.h"
27 #include "a68g-torrix.h"
28
29 #if defined (HAVE_GSL)
30
31 //! @brief Set permutation vector element - function fails in gsl.
32
33 void gsl_permutation_set (const gsl_permutation * p, const size_t i, const size_t j)
34 {
35 DATA (p)[i] = j;
36 }
37
38 //! @brief Map GSL error handler onto a68g error handler.
39
40 void torrix_error_handler (const char *reason, const char *file, int line, int gsl_errno)
41 {
42 if (line != 0) {
43 ASSERT (a68_bufprt (A68 (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
44 } else {
45 ASSERT (a68_bufprt (A68 (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
46 }
47 diagnostic (A68_RUNTIME_ERROR, A68 (f_entry), ERROR_TORRIX, A68 (edit_line), gsl_strerror (gsl_errno));
48 exit_genie (A68 (f_entry), A68_RUNTIME_ERROR);
49 }
50
51 //! @brief Pop [] INT on the stack as gsl_permutation.
52
53 gsl_permutation *pop_permutation (NODE_T * p, BOOL_T get)
54 {
55 A68_REF desc; A68_ARRAY *arr; A68_TUPLE *tup;
56 POP_REF (p, &desc);
57 CHECK_REF (p, desc, M_ROW_INT);
58 GET_DESCRIPTOR (arr, tup, &desc);
59 int len = ROW_SIZE (tup);
60 gsl_permutation *v = gsl_permutation_calloc ((size_t) len);
61 if (get && len > 0) {
62 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
63 int idx = VECTOR_OFFSET (arr, tup);
64 int inc = SPAN (tup) * ELEM_SIZE (arr);
65 for (int k = 0; k < len; k++, idx += inc) {
66 A68_INT *x = (A68_INT *) (base + idx);
67 CHECK_INIT (p, INITIALISED (x), M_INT);
68 gsl_permutation_set (v, (size_t) k, (size_t) VALUE (x));
69 }
70 }
71 return v;
72 }
73
74 //! @brief Push gsl_permutation on the stack as [] INT.
75
76 void push_permutation (NODE_T * p, gsl_permutation * v)
77 {
78 int len = (int) (SIZE (v));
79 A68_REF desc, row; A68_ARRAY arr; A68_TUPLE tup;
80 NEW_ROW_1D (desc, row, arr, tup, M_ROW_INT, M_INT, len);
81 BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
82 int idx = VECTOR_OFFSET (&arr, &tup);
83 int inc = SPAN (&tup) * ELEM_SIZE (&arr);
84 for (int k = 0; k < len; k++, idx += inc) {
85 A68_INT *x = (A68_INT *) (base + idx);
86 STATUS (x) = INIT_MASK;
87 VALUE (x) = (int) gsl_permutation_get (v, (size_t) k);
88 }
89 PUSH_REF (p, desc);
90 }
91
92 //! @brief Pop [] REAL on the stack as gsl_vector.
93
94 gsl_vector *pop_vector (NODE_T * p, BOOL_T get)
95 {
96 A68_REF desc; A68_ARRAY *arr; A68_TUPLE *tup;
97 POP_REF (p, &desc);
98 CHECK_REF (p, desc, M_ROW_REAL);
99 GET_DESCRIPTOR (arr, tup, &desc);
100 int len = ROW_SIZE (tup);
101 gsl_vector *v = gsl_vector_calloc ((size_t) len);
102 if (get && len > 0) {
103 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
104 int idx = VECTOR_OFFSET (arr, tup);
105 int inc = SPAN (tup) * ELEM_SIZE (arr);
106 for (int k = 0; k < len; k++, idx += inc) {
107 A68_REAL *x = (A68_REAL *) (base + idx);
108 CHECK_INIT (p, INITIALISED (x), M_REAL);
109 gsl_vector_set (v, (size_t) k, VALUE (x));
110 }
111 }
112 return v;
113 }
114
115 //! @brief Push gsl_vector on the stack as [] REAL.
116
117 void push_vector (NODE_T * p, gsl_vector * v)
118 {
119 PUSH_REF (p, vector_to_row (p, v));
120 }
121
122 //! @brief Pop [, ] REAL on the stack as gsl_matrix.
123
124 gsl_matrix *pop_matrix (NODE_T * p, BOOL_T get)
125 {
126 A68_REF desc; A68_ARRAY *arr; A68_TUPLE *tup1, *tup2;
127 POP_REF (p, &desc);
128 CHECK_REF (p, desc, M_ROW_ROW_REAL);
129 GET_DESCRIPTOR (arr, tup1, &desc);
130 tup2 = &(tup1[1]);
131 int len1 = ROW_SIZE (tup1), len2 = ROW_SIZE (tup2);
132 gsl_matrix *a = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
133 if (get && (len1 * len2 > 0)) {
134 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
135 int idx1 = MATRIX_OFFSET (arr, tup1, tup2);
136 int inc1 = SPAN (tup1) * ELEM_SIZE (arr), inc2 = SPAN (tup2) * ELEM_SIZE (arr);
137 for (int k1 = 0; k1 < len1; k1++, idx1 += inc1) {
138 for (int k2 = 0, idx2 = idx1; k2 < len2; k2++, idx2 += inc2) {
139 A68_REAL *x = (A68_REAL *) (base + idx2);
140 CHECK_INIT (p, INITIALISED (x), M_REAL);
141 gsl_matrix_set (a, (size_t) k1, (size_t) k2, VALUE (x));
142 }
143 }
144 }
145 return a;
146 }
147
148 //! @brief Push gsl_matrix on the stack as [, ] REAL.
149
150 void push_matrix (NODE_T * p, gsl_matrix * a)
151 {
152 PUSH_REF (p, matrix_to_row (p, a));
153 }
154
155 //! @brief Pop [] COMPLEX on the stack as gsl_vector_complex.
156
157 gsl_vector_complex *pop_vector_complex (NODE_T * p, BOOL_T get)
158 {
159 A68_REF desc; A68_ARRAY *arr; A68_TUPLE *tup;
160 POP_REF (p, &desc);
161 CHECK_REF (p, desc, M_ROW_COMPLEX);
162 GET_DESCRIPTOR (arr, tup, &desc);
163 int len = ROW_SIZE (tup);
164 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
165 if (get && len > 0) {
166 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
167 int idx = VECTOR_OFFSET (arr, tup);
168 int inc = SPAN (tup) * ELEM_SIZE (arr);
169 for (int k = 0; k < len; k++, idx += inc) {
170 A68_REAL *re = (A68_REAL *) (base + idx);
171 A68_REAL *im = (A68_REAL *) (base + idx + SIZE (M_REAL));
172 gsl_complex z;
173 CHECK_INIT (p, INITIALISED (re), M_COMPLEX);
174 CHECK_INIT (p, INITIALISED (im), M_COMPLEX);
175 GSL_SET_COMPLEX (&z, VALUE (re), VALUE (im));
176 gsl_vector_complex_set (v, (size_t) k, z);
177 }
178 }
179 return v;
180 }
181
182 //! @brief Push gsl_vector_complex on the stack as [] COMPLEX.
183
184 void push_vector_complex (NODE_T * p, gsl_vector_complex * v)
185 {
186 int len = (int) (SIZE (v));
187 A68_REF desc, row; A68_ARRAY arr; A68_TUPLE tup;
188 NEW_ROW_1D (desc, row, arr, tup, M_ROW_COMPLEX, M_COMPLEX, len);
189 BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
190 int idx = VECTOR_OFFSET (&arr, &tup);
191 int inc = SPAN (&tup) * ELEM_SIZE (&arr);
192 for (int k = 0; k < len; k++, idx += inc) {
193 A68_REAL *re = (A68_REAL *) (base + idx);
194 A68_REAL *im = (A68_REAL *) (base + idx + SIZE (M_REAL));
195 gsl_complex z = gsl_vector_complex_get (v, (size_t) k);
196 STATUS (re) = INIT_MASK;
197 VALUE (re) = GSL_REAL (z);
198 STATUS (im) = INIT_MASK;
199 VALUE (im) = GSL_IMAG (z);
200 CHECK_COMPLEX (p, VALUE (re), VALUE (im));
201 }
202 PUSH_REF (p, desc);
203 }
204
205 //! @brief Pop [, ] COMPLEX on the stack as gsl_matrix_complex.
206
207 gsl_matrix_complex *pop_matrix_complex (NODE_T * p, BOOL_T get)
208 {
209 A68_REF desc; A68_ARRAY *arr; A68_TUPLE *tup1, *tup2;
210 POP_REF (p, &desc);
211 CHECK_REF (p, desc, M_ROW_ROW_COMPLEX);
212 GET_DESCRIPTOR (arr, tup1, &desc);
213 tup2 = &(tup1[1]);
214 int len1 = ROW_SIZE (tup1), len2 = ROW_SIZE (tup2);
215 gsl_matrix_complex *a = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
216 if (get && (len1 * len2 > 0)) {
217 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
218 int idx1 = MATRIX_OFFSET (arr, tup1, tup2);
219 int inc1 = SPAN (tup1) * ELEM_SIZE (arr), inc2 = SPAN (tup2) * ELEM_SIZE (arr);
220 for (int k1 = 0; k1 < len1; k1++, idx1 += inc1) {
221 for (int k2 = 0, idx2 = idx1; k2 < len2; k2++, idx2 += inc2) {
222 A68_REAL *re = (A68_REAL *) (base + idx2);
223 A68_REAL *im = (A68_REAL *) (base + idx2 + SIZE (M_REAL));
224 CHECK_INIT (p, INITIALISED (re), M_COMPLEX);
225 CHECK_INIT (p, INITIALISED (im), M_COMPLEX);
226 gsl_complex z;
227 GSL_SET_COMPLEX (&z, VALUE (re), VALUE (im));
228 gsl_matrix_complex_set (a, (size_t) k1, (size_t) k2, z);
229 }
230 }
231 }
232 return a;
233 }
234
235 //! @brief Push gsl_matrix_complex on the stack as [, ] COMPLEX.
236
237 void push_matrix_complex (NODE_T * p, gsl_matrix_complex * a)
238 {
239 int len1 = (int) (SIZE1 (a)), len2 = (int) (SIZE2 (a));
240 A68_REF desc, row; A68_ARRAY arr; A68_TUPLE tup1, tup2;
241 desc = heap_generator (p, M_ROW_ROW_COMPLEX, DESCRIPTOR_SIZE (2));
242 row = heap_generator_3 (p, M_ROW_ROW_COMPLEX, len1, len2, 2 * SIZE (M_REAL));
243 DIM (&arr) = 2;
244 MOID (&arr) = M_COMPLEX;
245 ELEM_SIZE (&arr) = 2 * SIZE (M_REAL);
246 SLICE_OFFSET (&arr) = FIELD_OFFSET (&arr) = 0;
247 ARRAY (&arr) = row;
248 LWB (&tup1) = 1; UPB (&tup1) = len1; SPAN (&tup1) = 1;
249 SHIFT (&tup1) = LWB (&tup1); K (&tup1) = 0;
250 LWB (&tup2) = 1; UPB (&tup2) = len2; SPAN (&tup2) = ROW_SIZE (&tup1);
251 SHIFT (&tup2) = LWB (&tup2) * SPAN (&tup2); K (&tup2) = 0;
252 PUT_DESCRIPTOR2 (arr, tup1, tup2, &desc);
253 BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
254 int idx1 = MATRIX_OFFSET (&arr, &tup1, &tup2);
255 int inc1 = SPAN (&tup1) * ELEM_SIZE (&arr), inc2 = SPAN (&tup2) * ELEM_SIZE (&arr);
256 for (int k1 = 0; k1 < len1; k1++, idx1 += inc1) {
257 for (int k2 = 0, idx2 = idx1; k2 < len2; k2++, idx2 += inc2) {
258 A68_REAL *re = (A68_REAL *) (base + idx2);
259 A68_REAL *im = (A68_REAL *) (base + idx2 + SIZE (M_REAL));
260 gsl_complex z = gsl_matrix_complex_get (a, (size_t) k1, (size_t) k2);
261 STATUS (re) = INIT_MASK;
262 VALUE (re) = GSL_REAL (z);
263 STATUS (im) = INIT_MASK;
264 VALUE (im) = GSL_IMAG (z);
265 CHECK_COMPLEX (p, VALUE (re), VALUE (im));
266 }
267 }
268 PUSH_REF (p, desc);
269 }
270
271 //! @brief Generically perform operation and assign result (+:=, -:=, ...) .
272
273 void op_ab_torrix (NODE_T * p, MOID_T * m, MOID_T * n, GPROC * op)
274 {
275 ADDR_T parm_size = SIZE (m) + SIZE (n);
276 A68_REF dst, src, *save = (A68_REF *) STACK_OFFSET (-parm_size);
277 dst = *save;
278 CHECK_REF (p, dst, m);
279 *save = *DEREF (A68_ROW, &dst);
280 STATUS (&src) = (STATUS_MASK_T) (INIT_MASK | IN_STACK_MASK);
281 OFFSET (&src) = A68_SP - parm_size;
282 (*op) (p);
283 if (IS_REF (m)) {
284 genie_store (p, SUB (m), &dst, &src);
285 } else {
286 ABEND (A68_TRUE, ERROR_INTERNAL_CONSISTENCY, __func__);
287 }
288 *save = dst;
289 }
290
291 //! @brief PROC vector echo = ([] REAL) [] REAL
292
293 void genie_vector_echo (NODE_T * p)
294 {
295 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
296 gsl_vector *u = pop_vector (p, A68_TRUE);
297 push_vector (p, u);
298 gsl_vector_free (u);
299 (void) gsl_set_error_handler (save_handler);
300 }
301
302 //! @brief PROC matrix echo = ([, ] REAL) [, ] REAL
303
304 void genie_matrix_echo (NODE_T * p)
305 {
306 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
307 gsl_matrix *a = pop_matrix (p, A68_TRUE);
308 push_matrix (p, a);
309 gsl_matrix_free (a);
310 (void) gsl_set_error_handler (save_handler);
311 }
312
313 //! @brief PROC complex vector echo = ([] COMPLEX) [] COMPLEX
314
315 void genie_vector_complex_echo (NODE_T * p)
316 {
317 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
318 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
319 push_vector_complex (p, u);
320 gsl_vector_complex_free (u);
321 (void) gsl_set_error_handler (save_handler);
322 }
323
324 //! @brief PROC complex matrix echo = ([, ] COMPLEX) [, ] COMPLEX
325
326 void genie_matrix_complex_echo (NODE_T * p)
327 {
328 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
329 gsl_matrix_complex *a = pop_matrix_complex (p, A68_TRUE);
330 push_matrix_complex (p, a);
331 gsl_matrix_complex_free (a);
332 (void) gsl_set_error_handler (save_handler);
333 }
334
335 //! @brief OP ROW = ([] REAL) [, ] REAL
336
337 void genie_vector_row (NODE_T * p)
338 {
339 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
340 gsl_vector *u = pop_vector (p, A68_TRUE);
341 gsl_matrix *v = gsl_matrix_calloc (1, SIZE (u));
342 ASSERT_GSL (gsl_matrix_set_row (v, 0, u));
343 push_matrix (p, v);
344 gsl_vector_free (u);
345 gsl_matrix_free (v);
346 (void) gsl_set_error_handler (save_handler);
347 }
348
349 //! @brief OP COL = ([] REAL) [, ] REAL
350
351 void genie_vector_col (NODE_T * p)
352 {
353 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
354 gsl_vector *u = pop_vector (p, A68_TRUE);
355 gsl_matrix *v = gsl_matrix_calloc (SIZE (u), 1);
356 ASSERT_GSL (gsl_matrix_set_col (v, 0, u));
357 push_matrix (p, v);
358 gsl_vector_free (u);
359 gsl_matrix_free (v);
360 (void) gsl_set_error_handler (save_handler);
361 }
362
363 //! @brief OP - = ([] REAL) [] REAL
364
365 void genie_vector_minus (NODE_T * p)
366 {
367 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
368 gsl_vector *u = pop_vector (p, A68_TRUE);
369 ASSERT_GSL (gsl_vector_scale (u, -1));
370 push_vector (p, u);
371 gsl_vector_free (u);
372 (void) gsl_set_error_handler (save_handler);
373 }
374
375 //! @brief OP - = ([, ] REAL) [, ] REAL
376
377 void genie_matrix_minus (NODE_T * p)
378 {
379 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
380 gsl_matrix * a = pop_matrix (p, A68_TRUE);
381 ASSERT_GSL (gsl_matrix_scale (a, -1));
382 push_matrix (p, a);
383 gsl_matrix_free (a);
384 (void) gsl_set_error_handler (save_handler);
385 }
386
387 //! @brief OP T = ([, ] REAL) [, ] REAL
388
389 void genie_matrix_transpose (NODE_T * p)
390 {
391 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
392 gsl_matrix *a = pop_matrix (p, A68_TRUE);
393 gsl_matrix *t = gsl_matrix_calloc (SIZE2(a), SIZE1(a));
394 ASSERT_GSL (gsl_matrix_transpose_memcpy (t, a));
395 push_matrix (p, t);
396 gsl_matrix_free (a);
397 gsl_matrix_free (t);
398 (void) gsl_set_error_handler (save_handler);
399 }
400
401 //! @brief OP T = ([, ] COMPLEX) [, ] COMPLEX
402
403 void genie_matrix_complex_transpose (NODE_T * p)
404 {
405 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
406 gsl_matrix_complex *a = pop_matrix_complex (p, A68_TRUE);
407 gsl_matrix_complex *t = gsl_matrix_complex_calloc (SIZE2(a), SIZE1(a));
408 ASSERT_GSL ( gsl_matrix_complex_transpose_memcpy (t, a));
409 push_matrix_complex (p, a);
410 gsl_matrix_complex_free (a);
411 gsl_matrix_complex_free (t);
412 (void) gsl_set_error_handler (save_handler);
413 }
414
415 //! @brief OP INV = ([, ] REAL) [, ] REAL
416
417 void genie_matrix_inv (NODE_T * p)
418 {
419 // Avoid direct use of the inverse whenever possible; linear solver functions
420 // can obtain the same result more efficiently and reliably.
421 // Consult any introductory textbook on numerical linear algebra for details.
422 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
423 gsl_matrix *u = pop_matrix (p, A68_TRUE);
424 unt M = SIZE1 (u), N =SIZE2 (u);
425 MATH_RTE (p, M != N, M_ROW_ROW_REAL, "matrix is not square");
426 gsl_matrix *inv = NO_REAL_MATRIX;
427 compute_pseudo_inverse (p, &inv, u, 0); // Pseudo inverse equals inverse for a square matrix.
428 push_matrix (p, inv);
429 gsl_matrix_free (inv);
430 gsl_matrix_free (u);
431 (void) gsl_set_error_handler (save_handler);
432 }
433
434 //! @brief OP INV = ([, ] COMPLEX) [, ] COMPLEX
435
436 void genie_matrix_complex_inv (NODE_T * p)
437 {
438 // Avoid direct use of the inverse whenever possible; linear solver functions
439 // can obtain the same result more efficiently and reliably.
440 // Consult any introductory textbook on numerical linear algebra for details.
441 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
442 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
443 unt M = SIZE1 (u), N =SIZE2 (u);
444 MATH_RTE (p, M != N, M_ROW_ROW_COMPLEX, "matrix is not square");
445 gsl_permutation *q = gsl_permutation_calloc (SIZE1 (u));
446 int sign;
447 ASSERT_GSL (gsl_linalg_complex_LU_decomp (u, q, &sign));
448 gsl_matrix_complex *inv = gsl_matrix_complex_calloc (SIZE1 (u), SIZE2 (u));
449 ASSERT_GSL (gsl_linalg_complex_LU_invert (u, q, inv));
450 push_matrix_complex (p, inv);
451 gsl_matrix_complex_free (inv);
452 gsl_matrix_complex_free (u);
453 gsl_permutation_free (q);
454 (void) gsl_set_error_handler (save_handler);
455 }
456
457 //! @brief OP DET = ([, ] REAL) REAL
458
459 void genie_matrix_det (NODE_T * p)
460 {
461 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
462 gsl_matrix *u = pop_matrix (p, A68_TRUE);
463 gsl_permutation *q = gsl_permutation_calloc (SIZE1 (u));
464 int sign;
465 ASSERT_GSL (gsl_linalg_LU_decomp (u, q, &sign));
466 PUSH_VALUE (p, gsl_linalg_LU_det (u, sign), A68_REAL);
467 gsl_matrix_free (u);
468 gsl_permutation_free (q);
469 (void) gsl_set_error_handler (save_handler);
470 }
471
472 //! @brief OP DET = ([, ] COMPLEX) COMPLEX
473
474 void genie_matrix_complex_det (NODE_T * p)
475 {
476 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
477 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
478 gsl_permutation *q = gsl_permutation_calloc (SIZE1 (u));
479 int sign;
480 ASSERT_GSL (gsl_linalg_complex_LU_decomp (u, q, &sign));
481 gsl_complex det = gsl_linalg_complex_LU_det (u, sign);
482 PUSH_VALUE (p, GSL_REAL (det), A68_REAL);
483 PUSH_VALUE (p, GSL_IMAG (det), A68_REAL);
484 gsl_matrix_complex_free (u);
485 gsl_permutation_free (q);
486 (void) gsl_set_error_handler (save_handler);
487 }
488
489 //! @brief OP TRACE = ([, ] REAL) REAL
490
491 void genie_matrix_trace (NODE_T * p)
492 {
493 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
494 gsl_matrix *a = pop_matrix (p, A68_TRUE);
495 int len1 = (int) (SIZE1 (a)), len2 = (int) (SIZE2 (a));
496 if (len1 != len2) {
497 torrix_error_handler ("cannot calculate trace", __FILE__, __LINE__, GSL_ENOTSQR);
498 }
499 REAL_T sum = 0.0;
500 for (int k = 0; k < len1; k++) {
501 sum += gsl_matrix_get (a, (size_t) k, (size_t) k);
502 }
503 PUSH_VALUE (p, sum, A68_REAL);
504 gsl_matrix_free (a);
505 (void) gsl_set_error_handler (save_handler);
506 }
507
508 //! @brief OP TRACE = ([, ] COMPLEX) COMPLEX
509
510 void genie_matrix_complex_trace (NODE_T * p)
511 {
512 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
513 gsl_matrix_complex *a = pop_matrix_complex (p, A68_TRUE);
514 int len1 = (int) (SIZE1 (a)), len2 = (int) (SIZE2 (a));
515 if (len1 != len2) {
516 torrix_error_handler ("cannot calculate trace", __FILE__, __LINE__, GSL_ENOTSQR);
517 }
518 gsl_complex sum;
519 GSL_SET_COMPLEX (&sum, 0.0, 0.0);
520 for (int k = 0; k < len1; k++) {
521 sum = gsl_complex_add (sum, gsl_matrix_complex_get (a, (size_t) k, (size_t) k));
522 }
523 PUSH_VALUE (p, GSL_REAL (sum), A68_REAL);
524 PUSH_VALUE (p, GSL_IMAG (sum), A68_REAL);
525 gsl_matrix_complex_free (a);
526 (void) gsl_set_error_handler (save_handler);
527 }
528
529 //! @brief OP - = ([] COMPLEX) [] COMPLEX
530
531 void genie_vector_complex_minus (NODE_T * p)
532 {
533 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
534 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
535 gsl_blas_zdscal (-1, u);
536 push_vector_complex (p, u);
537 gsl_vector_complex_free (u);
538 (void) gsl_set_error_handler (save_handler);
539 }
540
541 //! @brief OP - = ([, ] COMPLEX) [, ] COMPLEX
542
543 void genie_matrix_complex_minus (NODE_T * p)
544 {
545 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
546 gsl_complex one;
547 GSL_SET_COMPLEX (&one, -1.0, 0.0);
548 gsl_matrix_complex *a = pop_matrix_complex (p, A68_TRUE);
549 ASSERT_GSL (gsl_matrix_complex_scale (a, one));
550 push_matrix_complex (p, a);
551 gsl_matrix_complex_free (a);
552 (void) gsl_set_error_handler (save_handler);
553 }
554
555 //! @brief OP + = ([] REAL, [] REAL) [] REAL
556
557 void genie_vector_add (NODE_T * p)
558 {
559 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
560 gsl_vector *v = pop_vector (p, A68_TRUE);
561 gsl_vector *u = pop_vector (p, A68_TRUE);
562 ASSERT_GSL (gsl_vector_add (u, v));
563 push_vector (p, u);
564 gsl_vector_free (u);
565 gsl_vector_free (v);
566 (void) gsl_set_error_handler (save_handler);
567 }
568
569 //! @brief OP - = ([] REAL, [] REAL) [] REAL
570
571 void genie_vector_sub (NODE_T * p)
572 {
573 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
574 gsl_vector *v = pop_vector (p, A68_TRUE);
575 gsl_vector *u = pop_vector (p, A68_TRUE);
576 ASSERT_GSL (gsl_vector_sub (u, v));
577 push_vector (p, u);
578 gsl_vector_free (u);
579 gsl_vector_free (v);
580 (void) gsl_set_error_handler (save_handler);
581 }
582
583 //! @brief OP = = ([] REAL, [] REAL) BOOL
584
585 void genie_vector_eq (NODE_T * p)
586 {
587 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
588 gsl_vector *v = pop_vector (p, A68_TRUE);
589 gsl_vector *u = pop_vector (p, A68_TRUE);
590 ASSERT_GSL (gsl_vector_sub (u, v));
591 PUSH_VALUE (p, (BOOL_T) (gsl_vector_isnull (u) ? A68_TRUE : A68_FALSE), A68_BOOL);
592 gsl_vector_free (u);
593 gsl_vector_free (v);
594 (void) gsl_set_error_handler (save_handler);
595 }
596
597 //! @brief OP /= = ([] REAL, [] REAL) BOOL
598
599 void genie_vector_ne (NODE_T * p)
600 {
601 genie_vector_eq (p);
602 genie_not_bool (p);
603 }
604
605 //! @brief OP +:= = (REF [] REAL, [] REAL) REF [] REAL
606
607 void genie_vector_plusab (NODE_T * p)
608 {
609 op_ab_torrix (p, M_REF_ROW_REAL, M_ROW_REAL, genie_vector_add);
610 }
611
612 //! @brief OP -:= = (REF [] REAL, [] REAL) REF [] REAL
613
614 void genie_vector_minusab (NODE_T * p)
615 {
616 op_ab_torrix (p, M_REF_ROW_REAL, M_ROW_REAL, genie_vector_sub);
617 }
618
619 //! @brief OP + = ([, ] REAL, [, ] REAL) [, ] REAL
620
621 void genie_matrix_add (NODE_T * p)
622 {
623 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
624 gsl_matrix *v = pop_matrix (p, A68_TRUE);
625 gsl_matrix *u = pop_matrix (p, A68_TRUE);
626 ASSERT_GSL (gsl_matrix_add (u, v));
627 push_matrix (p, u);
628 gsl_matrix_free (u);
629 gsl_matrix_free (v);
630 (void) gsl_set_error_handler (save_handler);
631 }
632
633 //! @brief OP - = ([, ] REAL, [, ] REAL) [, ] REAL
634
635 void genie_matrix_sub (NODE_T * p)
636 {
637 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
638 gsl_matrix *v = pop_matrix (p, A68_TRUE);
639 gsl_matrix *u = pop_matrix (p, A68_TRUE);
640 ASSERT_GSL (gsl_matrix_sub (u, v));
641 push_matrix (p, u);
642 gsl_matrix_free (u);
643 gsl_matrix_free (v);
644 (void) gsl_set_error_handler (save_handler);
645 }
646
647 //! @brief OP = = ([, ] REAL, [, ] REAL) [, ] BOOL
648
649 void genie_matrix_eq (NODE_T * p)
650 {
651 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
652 gsl_matrix *v = pop_matrix (p, A68_TRUE);
653 gsl_matrix *u = pop_matrix (p, A68_TRUE);
654 ASSERT_GSL (gsl_matrix_sub (u, v));
655 PUSH_VALUE (p, (BOOL_T) (gsl_matrix_isnull (u) ? A68_TRUE : A68_FALSE), A68_BOOL);
656 gsl_matrix_free (u);
657 gsl_matrix_free (v);
658 (void) gsl_set_error_handler (save_handler);
659 }
660
661 //! @brief OP /= = ([, ] REAL, [, ] REAL) [, ] BOOL
662
663 void genie_matrix_ne (NODE_T * p)
664 {
665 genie_matrix_eq (p);
666 genie_not_bool (p);
667 }
668
669 //! @brief OP +:= = (REF [, ] REAL, [, ] REAL) [, ] REAL
670
671 void genie_matrix_plusab (NODE_T * p)
672 {
673 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_ROW_ROW_REAL, genie_matrix_add);
674 }
675
676 //! @brief OP -:= = (REF [, ] REAL, [, ] REAL) [, ] REAL
677
678 void genie_matrix_minusab (NODE_T * p)
679 {
680 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_ROW_ROW_REAL, genie_matrix_sub);
681 }
682
683 //! @brief OP + = ([] COMPLEX, [] COMPLEX) [] COMPLEX
684
685 void genie_vector_complex_add (NODE_T * p)
686 {
687 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
688 gsl_complex one;
689 GSL_SET_COMPLEX (&one, 1.0, 0.0);
690 gsl_vector_complex *v = pop_vector_complex (p, A68_TRUE);
691 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
692 ASSERT_GSL (gsl_blas_zaxpy (one, u, v));
693 push_vector_complex (p, v);
694 gsl_vector_complex_free (u);
695 gsl_vector_complex_free (v);
696 (void) gsl_set_error_handler (save_handler);
697 }
698
699 //! @brief OP - = ([] COMPLEX, [] COMPLEX) [] COMPLEX
700
701 void genie_vector_complex_sub (NODE_T * p)
702 {
703 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
704 gsl_complex one;
705 GSL_SET_COMPLEX (&one, -1.0, 0.0);
706 gsl_vector_complex *v = pop_vector_complex (p, A68_TRUE);
707 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
708 ASSERT_GSL (gsl_blas_zaxpy (one, v, u));
709 push_vector_complex (p, u);
710 gsl_vector_complex_free (u);
711 gsl_vector_complex_free (v);
712 (void) gsl_set_error_handler (save_handler);
713 }
714
715 //! @brief OP = = ([] COMPLEX, [] COMPLEX) BOOL
716
717 void genie_vector_complex_eq (NODE_T * p)
718 {
719 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
720 gsl_complex one;
721 GSL_SET_COMPLEX (&one, -1.0, 0.0);
722 gsl_vector_complex *v = pop_vector_complex (p, A68_TRUE);
723 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
724 ASSERT_GSL (gsl_blas_zaxpy (one, v, u));
725 PUSH_VALUE (p, (BOOL_T) (gsl_vector_complex_isnull (u) ? A68_TRUE : A68_FALSE), A68_BOOL);
726 gsl_vector_complex_free (u);
727 gsl_vector_complex_free (v);
728 (void) gsl_set_error_handler (save_handler);
729 }
730
731 //! @brief OP /= = ([] COMPLEX, [] COMPLEX) BOOL
732
733 void genie_vector_complex_ne (NODE_T * p)
734 {
735 genie_vector_complex_eq (p);
736 genie_not_bool (p);
737 }
738
739 //! @brief OP +:= = (REF [] COMPLEX, [] COMPLEX) [] COMPLEX
740
741 void genie_vector_complex_plusab (NODE_T * p)
742 {
743 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_ROW_COMPLEX, genie_vector_complex_add);
744 }
745
746 //! @brief OP -:= = (REF [] COMPLEX, [] COMPLEX) [] COMPLEX
747
748 void genie_vector_complex_minusab (NODE_T * p)
749 {
750 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_ROW_COMPLEX, genie_vector_complex_sub);
751 }
752
753 //! @brief OP + = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
754
755 void genie_matrix_complex_add (NODE_T * p)
756 {
757 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
758 gsl_matrix_complex *v = pop_matrix_complex (p, A68_TRUE);
759 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
760 ASSERT_GSL (gsl_matrix_complex_add (u, v));
761 push_matrix_complex (p, u);
762 gsl_matrix_complex_free (u);
763 gsl_matrix_complex_free (v);
764 (void) gsl_set_error_handler (save_handler);
765 }
766
767 //! @brief OP - = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
768
769 void genie_matrix_complex_sub (NODE_T * p)
770 {
771 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
772 gsl_matrix_complex *v = pop_matrix_complex (p, A68_TRUE);
773 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
774 ASSERT_GSL (gsl_matrix_complex_sub (u, v));
775 push_matrix_complex (p, u);
776 gsl_matrix_complex_free (u);
777 gsl_matrix_complex_free (v);
778 (void) gsl_set_error_handler (save_handler);
779 }
780
781 //! @brief OP = = ([, ] COMPLEX, [, ] COMPLEX) BOOL
782
783 void genie_matrix_complex_eq (NODE_T * p)
784 {
785 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
786 gsl_matrix_complex *v = pop_matrix_complex (p, A68_TRUE);
787 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
788 ASSERT_GSL (gsl_matrix_complex_sub (u, v));
789 PUSH_VALUE (p, (BOOL_T) (gsl_matrix_complex_isnull (u) ? A68_TRUE : A68_FALSE), A68_BOOL);
790 gsl_matrix_complex_free (u);
791 gsl_matrix_complex_free (v);
792 (void) gsl_set_error_handler (save_handler);
793 }
794
795 //! @brief OP /= = ([, ] COMPLEX, [, ] COMPLEX) BOOL
796
797 void genie_matrix_complex_ne (NODE_T * p)
798 {
799 genie_matrix_complex_eq (p);
800 genie_not_bool (p);
801 }
802
803 //! @brief OP +:= = (REF [, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
804
805 void genie_matrix_complex_plusab (NODE_T * p)
806 {
807 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_ROW_ROW_COMPLEX, genie_matrix_complex_add);
808 }
809
810 //! @brief OP -:= = (REF [, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
811
812 void genie_matrix_complex_minusab (NODE_T * p)
813 {
814 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_ROW_ROW_COMPLEX, genie_matrix_complex_sub);
815 }
816
817 //! @brief OP * = ([] REAL, REAL) [] REAL
818
819 void genie_vector_scale_real (NODE_T * p)
820 {
821 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
822 A68_REAL v;
823 POP_OBJECT (p, &v, A68_REAL);
824 gsl_vector *u = pop_vector (p, A68_TRUE);
825 ASSERT_GSL (gsl_vector_scale (u, VALUE (&v)));
826 push_vector (p, u);
827 gsl_vector_free (u);
828 (void) gsl_set_error_handler (save_handler);
829 }
830
831 //! @brief OP * = (REAL, [] REAL) [] REAL
832
833 void genie_real_scale_vector (NODE_T * p)
834 {
835 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
836 gsl_vector *u = pop_vector (p, A68_TRUE);
837 A68_REAL v;
838 POP_OBJECT (p, &v, A68_REAL);
839 ASSERT_GSL (gsl_vector_scale (u, VALUE (&v)));
840 push_vector (p, u);
841 gsl_vector_free (u);
842 (void) gsl_set_error_handler (save_handler);
843 }
844
845 //! @brief OP * = ([, ] REAL, REAL) [, ] REAL
846
847 void genie_matrix_scale_real (NODE_T * p)
848 {
849 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
850 A68_REAL v;
851 POP_OBJECT (p, &v, A68_REAL);
852 gsl_matrix *u = pop_matrix (p, A68_TRUE);
853 ASSERT_GSL (gsl_matrix_scale (u, VALUE (&v)));
854 push_matrix (p, u);
855 gsl_matrix_free (u);
856 (void) gsl_set_error_handler (save_handler);
857 }
858
859 //! @brief OP * = (REAL, [, ] REAL) [, ] REAL
860
861 void genie_real_scale_matrix (NODE_T * p)
862 {
863 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
864 gsl_matrix *u = pop_matrix (p, A68_TRUE);
865 A68_REAL v;
866 POP_OBJECT (p, &v, A68_REAL);
867 ASSERT_GSL (gsl_matrix_scale (u, VALUE (&v)));
868 push_matrix (p, u);
869 gsl_matrix_free (u);
870 (void) gsl_set_error_handler (save_handler);
871 }
872
873 //! @brief OP * = ([] COMPLEX, COMPLEX) [] COMPLEX
874
875 void genie_vector_complex_scale_complex (NODE_T * p)
876 {
877 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
878 A68_REAL re, im; gsl_complex v;
879 POP_OBJECT (p, &im, A68_REAL);
880 POP_OBJECT (p, &re, A68_REAL);
881 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
882 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
883 gsl_blas_zscal (v, u);
884 push_vector_complex (p, u);
885 gsl_vector_complex_free (u);
886 (void) gsl_set_error_handler (save_handler);
887 }
888
889 //! @brief OP * = (COMPLEX, [] COMPLEX) [] COMPLEX
890
891 void genie_complex_scale_vector_complex (NODE_T * p)
892 {
893 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
894 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
895 A68_REAL re, im; gsl_complex v;
896 POP_OBJECT (p, &im, A68_REAL);
897 POP_OBJECT (p, &re, A68_REAL);
898 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
899 gsl_blas_zscal (v, u);
900 push_vector_complex (p, u);
901 gsl_vector_complex_free (u);
902 (void) gsl_set_error_handler (save_handler);
903 }
904
905 //! @brief OP * = ([, ] COMPLEX, COMPLEX) [, ] COMPLEX
906
907 void genie_matrix_complex_scale_complex (NODE_T * p)
908 {
909 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
910 A68_REAL re, im; gsl_complex v;
911 POP_OBJECT (p, &im, A68_REAL);
912 POP_OBJECT (p, &re, A68_REAL);
913 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
914 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
915 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
916 push_matrix_complex (p, u);
917 gsl_matrix_complex_free (u);
918 (void) gsl_set_error_handler (save_handler);
919 }
920
921 //! @brief OP * = (COMPLEX, [, ] COMPLEX) [, ] COMPLEX
922
923 void genie_complex_scale_matrix_complex (NODE_T * p)
924 {
925 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
926 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
927 A68_REAL re, im; gsl_complex v;
928 POP_OBJECT (p, &im, A68_REAL);
929 POP_OBJECT (p, &re, A68_REAL);
930 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
931 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
932 push_matrix_complex (p, u);
933 gsl_matrix_complex_free (u);
934 (void) gsl_set_error_handler (save_handler);
935 }
936
937 //! @brief OP *:= (REF [] REAL, REAL) REF [] REAL
938
939 void genie_vector_scale_real_ab (NODE_T * p)
940 {
941 op_ab_torrix (p, M_REF_ROW_REAL, M_REAL, genie_vector_scale_real);
942 }
943
944 //! @brief OP *:= (REF [, ] REAL, REAL) REF [, ] REAL
945
946 void genie_matrix_scale_real_ab (NODE_T * p)
947 {
948 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_REAL, genie_matrix_scale_real);
949 }
950
951 //! @brief OP *:= (REF [] COMPLEX, COMPLEX) REF [] COMPLEX
952
953 void genie_vector_complex_scale_complex_ab (NODE_T * p)
954 {
955 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_COMPLEX, genie_vector_complex_scale_complex);
956 }
957
958 //! @brief OP *:= (REF [, ] COMPLEX, COMPLEX) REF [, ] COMPLEX
959
960 void genie_matrix_complex_scale_complex_ab (NODE_T * p)
961 {
962 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_COMPLEX, genie_matrix_complex_scale_complex);
963 }
964
965 //! @brief OP / = ([] REAL, REAL) [] REAL
966
967 void genie_vector_div_real (NODE_T * p)
968 {
969 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
970 A68_REAL v;
971 POP_OBJECT (p, &v, A68_REAL);
972 if (VALUE (&v) == 0.0) {
973 diagnostic (A68_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_REAL);
974 exit_genie (p, A68_RUNTIME_ERROR);
975 }
976 gsl_vector *u = pop_vector (p, A68_TRUE);
977 ASSERT_GSL (gsl_vector_scale (u, 1.0 / VALUE (&v)));
978 push_vector (p, u);
979 gsl_vector_free (u);
980 (void) gsl_set_error_handler (save_handler);
981 }
982
983 //! @brief OP / = ([, ] REAL, REAL) [, ] REAL
984
985 void genie_matrix_div_real (NODE_T * p)
986 {
987 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
988 A68_REAL v;
989 POP_OBJECT (p, &v, A68_REAL);
990 if (VALUE (&v) == 0.0) {
991 diagnostic (A68_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_ROW_REAL);
992 exit_genie (p, A68_RUNTIME_ERROR);
993 }
994 gsl_matrix *u = pop_matrix (p, A68_TRUE);
995 ASSERT_GSL (gsl_matrix_scale (u, 1.0 / VALUE (&v)));
996 push_matrix (p, u);
997 gsl_matrix_free (u);
998 (void) gsl_set_error_handler (save_handler);
999 }
1000
1001 //! @brief OP / = ([] COMPLEX, COMPLEX) [] COMPLEX
1002
1003 void genie_vector_complex_div_complex (NODE_T * p)
1004 {
1005 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1006 A68_REAL re, im; gsl_complex v;
1007 POP_OBJECT (p, &im, A68_REAL);
1008 POP_OBJECT (p, &re, A68_REAL);
1009 if (VALUE (&re) == 0.0 && VALUE (&im) == 0.0) {
1010 diagnostic (A68_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_COMPLEX);
1011 exit_genie (p, A68_RUNTIME_ERROR);
1012 }
1013 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
1014 v = gsl_complex_inverse (v);
1015 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
1016 gsl_blas_zscal (v, u);
1017 push_vector_complex (p, u);
1018 gsl_vector_complex_free (u);
1019 (void) gsl_set_error_handler (save_handler);
1020 }
1021
1022 //! @brief OP / = ([, ] COMPLEX, COMPLEX) [, ] COMPLEX
1023
1024 void genie_matrix_complex_div_complex (NODE_T * p)
1025 {
1026 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1027 A68_REAL re, im; gsl_complex v;
1028 POP_OBJECT (p, &im, A68_REAL);
1029 POP_OBJECT (p, &re, A68_REAL);
1030 if (VALUE (&re) == 0.0 && VALUE (&im) == 0.0) {
1031 diagnostic (A68_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_ROW_COMPLEX);
1032 exit_genie (p, A68_RUNTIME_ERROR);
1033 }
1034 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
1035 v = gsl_complex_inverse (v);
1036 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
1037 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
1038 push_matrix_complex (p, u);
1039 gsl_matrix_complex_free (u);
1040 (void) gsl_set_error_handler (save_handler);
1041 }
1042
1043 //! @brief OP /:= (REF [] REAL, REAL) REF [] REAL
1044
1045 void genie_vector_div_real_ab (NODE_T * p)
1046 {
1047 op_ab_torrix (p, M_REF_ROW_REAL, M_REAL, genie_vector_div_real);
1048 }
1049
1050 //! @brief OP /:= (REF [, ] REAL, REAL) REF [, ] REAL
1051
1052 void genie_matrix_div_real_ab (NODE_T * p)
1053 {
1054 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_REAL, genie_matrix_div_real);
1055 }
1056
1057 //! @brief OP /:= (REF [] COMPLEX, COMPLEX) REF [] COMPLEX
1058
1059 void genie_vector_complex_div_complex_ab (NODE_T * p)
1060 {
1061 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_COMPLEX, genie_vector_complex_div_complex);
1062 }
1063
1064 //! @brief OP /:= (REF [, ] COMPLEX, COMPLEX) REF [, ] COMPLEX
1065
1066 void genie_matrix_complex_div_complex_ab (NODE_T * p)
1067 {
1068 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_COMPLEX, genie_matrix_complex_div_complex);
1069 }
1070
1071 //! @brief OP * = ([] REAL, [] REAL) REAL
1072
1073 void genie_vector_dot (NODE_T * p)
1074 {
1075 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1076 gsl_vector *v = pop_vector (p, A68_TRUE);
1077 gsl_vector *u = pop_vector (p, A68_TRUE);
1078 REAL_T w;
1079 ASSERT_GSL (gsl_blas_ddot (u, v, &w));
1080 PUSH_VALUE (p, w, A68_REAL);
1081 gsl_vector_free (u);
1082 gsl_vector_free (v);
1083 (void) gsl_set_error_handler (save_handler);
1084 }
1085
1086 //! @brief OP * = ([] COMPLEX, [] COMPLEX) COMPLEX
1087
1088 void genie_vector_complex_dot (NODE_T * p)
1089 {
1090 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1091 gsl_vector_complex *v = pop_vector_complex (p, A68_TRUE);
1092 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
1093 gsl_complex w;
1094 ASSERT_GSL (gsl_blas_zdotc (u, v, &w));
1095 PUSH_VALUE (p, GSL_REAL (w), A68_REAL);
1096 PUSH_VALUE (p, GSL_IMAG (w), A68_REAL);
1097 gsl_vector_complex_free (u);
1098 gsl_vector_complex_free (v);
1099 (void) gsl_set_error_handler (save_handler);
1100 }
1101
1102 //! @brief OP NORM = ([] REAL) REAL
1103
1104 void genie_vector_norm (NODE_T * p)
1105 {
1106 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1107 gsl_vector *u = pop_vector (p, A68_TRUE);
1108 PUSH_VALUE (p, gsl_blas_dnrm2 (u), A68_REAL);
1109 gsl_vector_free (u);
1110 (void) gsl_set_error_handler (save_handler);
1111 }
1112
1113 //! @brief OP NORM = ([] COMPLEX) COMPLEX
1114
1115 void genie_vector_complex_norm (NODE_T * p)
1116 {
1117 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1118 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
1119 PUSH_VALUE (p, gsl_blas_dznrm2 (u), A68_REAL);
1120 gsl_vector_complex_free (u);
1121 (void) gsl_set_error_handler (save_handler);
1122 }
1123
1124 //! @brief OP DYAD = ([] REAL, [] REAL) [, ] REAL
1125
1126 void genie_vector_dyad (NODE_T * p)
1127 {
1128 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1129 gsl_vector *v = pop_vector (p, A68_TRUE);
1130 gsl_vector *u = pop_vector (p, A68_TRUE);
1131 int len1 = (int) (SIZE (u)), len2 = (int) (SIZE (v));
1132 gsl_matrix *w = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
1133 for (int j = 0; j < len1; j++) {
1134 REAL_T uj = gsl_vector_get (u, (size_t) j);
1135 for (int k = 0; k < len2; k++) {
1136 REAL_T vk = gsl_vector_get (v, (size_t) k);
1137 gsl_matrix_set (w, (size_t) j, (size_t) k, uj * vk);
1138 }
1139 }
1140 push_matrix (p, w);
1141 gsl_vector_free (u);
1142 gsl_vector_free (v);
1143 gsl_matrix_free (w);
1144 (void) gsl_set_error_handler (save_handler);
1145 }
1146
1147 //! @brief OP DYAD = ([] COMPLEX, [] COMPLEX) [, ] COMPLEX
1148
1149 void genie_vector_complex_dyad (NODE_T * p)
1150 {
1151 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1152 gsl_vector_complex *v = pop_vector_complex (p, A68_TRUE);
1153 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
1154 int len1 = (int) (SIZE (u)), len2 = (int) (SIZE (v));
1155 gsl_matrix_complex *w = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
1156 for (int j = 0; j < len1; j++) {
1157 gsl_complex uj = gsl_vector_complex_get (u, (size_t) j);
1158 for (int k = 0; k < len2; k++) {
1159 gsl_complex vk = gsl_vector_complex_get (u, (size_t) k);
1160 gsl_matrix_complex_set (w, (size_t) j, (size_t) k, gsl_complex_mul (uj, vk));
1161 }
1162 }
1163 push_matrix_complex (p, w);
1164 gsl_vector_complex_free (u);
1165 gsl_vector_complex_free (v);
1166 gsl_matrix_complex_free (w);
1167 (void) gsl_set_error_handler (save_handler);
1168 }
1169
1170 //! @brief OP * = ([, ] REAL, [] REAL) [] REAL
1171
1172 void genie_matrix_times_vector (NODE_T * p)
1173 {
1174 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1175 gsl_vector *u = pop_vector (p, A68_TRUE);
1176 gsl_matrix *w = pop_matrix (p, A68_TRUE);
1177 int len = (int) (SIZE (u));
1178 gsl_vector *v = gsl_vector_calloc ((size_t) len);
1179 gsl_vector_set_zero (v);
1180 ASSERT_GSL (gsl_blas_dgemv (CblasNoTrans, 1.0, w, u, 0.0, v));
1181 push_vector (p, v);
1182 gsl_vector_free (u);
1183 gsl_vector_free (v);
1184 gsl_matrix_free (w);
1185 (void) gsl_set_error_handler (save_handler);
1186 }
1187
1188 //! @brief OP * = ([] REAL, [, ] REAL) [] REAL
1189
1190 void genie_vector_times_matrix (NODE_T * p)
1191 {
1192 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1193 gsl_matrix *w = pop_matrix (p, A68_TRUE);
1194 ASSERT_GSL (gsl_matrix_transpose (w));
1195 gsl_vector *u = pop_vector (p, A68_TRUE);
1196 int len = (int) (SIZE (u));
1197 gsl_vector *v = gsl_vector_calloc ((size_t) len);
1198 gsl_vector_set_zero (v);
1199 ASSERT_GSL (gsl_blas_dgemv (CblasNoTrans, 1.0, w, u, 0.0, v));
1200 push_vector (p, v);
1201 gsl_vector_free (u);
1202 gsl_vector_free (v);
1203 gsl_matrix_free (w);
1204 (void) gsl_set_error_handler (save_handler);
1205 }
1206
1207 //! @brief OP * = ([, ] REAL, [, ] REAL) [, ] REAL
1208
1209 void genie_matrix_times_matrix (NODE_T * p)
1210 {
1211 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1212 gsl_matrix *v = pop_matrix (p, A68_TRUE);
1213 gsl_matrix *u = pop_matrix (p, A68_TRUE);
1214 int len2 = (int) (SIZE2 (v)), len1 = (int) (SIZE1 (u));
1215 gsl_matrix *w = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
1216 ASSERT_GSL (gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, u, v, 0.0, w));
1217 push_matrix (p, w);
1218 gsl_matrix_free (u);
1219 gsl_matrix_free (v);
1220 gsl_matrix_free (w);
1221 (void) gsl_set_error_handler (save_handler);
1222 }
1223
1224 //! @brief OP * = ([, ] COMPLEX, [] COMPLEX) [] COMPLEX
1225
1226 void genie_matrix_complex_times_vector (NODE_T * p)
1227 {
1228 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1229 gsl_complex zero, one;
1230 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1231 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1232 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
1233 gsl_matrix_complex *w = pop_matrix_complex (p, A68_TRUE);
1234 int len = (int) (SIZE (u));
1235 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
1236 ASSERT_GSL (gsl_blas_zgemv (CblasNoTrans, one, w, u, zero, v));
1237 push_vector_complex (p, v);
1238 gsl_vector_complex_free (u);
1239 gsl_vector_complex_free (v);
1240 gsl_matrix_complex_free (w);
1241 (void) gsl_set_error_handler (save_handler);
1242 }
1243
1244 //! @brief OP * = ([] COMPLEX, [, ] COMPLEX) [] COMPLEX
1245
1246 void genie_vector_complex_times_matrix (NODE_T * p)
1247 {
1248 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1249 gsl_complex zero, one;
1250 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1251 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1252 gsl_matrix_complex *w = pop_matrix_complex (p, A68_TRUE);
1253 ASSERT_GSL (gsl_matrix_complex_transpose (w));
1254 gsl_vector_complex *u = pop_vector_complex (p, A68_TRUE);
1255 int len = (int) (SIZE (u));
1256 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
1257 ASSERT_GSL (gsl_blas_zgemv (CblasNoTrans, one, w, u, zero, v));
1258 push_vector_complex (p, v);
1259 gsl_vector_complex_free (u);
1260 gsl_vector_complex_free (v);
1261 gsl_matrix_complex_free (w);
1262 (void) gsl_set_error_handler (save_handler);
1263 }
1264
1265 //! @brief OP * = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
1266
1267 void genie_matrix_complex_times_matrix (NODE_T * p)
1268 {
1269 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1270 gsl_complex zero, one;
1271 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1272 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1273 gsl_matrix_complex *v = pop_matrix_complex (p, A68_TRUE);
1274 gsl_matrix_complex *u = pop_matrix_complex (p, A68_TRUE);
1275 int len1 = (int) (SIZE2 (v)), len2 = (int) (SIZE1 (u));
1276 gsl_matrix_complex *w = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
1277 ASSERT_GSL (gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, one, u, v, zero, w));
1278 gsl_matrix_complex_free (u);
1279 gsl_matrix_complex_free (v);
1280 push_matrix_complex (p, w);
1281 gsl_matrix_complex_free (w);
1282 (void) gsl_set_error_handler (save_handler);
1283 }
1284
1285 #endif
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
