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 (a68g_bufprt (A68G (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
44 } else {
45 ASSERT (a68g_bufprt (A68G (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
46 }
47 diagnostic (A68G_RUNTIME_ERROR, A68G (f_entry), ERROR_TORRIX, A68G (edit_line), gsl_strerror (gsl_errno));
48 exit_genie (A68G (f_entry), A68G_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 A68G_REF desc; A68G_ARRAY *arr; A68G_TUPLE *tup;
56 POP_REF (p, &desc);
57 CHECK_REF (p, desc, M_ROW_INT);
58 GET_DESCRIPTOR (arr, tup, &desc);
59 size_t 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 A68G_INT *x = (A68G_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 size_t len = SIZE (v);
79 A68G_REF desc, row; A68G_ARRAY arr; A68G_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 (size_t k = 0; k < len; k++, idx += inc) {
85 A68G_INT *x = (A68G_INT *) (base + idx);
86 STATUS (x) = INIT_MASK;
87 VALUE (x) = (INT_T) gsl_permutation_get (v, 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 A68G_REF desc; A68G_ARRAY *arr; A68G_TUPLE *tup;
97 POP_REF (p, &desc);
98 CHECK_REF (p, desc, M_ROW_REAL);
99 GET_DESCRIPTOR (arr, tup, &desc);
100 size_t 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 A68G_REAL *x = (A68G_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 A68G_REF desc; A68G_ARRAY *arr; A68G_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 size_t 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 A68G_REAL *x = (A68G_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 A68G_REF desc; A68G_ARRAY *arr; A68G_TUPLE *tup;
160 POP_REF (p, &desc);
161 CHECK_REF (p, desc, M_ROW_COMPLEX);
162 GET_DESCRIPTOR (arr, tup, &desc);
163 size_t 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 A68G_REAL *re = (A68G_REAL *) (base + idx);
171 A68G_REAL *im = (A68G_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 size_t len = SIZE (v);
187 A68G_REF desc, row; A68G_ARRAY arr; A68G_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 A68G_REAL *re = (A68G_REAL *) (base + idx);
194 A68G_REAL *im = (A68G_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 A68G_REF desc; A68G_ARRAY *arr; A68G_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 size_t 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 A68G_REAL *re = (A68G_REAL *) (base + idx2);
223 A68G_REAL *im = (A68G_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 size_t len1 = SIZE1 (a), len2 = SIZE2 (a);
240 A68G_REF desc, row; A68G_ARRAY arr; A68G_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 A68G_REAL *re = (A68G_REAL *) (base + idx2);
259 A68G_REAL *im = (A68G_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 A68G_REF dst, src, *save = (A68G_REF *) STACK_OFFSET (-parm_size);
277 dst = *save;
278 CHECK_REF (p, dst, m);
279 *save = *DEREF (A68G_ROW, &dst);
280 STATUS (&src) = (STATUS_MASK_T) (INIT_MASK | IN_STACK_MASK);
281 OFFSET (&src) = A68G_SP - parm_size;
282 (*op) (p);
283 if (IS_REF (m)) {
284 genie_store (p, SUB (m), &dst, &src);
285 } else {
286 ABEND (A68G_TRUE, ERROR_INTERNAL_CONSISTENCY, NO_TEXT);
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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_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, A68G_TRUE);
463 gsl_permutation *q = gsl_permutation_calloc (SIZE1 (u));
464 int sign, rc = gsl_linalg_LU_decomp (u, q, &sign);
465 if (rc != GSL_SUCCESS) {
466 PUSH_VALUE (p, 0, A68G_REAL);
467 } else {
468 REAL_T det = gsl_linalg_LU_det (u, sign);
469 if (a68g_isnan_real (det)) {
470 PUSH_VALUE (p, 0, A68G_REAL);
471 } else {
472 PUSH_VALUE (p, det, A68G_REAL);
473 }
474 }
475 gsl_matrix_free (u);
476 gsl_permutation_free (q);
477 (void) gsl_set_error_handler (save_handler);
478 }
479
480 //! @brief OP DET = ([, ] COMPLEX) COMPLEX
481
482 void genie_matrix_complex_det (NODE_T * p)
483 {
484 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
485 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
486 gsl_permutation *q = gsl_permutation_calloc (SIZE1 (u));
487 int sign, rc = gsl_linalg_complex_LU_decomp (u, q, &sign);
488 if (rc != GSL_SUCCESS) {
489 PUSH_COMPLEX_VALUE (p, 0, 0);
490 } else {
491 gsl_complex det = gsl_linalg_complex_LU_det (u, sign);
492 REAL_T re_det = GSL_REAL (det), im_det = GSL_IMAG (det);
493 if (a68g_isnan_real (re_det) || a68g_isnan_real (im_det)) {
494 PUSH_COMPLEX_VALUE (p, 0, 0);
495 } else {
496 PUSH_COMPLEX_VALUE (p, re_det, im_det);
497 }
498 }
499 gsl_matrix_complex_free (u);
500 gsl_permutation_free (q);
501 (void) gsl_set_error_handler (save_handler);
502 }
503
504 //! @brief OP TRACE = ([, ] REAL) REAL
505
506 void genie_matrix_trace (NODE_T * p)
507 {
508 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
509 gsl_matrix *a = pop_matrix (p, A68G_TRUE);
510 size_t len1 = SIZE1 (a), len2 = SIZE2 (a);
511 if (len1 != len2) {
512 torrix_error_handler ("cannot calculate trace", __FILE__, __LINE__, GSL_ENOTSQR);
513 }
514 REAL_T sum = 0.0;
515 for (int k = 0; k < len1; k++) {
516 sum += gsl_matrix_get (a, (size_t) k, (size_t) k);
517 }
518 PUSH_VALUE (p, sum, A68G_REAL);
519 gsl_matrix_free (a);
520 (void) gsl_set_error_handler (save_handler);
521 }
522
523 //! @brief OP TRACE = ([, ] COMPLEX) COMPLEX
524
525 void genie_matrix_complex_trace (NODE_T * p)
526 {
527 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
528 gsl_matrix_complex *a = pop_matrix_complex (p, A68G_TRUE);
529 size_t len1 = SIZE1 (a), len2 = SIZE2 (a);
530 if (len1 != len2) {
531 torrix_error_handler ("cannot calculate trace", __FILE__, __LINE__, GSL_ENOTSQR);
532 }
533 gsl_complex sum;
534 GSL_SET_COMPLEX (&sum, 0.0, 0.0);
535 for (int k = 0; k < len1; k++) {
536 sum = gsl_complex_add (sum, gsl_matrix_complex_get (a, (size_t) k, (size_t) k));
537 }
538 PUSH_COMPLEX_VALUE (p, GSL_REAL (sum), GSL_IMAG (sum));
539 gsl_matrix_complex_free (a);
540 (void) gsl_set_error_handler (save_handler);
541 }
542
543 //! @brief OP - = ([] COMPLEX) [] COMPLEX
544
545 void genie_vector_complex_minus (NODE_T * p)
546 {
547 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
548 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
549 gsl_blas_zdscal (-1, u);
550 push_vector_complex (p, u);
551 gsl_vector_complex_free (u);
552 (void) gsl_set_error_handler (save_handler);
553 }
554
555 //! @brief OP - = ([, ] COMPLEX) [, ] COMPLEX
556
557 void genie_matrix_complex_minus (NODE_T * p)
558 {
559 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
560 gsl_complex one;
561 GSL_SET_COMPLEX (&one, -1.0, 0.0);
562 gsl_matrix_complex *a = pop_matrix_complex (p, A68G_TRUE);
563 ASSERT_GSL (gsl_matrix_complex_scale (a, one));
564 push_matrix_complex (p, a);
565 gsl_matrix_complex_free (a);
566 (void) gsl_set_error_handler (save_handler);
567 }
568
569 //! @brief OP + = ([] REAL, [] REAL) [] REAL
570
571 void genie_vector_add (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, A68G_TRUE);
575 gsl_vector *u = pop_vector (p, A68G_TRUE);
576 ASSERT_GSL (gsl_vector_add (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) [] REAL
584
585 void genie_vector_sub (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, A68G_TRUE);
589 gsl_vector *u = pop_vector (p, A68G_TRUE);
590 ASSERT_GSL (gsl_vector_sub (u, v));
591 push_vector (p, u);
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_eq (NODE_T * p)
600 {
601 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
602 gsl_vector *v = pop_vector (p, A68G_TRUE);
603 gsl_vector *u = pop_vector (p, A68G_TRUE);
604 ASSERT_GSL (gsl_vector_sub (u, v));
605 PUSH_VALUE (p, (BOOL_T) (gsl_vector_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
606 gsl_vector_free (u);
607 gsl_vector_free (v);
608 (void) gsl_set_error_handler (save_handler);
609 }
610
611 //! @brief OP /= = ([] REAL, [] REAL) BOOL
612
613 void genie_vector_ne (NODE_T * p)
614 {
615 genie_vector_eq (p);
616 genie_not_bool (p);
617 }
618
619 //! @brief OP +:= = (REF [] REAL, [] REAL) REF [] REAL
620
621 void genie_vector_plusab (NODE_T * p)
622 {
623 op_ab_torrix (p, M_REF_ROW_REAL, M_ROW_REAL, genie_vector_add);
624 }
625
626 //! @brief OP -:= = (REF [] REAL, [] REAL) REF [] REAL
627
628 void genie_vector_minusab (NODE_T * p)
629 {
630 op_ab_torrix (p, M_REF_ROW_REAL, M_ROW_REAL, genie_vector_sub);
631 }
632
633 //! @brief OP + = ([, ] REAL, [, ] REAL) [, ] REAL
634
635 void genie_matrix_add (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, A68G_TRUE);
639 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
640 ASSERT_GSL (gsl_matrix_add (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) [, ] REAL
648
649 void genie_matrix_sub (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, A68G_TRUE);
653 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
654 ASSERT_GSL (gsl_matrix_sub (u, v));
655 push_matrix (p, u);
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_eq (NODE_T * p)
664 {
665 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
666 gsl_matrix *v = pop_matrix (p, A68G_TRUE);
667 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
668 ASSERT_GSL (gsl_matrix_sub (u, v));
669 PUSH_VALUE (p, (BOOL_T) (gsl_matrix_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
670 gsl_matrix_free (u);
671 gsl_matrix_free (v);
672 (void) gsl_set_error_handler (save_handler);
673 }
674
675 //! @brief OP /= = ([, ] REAL, [, ] REAL) [, ] BOOL
676
677 void genie_matrix_ne (NODE_T * p)
678 {
679 genie_matrix_eq (p);
680 genie_not_bool (p);
681 }
682
683 //! @brief OP +:= = (REF [, ] REAL, [, ] REAL) [, ] REAL
684
685 void genie_matrix_plusab (NODE_T * p)
686 {
687 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_ROW_ROW_REAL, genie_matrix_add);
688 }
689
690 //! @brief OP -:= = (REF [, ] REAL, [, ] REAL) [, ] REAL
691
692 void genie_matrix_minusab (NODE_T * p)
693 {
694 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_ROW_ROW_REAL, genie_matrix_sub);
695 }
696
697 //! @brief OP + = ([] COMPLEX, [] COMPLEX) [] COMPLEX
698
699 void genie_vector_complex_add (NODE_T * p)
700 {
701 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
702 gsl_complex one;
703 GSL_SET_COMPLEX (&one, 1.0, 0.0);
704 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
705 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
706 ASSERT_GSL (gsl_blas_zaxpy (one, u, v));
707 push_vector_complex (p, v);
708 gsl_vector_complex_free (u);
709 gsl_vector_complex_free (v);
710 (void) gsl_set_error_handler (save_handler);
711 }
712
713 //! @brief OP - = ([] COMPLEX, [] COMPLEX) [] COMPLEX
714
715 void genie_vector_complex_sub (NODE_T * p)
716 {
717 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
718 gsl_complex one;
719 GSL_SET_COMPLEX (&one, -1.0, 0.0);
720 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
721 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
722 ASSERT_GSL (gsl_blas_zaxpy (one, v, u));
723 push_vector_complex (p, u);
724 gsl_vector_complex_free (u);
725 gsl_vector_complex_free (v);
726 (void) gsl_set_error_handler (save_handler);
727 }
728
729 //! @brief OP = = ([] COMPLEX, [] COMPLEX) BOOL
730
731 void genie_vector_complex_eq (NODE_T * p)
732 {
733 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
734 gsl_complex one;
735 GSL_SET_COMPLEX (&one, -1.0, 0.0);
736 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
737 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
738 ASSERT_GSL (gsl_blas_zaxpy (one, v, u));
739 PUSH_VALUE (p, (BOOL_T) (gsl_vector_complex_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
740 gsl_vector_complex_free (u);
741 gsl_vector_complex_free (v);
742 (void) gsl_set_error_handler (save_handler);
743 }
744
745 //! @brief OP /= = ([] COMPLEX, [] COMPLEX) BOOL
746
747 void genie_vector_complex_ne (NODE_T * p)
748 {
749 genie_vector_complex_eq (p);
750 genie_not_bool (p);
751 }
752
753 //! @brief OP +:= = (REF [] COMPLEX, [] COMPLEX) [] COMPLEX
754
755 void genie_vector_complex_plusab (NODE_T * p)
756 {
757 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_ROW_COMPLEX, genie_vector_complex_add);
758 }
759
760 //! @brief OP -:= = (REF [] COMPLEX, [] COMPLEX) [] COMPLEX
761
762 void genie_vector_complex_minusab (NODE_T * p)
763 {
764 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_ROW_COMPLEX, genie_vector_complex_sub);
765 }
766
767 //! @brief OP + = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
768
769 void genie_matrix_complex_add (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, A68G_TRUE);
773 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
774 ASSERT_GSL (gsl_matrix_complex_add (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) [, ] COMPLEX
782
783 void genie_matrix_complex_sub (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, A68G_TRUE);
787 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
788 ASSERT_GSL (gsl_matrix_complex_sub (u, v));
789 push_matrix_complex (p, u);
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_eq (NODE_T * p)
798 {
799 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
800 gsl_matrix_complex *v = pop_matrix_complex (p, A68G_TRUE);
801 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
802 ASSERT_GSL (gsl_matrix_complex_sub (u, v));
803 PUSH_VALUE (p, (BOOL_T) (gsl_matrix_complex_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
804 gsl_matrix_complex_free (u);
805 gsl_matrix_complex_free (v);
806 (void) gsl_set_error_handler (save_handler);
807 }
808
809 //! @brief OP /= = ([, ] COMPLEX, [, ] COMPLEX) BOOL
810
811 void genie_matrix_complex_ne (NODE_T * p)
812 {
813 genie_matrix_complex_eq (p);
814 genie_not_bool (p);
815 }
816
817 //! @brief OP +:= = (REF [, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
818
819 void genie_matrix_complex_plusab (NODE_T * p)
820 {
821 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_ROW_ROW_COMPLEX, genie_matrix_complex_add);
822 }
823
824 //! @brief OP -:= = (REF [, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
825
826 void genie_matrix_complex_minusab (NODE_T * p)
827 {
828 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_ROW_ROW_COMPLEX, genie_matrix_complex_sub);
829 }
830
831 //! @brief OP * = ([] REAL, REAL) [] REAL
832
833 void genie_vector_scale_real (NODE_T * p)
834 {
835 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
836 A68G_REAL v;
837 POP_OBJECT (p, &v, A68G_REAL);
838 gsl_vector *u = pop_vector (p, A68G_TRUE);
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_real_scale_vector (NODE_T * p)
848 {
849 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
850 gsl_vector *u = pop_vector (p, A68G_TRUE);
851 A68G_REAL v;
852 POP_OBJECT (p, &v, A68G_REAL);
853 ASSERT_GSL (gsl_vector_scale (u, VALUE (&v)));
854 push_vector (p, u);
855 gsl_vector_free (u);
856 (void) gsl_set_error_handler (save_handler);
857 }
858
859 //! @brief OP * = ([, ] REAL, REAL) [, ] REAL
860
861 void genie_matrix_scale_real (NODE_T * p)
862 {
863 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
864 A68G_REAL v;
865 POP_OBJECT (p, &v, A68G_REAL);
866 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
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 * = (REAL, [, ] REAL) [, ] REAL
874
875 void genie_real_scale_matrix (NODE_T * p)
876 {
877 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
878 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
879 A68G_REAL v;
880 POP_OBJECT (p, &v, A68G_REAL);
881 ASSERT_GSL (gsl_matrix_scale (u, VALUE (&v)));
882 push_matrix (p, u);
883 gsl_matrix_free (u);
884 (void) gsl_set_error_handler (save_handler);
885 }
886
887 //! @brief OP * = ([] COMPLEX, COMPLEX) [] COMPLEX
888
889 void genie_vector_complex_scale_complex (NODE_T * p)
890 {
891 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
892 A68G_REAL re, im; gsl_complex v;
893 POP_OBJECT (p, &im, A68G_REAL);
894 POP_OBJECT (p, &re, A68G_REAL);
895 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
896 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
897 gsl_blas_zscal (v, u);
898 push_vector_complex (p, u);
899 gsl_vector_complex_free (u);
900 (void) gsl_set_error_handler (save_handler);
901 }
902
903 //! @brief OP * = (COMPLEX, [] COMPLEX) [] COMPLEX
904
905 void genie_complex_scale_vector_complex (NODE_T * p)
906 {
907 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
908 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
909 A68G_REAL re, im; gsl_complex v;
910 POP_OBJECT (p, &im, A68G_REAL);
911 POP_OBJECT (p, &re, A68G_REAL);
912 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
913 gsl_blas_zscal (v, u);
914 push_vector_complex (p, u);
915 gsl_vector_complex_free (u);
916 (void) gsl_set_error_handler (save_handler);
917 }
918
919 //! @brief OP * = ([, ] COMPLEX, COMPLEX) [, ] COMPLEX
920
921 void genie_matrix_complex_scale_complex (NODE_T * p)
922 {
923 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
924 A68G_REAL re, im; gsl_complex v;
925 POP_OBJECT (p, &im, A68G_REAL);
926 POP_OBJECT (p, &re, A68G_REAL);
927 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
928 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
929 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
930 push_matrix_complex (p, u);
931 gsl_matrix_complex_free (u);
932 (void) gsl_set_error_handler (save_handler);
933 }
934
935 //! @brief OP * = (COMPLEX, [, ] COMPLEX) [, ] COMPLEX
936
937 void genie_complex_scale_matrix_complex (NODE_T * p)
938 {
939 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
940 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
941 A68G_REAL re, im; gsl_complex v;
942 POP_OBJECT (p, &im, A68G_REAL);
943 POP_OBJECT (p, &re, A68G_REAL);
944 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
945 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
946 push_matrix_complex (p, u);
947 gsl_matrix_complex_free (u);
948 (void) gsl_set_error_handler (save_handler);
949 }
950
951 //! @brief OP *:= (REF [] REAL, REAL) REF [] REAL
952
953 void genie_vector_scale_real_ab (NODE_T * p)
954 {
955 op_ab_torrix (p, M_REF_ROW_REAL, M_REAL, genie_vector_scale_real);
956 }
957
958 //! @brief OP *:= (REF [, ] REAL, REAL) REF [, ] REAL
959
960 void genie_matrix_scale_real_ab (NODE_T * p)
961 {
962 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_REAL, genie_matrix_scale_real);
963 }
964
965 //! @brief OP *:= (REF [] COMPLEX, COMPLEX) REF [] COMPLEX
966
967 void genie_vector_complex_scale_complex_ab (NODE_T * p)
968 {
969 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_COMPLEX, genie_vector_complex_scale_complex);
970 }
971
972 //! @brief OP *:= (REF [, ] COMPLEX, COMPLEX) REF [, ] COMPLEX
973
974 void genie_matrix_complex_scale_complex_ab (NODE_T * p)
975 {
976 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_COMPLEX, genie_matrix_complex_scale_complex);
977 }
978
979 //! @brief OP / = ([] REAL, REAL) [] REAL
980
981 void genie_vector_div_real (NODE_T * p)
982 {
983 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
984 A68G_REAL v;
985 POP_OBJECT (p, &v, A68G_REAL);
986 if (VALUE (&v) == 0.0) {
987 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_REAL);
988 exit_genie (p, A68G_RUNTIME_ERROR);
989 }
990 gsl_vector *u = pop_vector (p, A68G_TRUE);
991 ASSERT_GSL (gsl_vector_scale (u, 1.0 / VALUE (&v)));
992 push_vector (p, u);
993 gsl_vector_free (u);
994 (void) gsl_set_error_handler (save_handler);
995 }
996
997 //! @brief OP / = ([, ] REAL, REAL) [, ] REAL
998
999 void genie_matrix_div_real (NODE_T * p)
1000 {
1001 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1002 A68G_REAL v;
1003 POP_OBJECT (p, &v, A68G_REAL);
1004 if (VALUE (&v) == 0.0) {
1005 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_ROW_REAL);
1006 exit_genie (p, A68G_RUNTIME_ERROR);
1007 }
1008 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
1009 ASSERT_GSL (gsl_matrix_scale (u, 1.0 / VALUE (&v)));
1010 push_matrix (p, u);
1011 gsl_matrix_free (u);
1012 (void) gsl_set_error_handler (save_handler);
1013 }
1014
1015 //! @brief OP / = ([] COMPLEX, COMPLEX) [] COMPLEX
1016
1017 void genie_vector_complex_div_complex (NODE_T * p)
1018 {
1019 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1020 A68G_REAL re, im; gsl_complex v;
1021 POP_OBJECT (p, &im, A68G_REAL);
1022 POP_OBJECT (p, &re, A68G_REAL);
1023 if (VALUE (&re) == 0.0 && VALUE (&im) == 0.0) {
1024 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_COMPLEX);
1025 exit_genie (p, A68G_RUNTIME_ERROR);
1026 }
1027 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
1028 v = gsl_complex_inverse (v);
1029 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1030 gsl_blas_zscal (v, u);
1031 push_vector_complex (p, u);
1032 gsl_vector_complex_free (u);
1033 (void) gsl_set_error_handler (save_handler);
1034 }
1035
1036 //! @brief OP / = ([, ] COMPLEX, COMPLEX) [, ] COMPLEX
1037
1038 void genie_matrix_complex_div_complex (NODE_T * p)
1039 {
1040 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1041 A68G_REAL re, im; gsl_complex v;
1042 POP_OBJECT (p, &im, A68G_REAL);
1043 POP_OBJECT (p, &re, A68G_REAL);
1044 if (VALUE (&re) == 0.0 && VALUE (&im) == 0.0) {
1045 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_ROW_COMPLEX);
1046 exit_genie (p, A68G_RUNTIME_ERROR);
1047 }
1048 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
1049 v = gsl_complex_inverse (v);
1050 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
1051 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
1052 push_matrix_complex (p, u);
1053 gsl_matrix_complex_free (u);
1054 (void) gsl_set_error_handler (save_handler);
1055 }
1056
1057 //! @brief OP /:= (REF [] REAL, REAL) REF [] REAL
1058
1059 void genie_vector_div_real_ab (NODE_T * p)
1060 {
1061 op_ab_torrix (p, M_REF_ROW_REAL, M_REAL, genie_vector_div_real);
1062 }
1063
1064 //! @brief OP /:= (REF [, ] REAL, REAL) REF [, ] REAL
1065
1066 void genie_matrix_div_real_ab (NODE_T * p)
1067 {
1068 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_REAL, genie_matrix_div_real);
1069 }
1070
1071 //! @brief OP /:= (REF [] COMPLEX, COMPLEX) REF [] COMPLEX
1072
1073 void genie_vector_complex_div_complex_ab (NODE_T * p)
1074 {
1075 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_COMPLEX, genie_vector_complex_div_complex);
1076 }
1077
1078 //! @brief OP /:= (REF [, ] COMPLEX, COMPLEX) REF [, ] COMPLEX
1079
1080 void genie_matrix_complex_div_complex_ab (NODE_T * p)
1081 {
1082 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_COMPLEX, genie_matrix_complex_div_complex);
1083 }
1084
1085 //! @brief OP * = ([] REAL, [] REAL) REAL
1086
1087 void genie_vector_dot (NODE_T * p)
1088 {
1089 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1090 gsl_vector *v = pop_vector (p, A68G_TRUE);
1091 gsl_vector *u = pop_vector (p, A68G_TRUE);
1092 REAL_T w;
1093 ASSERT_GSL (gsl_blas_ddot (u, v, &w));
1094 PUSH_VALUE (p, w, A68G_REAL);
1095 gsl_vector_free (u);
1096 gsl_vector_free (v);
1097 (void) gsl_set_error_handler (save_handler);
1098 }
1099
1100 //! @brief OP * = ([] COMPLEX, [] COMPLEX) COMPLEX
1101
1102 void genie_vector_complex_dot (NODE_T * p)
1103 {
1104 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1105 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
1106 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1107 gsl_complex w;
1108 ASSERT_GSL (gsl_blas_zdotc (u, v, &w));
1109 PUSH_COMPLEX_VALUE (p, GSL_REAL (w), GSL_IMAG (w));
1110 gsl_vector_complex_free (u);
1111 gsl_vector_complex_free (v);
1112 (void) gsl_set_error_handler (save_handler);
1113 }
1114
1115 //! @brief OP NORM = ([] REAL) REAL
1116
1117 void genie_vector_norm (NODE_T * p)
1118 {
1119 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1120 gsl_vector *u = pop_vector (p, A68G_TRUE);
1121 PUSH_VALUE (p, gsl_blas_dnrm2 (u), A68G_REAL);
1122 gsl_vector_free (u);
1123 (void) gsl_set_error_handler (save_handler);
1124 }
1125
1126 //! @brief OP NORM = ([] COMPLEX) COMPLEX
1127
1128 void genie_vector_complex_norm (NODE_T * p)
1129 {
1130 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1131 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1132 PUSH_VALUE (p, gsl_blas_dznrm2 (u), A68G_REAL);
1133 gsl_vector_complex_free (u);
1134 (void) gsl_set_error_handler (save_handler);
1135 }
1136
1137 //! @brief OP DYAD = ([] REAL, [] REAL) [, ] REAL
1138
1139 void genie_vector_dyad (NODE_T * p)
1140 {
1141 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1142 gsl_vector *v = pop_vector (p, A68G_TRUE);
1143 gsl_vector *u = pop_vector (p, A68G_TRUE);
1144 size_t len1 = SIZE (u), len2 = SIZE (v);
1145 gsl_matrix *w = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
1146 for (int j = 0; j < len1; j++) {
1147 REAL_T uj = gsl_vector_get (u, (size_t) j);
1148 for (int k = 0; k < len2; k++) {
1149 REAL_T vk = gsl_vector_get (v, (size_t) k);
1150 gsl_matrix_set (w, (size_t) j, (size_t) k, uj * vk);
1151 }
1152 }
1153 push_matrix (p, w);
1154 gsl_vector_free (u);
1155 gsl_vector_free (v);
1156 gsl_matrix_free (w);
1157 (void) gsl_set_error_handler (save_handler);
1158 }
1159
1160 //! @brief OP DYAD = ([] COMPLEX, [] COMPLEX) [, ] COMPLEX
1161
1162 void genie_vector_complex_dyad (NODE_T * p)
1163 {
1164 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1165 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
1166 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1167 size_t len1 = SIZE (u), len2 = SIZE (v);
1168 gsl_matrix_complex *w = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
1169 for (int j = 0; j < len1; j++) {
1170 gsl_complex uj = gsl_vector_complex_get (u, (size_t) j);
1171 for (int k = 0; k < len2; k++) {
1172 gsl_complex vk = gsl_vector_complex_get (u, (size_t) k);
1173 gsl_matrix_complex_set (w, (size_t) j, (size_t) k, gsl_complex_mul (uj, vk));
1174 }
1175 }
1176 push_matrix_complex (p, w);
1177 gsl_vector_complex_free (u);
1178 gsl_vector_complex_free (v);
1179 gsl_matrix_complex_free (w);
1180 (void) gsl_set_error_handler (save_handler);
1181 }
1182
1183 //! @brief OP * = ([, ] REAL, [] REAL) [] REAL
1184
1185 void genie_matrix_times_vector (NODE_T * p)
1186 {
1187 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1188 gsl_vector *u = pop_vector (p, A68G_TRUE);
1189 gsl_matrix *w = pop_matrix (p, A68G_TRUE);
1190 size_t len = SIZE (u);
1191 gsl_vector *v = gsl_vector_calloc ((size_t) len);
1192 gsl_vector_set_zero (v);
1193 ASSERT_GSL (gsl_blas_dgemv (CblasNoTrans, 1.0, w, u, 0.0, v));
1194 push_vector (p, v);
1195 gsl_vector_free (u);
1196 gsl_vector_free (v);
1197 gsl_matrix_free (w);
1198 (void) gsl_set_error_handler (save_handler);
1199 }
1200
1201 //! @brief OP * = ([] REAL, [, ] REAL) [] REAL
1202
1203 void genie_vector_times_matrix (NODE_T * p)
1204 {
1205 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1206 gsl_matrix *w = pop_matrix (p, A68G_TRUE);
1207 ASSERT_GSL (gsl_matrix_transpose (w));
1208 gsl_vector *u = pop_vector (p, A68G_TRUE);
1209 size_t len = SIZE (u);
1210 gsl_vector *v = gsl_vector_calloc ((size_t) len);
1211 gsl_vector_set_zero (v);
1212 ASSERT_GSL (gsl_blas_dgemv (CblasNoTrans, 1.0, w, u, 0.0, v));
1213 push_vector (p, v);
1214 gsl_vector_free (u);
1215 gsl_vector_free (v);
1216 gsl_matrix_free (w);
1217 (void) gsl_set_error_handler (save_handler);
1218 }
1219
1220 //! @brief OP * = ([, ] REAL, [, ] REAL) [, ] REAL
1221
1222 void genie_matrix_times_matrix (NODE_T * p)
1223 {
1224 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1225 gsl_matrix *v = pop_matrix (p, A68G_TRUE);
1226 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
1227 size_t len2 = SIZE2 (v), len1 = SIZE1 (u);
1228 gsl_matrix *w = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
1229 ASSERT_GSL (gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, u, v, 0.0, w));
1230 push_matrix (p, w);
1231 gsl_matrix_free (u);
1232 gsl_matrix_free (v);
1233 gsl_matrix_free (w);
1234 (void) gsl_set_error_handler (save_handler);
1235 }
1236
1237 //! @brief OP * = ([, ] COMPLEX, [] COMPLEX) [] COMPLEX
1238
1239 void genie_matrix_complex_times_vector (NODE_T * p)
1240 {
1241 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1242 gsl_complex zero, one;
1243 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1244 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1245 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1246 gsl_matrix_complex *w = pop_matrix_complex (p, A68G_TRUE);
1247 size_t len = SIZE (u);
1248 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
1249 ASSERT_GSL (gsl_blas_zgemv (CblasNoTrans, one, w, u, zero, v));
1250 push_vector_complex (p, v);
1251 gsl_vector_complex_free (u);
1252 gsl_vector_complex_free (v);
1253 gsl_matrix_complex_free (w);
1254 (void) gsl_set_error_handler (save_handler);
1255 }
1256
1257 //! @brief OP * = ([] COMPLEX, [, ] COMPLEX) [] COMPLEX
1258
1259 void genie_vector_complex_times_matrix (NODE_T * p)
1260 {
1261 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1262 gsl_complex zero, one;
1263 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1264 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1265 gsl_matrix_complex *w = pop_matrix_complex (p, A68G_TRUE);
1266 ASSERT_GSL (gsl_matrix_complex_transpose (w));
1267 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1268 size_t len = SIZE (u);
1269 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
1270 ASSERT_GSL (gsl_blas_zgemv (CblasNoTrans, one, w, u, zero, v));
1271 push_vector_complex (p, v);
1272 gsl_vector_complex_free (u);
1273 gsl_vector_complex_free (v);
1274 gsl_matrix_complex_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_matrix (NODE_T * p)
1281 {
1282 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1283 gsl_complex zero, one;
1284 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1285 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1286 gsl_matrix_complex *v = pop_matrix_complex (p, A68G_TRUE);
1287 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
1288 size_t len1 = SIZE2 (v), len2 = SIZE1 (u);
1289 gsl_matrix_complex *w = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
1290 ASSERT_GSL (gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, one, u, v, zero, w));
1291 gsl_matrix_complex_free (u);
1292 gsl_matrix_complex_free (v);
1293 push_matrix_complex (p, w);
1294 gsl_matrix_complex_free (w);
1295 (void) gsl_set_error_handler (save_handler);
1296 }
1297
1298 #endif
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
