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