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, __func__);
287 }
288 *save = dst;
289 }
290
291 //! @brief PROC vector echo = ([] REAL) [] REAL
292
293 void genie_vector_echo (NODE_T * p)
294 {
295 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
296 gsl_vector *u = pop_vector (p, 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;
465 ASSERT_GSL (gsl_linalg_LU_decomp (u, q, &sign));
466 REAL_T det = gsl_linalg_LU_det (u, sign);
467 if (a68g_isnan_real (det)) { // Singularity.
468 PUSH_VALUE (p, 0, A68G_REAL);
469 } else {
470 PUSH_VALUE (p, gsl_linalg_LU_det (u, sign), A68G_REAL);
471 }
472 gsl_matrix_free (u);
473 gsl_permutation_free (q);
474 (void) gsl_set_error_handler (save_handler);
475 }
476
477 //! @brief OP DET = ([, ] COMPLEX) COMPLEX
478
479 void genie_matrix_complex_det (NODE_T * p)
480 {
481 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
482 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
483 gsl_permutation *q = gsl_permutation_calloc (SIZE1 (u));
484 int sign;
485 ASSERT_GSL (gsl_linalg_complex_LU_decomp (u, q, &sign));
486 gsl_complex det = gsl_linalg_complex_LU_det (u, sign);
487 REAL_T re_det = GSL_REAL (det);
488 if (a68g_isnan_real (re_det)) { // Singularity.
489 PUSH_VALUE (p, 0, A68G_REAL);
490 } else {
491 PUSH_VALUE (p, re_det, A68G_REAL);
492 }
493 REAL_T im_det = GSL_IMAG (det);
494 if (a68g_isnan_real (im_det)) { // Singularity.
495 PUSH_VALUE (p, 0, A68G_REAL);
496 } else {
497 PUSH_VALUE (p, im_det, A68G_REAL);
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_VALUE (p, GSL_REAL (sum), A68G_REAL);
539 PUSH_VALUE (p, GSL_IMAG (sum), A68G_REAL);
540 gsl_matrix_complex_free (a);
541 (void) gsl_set_error_handler (save_handler);
542 }
543
544 //! @brief OP - = ([] COMPLEX) [] COMPLEX
545
546 void genie_vector_complex_minus (NODE_T * p)
547 {
548 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
549 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
550 gsl_blas_zdscal (-1, u);
551 push_vector_complex (p, u);
552 gsl_vector_complex_free (u);
553 (void) gsl_set_error_handler (save_handler);
554 }
555
556 //! @brief OP - = ([, ] COMPLEX) [, ] COMPLEX
557
558 void genie_matrix_complex_minus (NODE_T * p)
559 {
560 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
561 gsl_complex one;
562 GSL_SET_COMPLEX (&one, -1.0, 0.0);
563 gsl_matrix_complex *a = pop_matrix_complex (p, A68G_TRUE);
564 ASSERT_GSL (gsl_matrix_complex_scale (a, one));
565 push_matrix_complex (p, a);
566 gsl_matrix_complex_free (a);
567 (void) gsl_set_error_handler (save_handler);
568 }
569
570 //! @brief OP + = ([] REAL, [] REAL) [] REAL
571
572 void genie_vector_add (NODE_T * p)
573 {
574 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
575 gsl_vector *v = pop_vector (p, A68G_TRUE);
576 gsl_vector *u = pop_vector (p, A68G_TRUE);
577 ASSERT_GSL (gsl_vector_add (u, v));
578 push_vector (p, u);
579 gsl_vector_free (u);
580 gsl_vector_free (v);
581 (void) gsl_set_error_handler (save_handler);
582 }
583
584 //! @brief OP - = ([] REAL, [] REAL) [] REAL
585
586 void genie_vector_sub (NODE_T * p)
587 {
588 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
589 gsl_vector *v = pop_vector (p, A68G_TRUE);
590 gsl_vector *u = pop_vector (p, A68G_TRUE);
591 ASSERT_GSL (gsl_vector_sub (u, v));
592 push_vector (p, u);
593 gsl_vector_free (u);
594 gsl_vector_free (v);
595 (void) gsl_set_error_handler (save_handler);
596 }
597
598 //! @brief OP = = ([] REAL, [] REAL) BOOL
599
600 void genie_vector_eq (NODE_T * p)
601 {
602 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
603 gsl_vector *v = pop_vector (p, A68G_TRUE);
604 gsl_vector *u = pop_vector (p, A68G_TRUE);
605 ASSERT_GSL (gsl_vector_sub (u, v));
606 PUSH_VALUE (p, (BOOL_T) (gsl_vector_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
607 gsl_vector_free (u);
608 gsl_vector_free (v);
609 (void) gsl_set_error_handler (save_handler);
610 }
611
612 //! @brief OP /= = ([] REAL, [] REAL) BOOL
613
614 void genie_vector_ne (NODE_T * p)
615 {
616 genie_vector_eq (p);
617 genie_not_bool (p);
618 }
619
620 //! @brief OP +:= = (REF [] REAL, [] REAL) REF [] REAL
621
622 void genie_vector_plusab (NODE_T * p)
623 {
624 op_ab_torrix (p, M_REF_ROW_REAL, M_ROW_REAL, genie_vector_add);
625 }
626
627 //! @brief OP -:= = (REF [] REAL, [] REAL) REF [] REAL
628
629 void genie_vector_minusab (NODE_T * p)
630 {
631 op_ab_torrix (p, M_REF_ROW_REAL, M_ROW_REAL, genie_vector_sub);
632 }
633
634 //! @brief OP + = ([, ] REAL, [, ] REAL) [, ] REAL
635
636 void genie_matrix_add (NODE_T * p)
637 {
638 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
639 gsl_matrix *v = pop_matrix (p, A68G_TRUE);
640 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
641 ASSERT_GSL (gsl_matrix_add (u, v));
642 push_matrix (p, u);
643 gsl_matrix_free (u);
644 gsl_matrix_free (v);
645 (void) gsl_set_error_handler (save_handler);
646 }
647
648 //! @brief OP - = ([, ] REAL, [, ] REAL) [, ] REAL
649
650 void genie_matrix_sub (NODE_T * p)
651 {
652 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
653 gsl_matrix *v = pop_matrix (p, A68G_TRUE);
654 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
655 ASSERT_GSL (gsl_matrix_sub (u, v));
656 push_matrix (p, u);
657 gsl_matrix_free (u);
658 gsl_matrix_free (v);
659 (void) gsl_set_error_handler (save_handler);
660 }
661
662 //! @brief OP = = ([, ] REAL, [, ] REAL) [, ] BOOL
663
664 void genie_matrix_eq (NODE_T * p)
665 {
666 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
667 gsl_matrix *v = pop_matrix (p, A68G_TRUE);
668 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
669 ASSERT_GSL (gsl_matrix_sub (u, v));
670 PUSH_VALUE (p, (BOOL_T) (gsl_matrix_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
671 gsl_matrix_free (u);
672 gsl_matrix_free (v);
673 (void) gsl_set_error_handler (save_handler);
674 }
675
676 //! @brief OP /= = ([, ] REAL, [, ] REAL) [, ] BOOL
677
678 void genie_matrix_ne (NODE_T * p)
679 {
680 genie_matrix_eq (p);
681 genie_not_bool (p);
682 }
683
684 //! @brief OP +:= = (REF [, ] REAL, [, ] REAL) [, ] REAL
685
686 void genie_matrix_plusab (NODE_T * p)
687 {
688 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_ROW_ROW_REAL, genie_matrix_add);
689 }
690
691 //! @brief OP -:= = (REF [, ] REAL, [, ] REAL) [, ] REAL
692
693 void genie_matrix_minusab (NODE_T * p)
694 {
695 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_ROW_ROW_REAL, genie_matrix_sub);
696 }
697
698 //! @brief OP + = ([] COMPLEX, [] COMPLEX) [] COMPLEX
699
700 void genie_vector_complex_add (NODE_T * p)
701 {
702 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
703 gsl_complex one;
704 GSL_SET_COMPLEX (&one, 1.0, 0.0);
705 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
706 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
707 ASSERT_GSL (gsl_blas_zaxpy (one, u, v));
708 push_vector_complex (p, v);
709 gsl_vector_complex_free (u);
710 gsl_vector_complex_free (v);
711 (void) gsl_set_error_handler (save_handler);
712 }
713
714 //! @brief OP - = ([] COMPLEX, [] COMPLEX) [] COMPLEX
715
716 void genie_vector_complex_sub (NODE_T * p)
717 {
718 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
719 gsl_complex one;
720 GSL_SET_COMPLEX (&one, -1.0, 0.0);
721 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
722 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
723 ASSERT_GSL (gsl_blas_zaxpy (one, v, u));
724 push_vector_complex (p, u);
725 gsl_vector_complex_free (u);
726 gsl_vector_complex_free (v);
727 (void) gsl_set_error_handler (save_handler);
728 }
729
730 //! @brief OP = = ([] COMPLEX, [] COMPLEX) BOOL
731
732 void genie_vector_complex_eq (NODE_T * p)
733 {
734 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
735 gsl_complex one;
736 GSL_SET_COMPLEX (&one, -1.0, 0.0);
737 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
738 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
739 ASSERT_GSL (gsl_blas_zaxpy (one, v, u));
740 PUSH_VALUE (p, (BOOL_T) (gsl_vector_complex_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
741 gsl_vector_complex_free (u);
742 gsl_vector_complex_free (v);
743 (void) gsl_set_error_handler (save_handler);
744 }
745
746 //! @brief OP /= = ([] COMPLEX, [] COMPLEX) BOOL
747
748 void genie_vector_complex_ne (NODE_T * p)
749 {
750 genie_vector_complex_eq (p);
751 genie_not_bool (p);
752 }
753
754 //! @brief OP +:= = (REF [] COMPLEX, [] COMPLEX) [] COMPLEX
755
756 void genie_vector_complex_plusab (NODE_T * p)
757 {
758 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_ROW_COMPLEX, genie_vector_complex_add);
759 }
760
761 //! @brief OP -:= = (REF [] COMPLEX, [] COMPLEX) [] COMPLEX
762
763 void genie_vector_complex_minusab (NODE_T * p)
764 {
765 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_ROW_COMPLEX, genie_vector_complex_sub);
766 }
767
768 //! @brief OP + = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
769
770 void genie_matrix_complex_add (NODE_T * p)
771 {
772 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
773 gsl_matrix_complex *v = pop_matrix_complex (p, A68G_TRUE);
774 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
775 ASSERT_GSL (gsl_matrix_complex_add (u, v));
776 push_matrix_complex (p, u);
777 gsl_matrix_complex_free (u);
778 gsl_matrix_complex_free (v);
779 (void) gsl_set_error_handler (save_handler);
780 }
781
782 //! @brief OP - = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
783
784 void genie_matrix_complex_sub (NODE_T * p)
785 {
786 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
787 gsl_matrix_complex *v = pop_matrix_complex (p, A68G_TRUE);
788 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
789 ASSERT_GSL (gsl_matrix_complex_sub (u, v));
790 push_matrix_complex (p, u);
791 gsl_matrix_complex_free (u);
792 gsl_matrix_complex_free (v);
793 (void) gsl_set_error_handler (save_handler);
794 }
795
796 //! @brief OP = = ([, ] COMPLEX, [, ] COMPLEX) BOOL
797
798 void genie_matrix_complex_eq (NODE_T * p)
799 {
800 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
801 gsl_matrix_complex *v = pop_matrix_complex (p, A68G_TRUE);
802 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
803 ASSERT_GSL (gsl_matrix_complex_sub (u, v));
804 PUSH_VALUE (p, (BOOL_T) (gsl_matrix_complex_isnull (u) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
805 gsl_matrix_complex_free (u);
806 gsl_matrix_complex_free (v);
807 (void) gsl_set_error_handler (save_handler);
808 }
809
810 //! @brief OP /= = ([, ] COMPLEX, [, ] COMPLEX) BOOL
811
812 void genie_matrix_complex_ne (NODE_T * p)
813 {
814 genie_matrix_complex_eq (p);
815 genie_not_bool (p);
816 }
817
818 //! @brief OP +:= = (REF [, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
819
820 void genie_matrix_complex_plusab (NODE_T * p)
821 {
822 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_ROW_ROW_COMPLEX, genie_matrix_complex_add);
823 }
824
825 //! @brief OP -:= = (REF [, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
826
827 void genie_matrix_complex_minusab (NODE_T * p)
828 {
829 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_ROW_ROW_COMPLEX, genie_matrix_complex_sub);
830 }
831
832 //! @brief OP * = ([] REAL, REAL) [] REAL
833
834 void genie_vector_scale_real (NODE_T * p)
835 {
836 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
837 A68G_REAL v;
838 POP_OBJECT (p, &v, A68G_REAL);
839 gsl_vector *u = pop_vector (p, A68G_TRUE);
840 ASSERT_GSL (gsl_vector_scale (u, VALUE (&v)));
841 push_vector (p, u);
842 gsl_vector_free (u);
843 (void) gsl_set_error_handler (save_handler);
844 }
845
846 //! @brief OP * = (REAL, [] REAL) [] REAL
847
848 void genie_real_scale_vector (NODE_T * p)
849 {
850 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
851 gsl_vector *u = pop_vector (p, A68G_TRUE);
852 A68G_REAL v;
853 POP_OBJECT (p, &v, A68G_REAL);
854 ASSERT_GSL (gsl_vector_scale (u, VALUE (&v)));
855 push_vector (p, u);
856 gsl_vector_free (u);
857 (void) gsl_set_error_handler (save_handler);
858 }
859
860 //! @brief OP * = ([, ] REAL, REAL) [, ] REAL
861
862 void genie_matrix_scale_real (NODE_T * p)
863 {
864 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
865 A68G_REAL v;
866 POP_OBJECT (p, &v, A68G_REAL);
867 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
868 ASSERT_GSL (gsl_matrix_scale (u, VALUE (&v)));
869 push_matrix (p, u);
870 gsl_matrix_free (u);
871 (void) gsl_set_error_handler (save_handler);
872 }
873
874 //! @brief OP * = (REAL, [, ] REAL) [, ] REAL
875
876 void genie_real_scale_matrix (NODE_T * p)
877 {
878 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
879 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
880 A68G_REAL v;
881 POP_OBJECT (p, &v, A68G_REAL);
882 ASSERT_GSL (gsl_matrix_scale (u, VALUE (&v)));
883 push_matrix (p, u);
884 gsl_matrix_free (u);
885 (void) gsl_set_error_handler (save_handler);
886 }
887
888 //! @brief OP * = ([] COMPLEX, COMPLEX) [] COMPLEX
889
890 void genie_vector_complex_scale_complex (NODE_T * p)
891 {
892 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
893 A68G_REAL re, im; gsl_complex v;
894 POP_OBJECT (p, &im, A68G_REAL);
895 POP_OBJECT (p, &re, A68G_REAL);
896 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
897 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
898 gsl_blas_zscal (v, u);
899 push_vector_complex (p, u);
900 gsl_vector_complex_free (u);
901 (void) gsl_set_error_handler (save_handler);
902 }
903
904 //! @brief OP * = (COMPLEX, [] COMPLEX) [] COMPLEX
905
906 void genie_complex_scale_vector_complex (NODE_T * p)
907 {
908 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
909 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
910 A68G_REAL re, im; gsl_complex v;
911 POP_OBJECT (p, &im, A68G_REAL);
912 POP_OBJECT (p, &re, A68G_REAL);
913 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
914 gsl_blas_zscal (v, u);
915 push_vector_complex (p, u);
916 gsl_vector_complex_free (u);
917 (void) gsl_set_error_handler (save_handler);
918 }
919
920 //! @brief OP * = ([, ] COMPLEX, COMPLEX) [, ] COMPLEX
921
922 void genie_matrix_complex_scale_complex (NODE_T * p)
923 {
924 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
925 A68G_REAL re, im; gsl_complex v;
926 POP_OBJECT (p, &im, A68G_REAL);
927 POP_OBJECT (p, &re, A68G_REAL);
928 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
929 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
930 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
931 push_matrix_complex (p, u);
932 gsl_matrix_complex_free (u);
933 (void) gsl_set_error_handler (save_handler);
934 }
935
936 //! @brief OP * = (COMPLEX, [, ] COMPLEX) [, ] COMPLEX
937
938 void genie_complex_scale_matrix_complex (NODE_T * p)
939 {
940 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
941 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
942 A68G_REAL re, im; gsl_complex v;
943 POP_OBJECT (p, &im, A68G_REAL);
944 POP_OBJECT (p, &re, A68G_REAL);
945 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
946 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
947 push_matrix_complex (p, u);
948 gsl_matrix_complex_free (u);
949 (void) gsl_set_error_handler (save_handler);
950 }
951
952 //! @brief OP *:= (REF [] REAL, REAL) REF [] REAL
953
954 void genie_vector_scale_real_ab (NODE_T * p)
955 {
956 op_ab_torrix (p, M_REF_ROW_REAL, M_REAL, genie_vector_scale_real);
957 }
958
959 //! @brief OP *:= (REF [, ] REAL, REAL) REF [, ] REAL
960
961 void genie_matrix_scale_real_ab (NODE_T * p)
962 {
963 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_REAL, genie_matrix_scale_real);
964 }
965
966 //! @brief OP *:= (REF [] COMPLEX, COMPLEX) REF [] COMPLEX
967
968 void genie_vector_complex_scale_complex_ab (NODE_T * p)
969 {
970 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_COMPLEX, genie_vector_complex_scale_complex);
971 }
972
973 //! @brief OP *:= (REF [, ] COMPLEX, COMPLEX) REF [, ] COMPLEX
974
975 void genie_matrix_complex_scale_complex_ab (NODE_T * p)
976 {
977 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_COMPLEX, genie_matrix_complex_scale_complex);
978 }
979
980 //! @brief OP / = ([] REAL, REAL) [] REAL
981
982 void genie_vector_div_real (NODE_T * p)
983 {
984 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
985 A68G_REAL v;
986 POP_OBJECT (p, &v, A68G_REAL);
987 if (VALUE (&v) == 0.0) {
988 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_REAL);
989 exit_genie (p, A68G_RUNTIME_ERROR);
990 }
991 gsl_vector *u = pop_vector (p, A68G_TRUE);
992 ASSERT_GSL (gsl_vector_scale (u, 1.0 / VALUE (&v)));
993 push_vector (p, u);
994 gsl_vector_free (u);
995 (void) gsl_set_error_handler (save_handler);
996 }
997
998 //! @brief OP / = ([, ] REAL, REAL) [, ] REAL
999
1000 void genie_matrix_div_real (NODE_T * p)
1001 {
1002 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1003 A68G_REAL v;
1004 POP_OBJECT (p, &v, A68G_REAL);
1005 if (VALUE (&v) == 0.0) {
1006 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_ROW_REAL);
1007 exit_genie (p, A68G_RUNTIME_ERROR);
1008 }
1009 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
1010 ASSERT_GSL (gsl_matrix_scale (u, 1.0 / VALUE (&v)));
1011 push_matrix (p, u);
1012 gsl_matrix_free (u);
1013 (void) gsl_set_error_handler (save_handler);
1014 }
1015
1016 //! @brief OP / = ([] COMPLEX, COMPLEX) [] COMPLEX
1017
1018 void genie_vector_complex_div_complex (NODE_T * p)
1019 {
1020 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1021 A68G_REAL re, im; gsl_complex v;
1022 POP_OBJECT (p, &im, A68G_REAL);
1023 POP_OBJECT (p, &re, A68G_REAL);
1024 if (VALUE (&re) == 0.0 && VALUE (&im) == 0.0) {
1025 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_COMPLEX);
1026 exit_genie (p, A68G_RUNTIME_ERROR);
1027 }
1028 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
1029 v = gsl_complex_inverse (v);
1030 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1031 gsl_blas_zscal (v, u);
1032 push_vector_complex (p, u);
1033 gsl_vector_complex_free (u);
1034 (void) gsl_set_error_handler (save_handler);
1035 }
1036
1037 //! @brief OP / = ([, ] COMPLEX, COMPLEX) [, ] COMPLEX
1038
1039 void genie_matrix_complex_div_complex (NODE_T * p)
1040 {
1041 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1042 A68G_REAL re, im; gsl_complex v;
1043 POP_OBJECT (p, &im, A68G_REAL);
1044 POP_OBJECT (p, &re, A68G_REAL);
1045 if (VALUE (&re) == 0.0 && VALUE (&im) == 0.0) {
1046 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_DIVISION_BY_ZERO, M_ROW_ROW_COMPLEX);
1047 exit_genie (p, A68G_RUNTIME_ERROR);
1048 }
1049 GSL_SET_COMPLEX (&v, VALUE (&re), VALUE (&im));
1050 v = gsl_complex_inverse (v);
1051 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
1052 ASSERT_GSL (gsl_matrix_complex_scale (u, v));
1053 push_matrix_complex (p, u);
1054 gsl_matrix_complex_free (u);
1055 (void) gsl_set_error_handler (save_handler);
1056 }
1057
1058 //! @brief OP /:= (REF [] REAL, REAL) REF [] REAL
1059
1060 void genie_vector_div_real_ab (NODE_T * p)
1061 {
1062 op_ab_torrix (p, M_REF_ROW_REAL, M_REAL, genie_vector_div_real);
1063 }
1064
1065 //! @brief OP /:= (REF [, ] REAL, REAL) REF [, ] REAL
1066
1067 void genie_matrix_div_real_ab (NODE_T * p)
1068 {
1069 op_ab_torrix (p, M_REF_ROW_ROW_REAL, M_REAL, genie_matrix_div_real);
1070 }
1071
1072 //! @brief OP /:= (REF [] COMPLEX, COMPLEX) REF [] COMPLEX
1073
1074 void genie_vector_complex_div_complex_ab (NODE_T * p)
1075 {
1076 op_ab_torrix (p, M_REF_ROW_COMPLEX, M_COMPLEX, genie_vector_complex_div_complex);
1077 }
1078
1079 //! @brief OP /:= (REF [, ] COMPLEX, COMPLEX) REF [, ] COMPLEX
1080
1081 void genie_matrix_complex_div_complex_ab (NODE_T * p)
1082 {
1083 op_ab_torrix (p, M_REF_ROW_ROW_COMPLEX, M_COMPLEX, genie_matrix_complex_div_complex);
1084 }
1085
1086 //! @brief OP * = ([] REAL, [] REAL) REAL
1087
1088 void genie_vector_dot (NODE_T * p)
1089 {
1090 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1091 gsl_vector *v = pop_vector (p, A68G_TRUE);
1092 gsl_vector *u = pop_vector (p, A68G_TRUE);
1093 REAL_T w;
1094 ASSERT_GSL (gsl_blas_ddot (u, v, &w));
1095 PUSH_VALUE (p, w, A68G_REAL);
1096 gsl_vector_free (u);
1097 gsl_vector_free (v);
1098 (void) gsl_set_error_handler (save_handler);
1099 }
1100
1101 //! @brief OP * = ([] COMPLEX, [] COMPLEX) COMPLEX
1102
1103 void genie_vector_complex_dot (NODE_T * p)
1104 {
1105 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1106 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
1107 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1108 gsl_complex w;
1109 ASSERT_GSL (gsl_blas_zdotc (u, v, &w));
1110 PUSH_VALUE (p, GSL_REAL (w), A68G_REAL);
1111 PUSH_VALUE (p, GSL_IMAG (w), A68G_REAL);
1112 gsl_vector_complex_free (u);
1113 gsl_vector_complex_free (v);
1114 (void) gsl_set_error_handler (save_handler);
1115 }
1116
1117 //! @brief OP NORM = ([] REAL) REAL
1118
1119 void genie_vector_norm (NODE_T * p)
1120 {
1121 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1122 gsl_vector *u = pop_vector (p, A68G_TRUE);
1123 PUSH_VALUE (p, gsl_blas_dnrm2 (u), A68G_REAL);
1124 gsl_vector_free (u);
1125 (void) gsl_set_error_handler (save_handler);
1126 }
1127
1128 //! @brief OP NORM = ([] COMPLEX) COMPLEX
1129
1130 void genie_vector_complex_norm (NODE_T * p)
1131 {
1132 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1133 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1134 PUSH_VALUE (p, gsl_blas_dznrm2 (u), A68G_REAL);
1135 gsl_vector_complex_free (u);
1136 (void) gsl_set_error_handler (save_handler);
1137 }
1138
1139 //! @brief OP DYAD = ([] REAL, [] REAL) [, ] REAL
1140
1141 void genie_vector_dyad (NODE_T * p)
1142 {
1143 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1144 gsl_vector *v = pop_vector (p, A68G_TRUE);
1145 gsl_vector *u = pop_vector (p, A68G_TRUE);
1146 size_t len1 = SIZE (u), len2 = SIZE (v);
1147 gsl_matrix *w = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
1148 for (int j = 0; j < len1; j++) {
1149 REAL_T uj = gsl_vector_get (u, (size_t) j);
1150 for (int k = 0; k < len2; k++) {
1151 REAL_T vk = gsl_vector_get (v, (size_t) k);
1152 gsl_matrix_set (w, (size_t) j, (size_t) k, uj * vk);
1153 }
1154 }
1155 push_matrix (p, w);
1156 gsl_vector_free (u);
1157 gsl_vector_free (v);
1158 gsl_matrix_free (w);
1159 (void) gsl_set_error_handler (save_handler);
1160 }
1161
1162 //! @brief OP DYAD = ([] COMPLEX, [] COMPLEX) [, ] COMPLEX
1163
1164 void genie_vector_complex_dyad (NODE_T * p)
1165 {
1166 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1167 gsl_vector_complex *v = pop_vector_complex (p, A68G_TRUE);
1168 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1169 size_t len1 = SIZE (u), len2 = SIZE (v);
1170 gsl_matrix_complex *w = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
1171 for (int j = 0; j < len1; j++) {
1172 gsl_complex uj = gsl_vector_complex_get (u, (size_t) j);
1173 for (int k = 0; k < len2; k++) {
1174 gsl_complex vk = gsl_vector_complex_get (u, (size_t) k);
1175 gsl_matrix_complex_set (w, (size_t) j, (size_t) k, gsl_complex_mul (uj, vk));
1176 }
1177 }
1178 push_matrix_complex (p, w);
1179 gsl_vector_complex_free (u);
1180 gsl_vector_complex_free (v);
1181 gsl_matrix_complex_free (w);
1182 (void) gsl_set_error_handler (save_handler);
1183 }
1184
1185 //! @brief OP * = ([, ] REAL, [] REAL) [] REAL
1186
1187 void genie_matrix_times_vector (NODE_T * p)
1188 {
1189 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1190 gsl_vector *u = pop_vector (p, A68G_TRUE);
1191 gsl_matrix *w = pop_matrix (p, A68G_TRUE);
1192 size_t len = SIZE (u);
1193 gsl_vector *v = gsl_vector_calloc ((size_t) len);
1194 gsl_vector_set_zero (v);
1195 ASSERT_GSL (gsl_blas_dgemv (CblasNoTrans, 1.0, w, u, 0.0, v));
1196 push_vector (p, v);
1197 gsl_vector_free (u);
1198 gsl_vector_free (v);
1199 gsl_matrix_free (w);
1200 (void) gsl_set_error_handler (save_handler);
1201 }
1202
1203 //! @brief OP * = ([] REAL, [, ] REAL) [] REAL
1204
1205 void genie_vector_times_matrix (NODE_T * p)
1206 {
1207 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1208 gsl_matrix *w = pop_matrix (p, A68G_TRUE);
1209 ASSERT_GSL (gsl_matrix_transpose (w));
1210 gsl_vector *u = pop_vector (p, A68G_TRUE);
1211 size_t len = SIZE (u);
1212 gsl_vector *v = gsl_vector_calloc ((size_t) len);
1213 gsl_vector_set_zero (v);
1214 ASSERT_GSL (gsl_blas_dgemv (CblasNoTrans, 1.0, w, u, 0.0, v));
1215 push_vector (p, v);
1216 gsl_vector_free (u);
1217 gsl_vector_free (v);
1218 gsl_matrix_free (w);
1219 (void) gsl_set_error_handler (save_handler);
1220 }
1221
1222 //! @brief OP * = ([, ] REAL, [, ] REAL) [, ] REAL
1223
1224 void genie_matrix_times_matrix (NODE_T * p)
1225 {
1226 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1227 gsl_matrix *v = pop_matrix (p, A68G_TRUE);
1228 gsl_matrix *u = pop_matrix (p, A68G_TRUE);
1229 size_t len2 = SIZE2 (v), len1 = SIZE1 (u);
1230 gsl_matrix *w = gsl_matrix_calloc ((size_t) len1, (size_t) len2);
1231 ASSERT_GSL (gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, u, v, 0.0, w));
1232 push_matrix (p, w);
1233 gsl_matrix_free (u);
1234 gsl_matrix_free (v);
1235 gsl_matrix_free (w);
1236 (void) gsl_set_error_handler (save_handler);
1237 }
1238
1239 //! @brief OP * = ([, ] COMPLEX, [] COMPLEX) [] COMPLEX
1240
1241 void genie_matrix_complex_times_vector (NODE_T * p)
1242 {
1243 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1244 gsl_complex zero, one;
1245 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1246 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1247 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1248 gsl_matrix_complex *w = pop_matrix_complex (p, A68G_TRUE);
1249 size_t len = SIZE (u);
1250 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
1251 ASSERT_GSL (gsl_blas_zgemv (CblasNoTrans, one, w, u, zero, v));
1252 push_vector_complex (p, v);
1253 gsl_vector_complex_free (u);
1254 gsl_vector_complex_free (v);
1255 gsl_matrix_complex_free (w);
1256 (void) gsl_set_error_handler (save_handler);
1257 }
1258
1259 //! @brief OP * = ([] COMPLEX, [, ] COMPLEX) [] COMPLEX
1260
1261 void genie_vector_complex_times_matrix (NODE_T * p)
1262 {
1263 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1264 gsl_complex zero, one;
1265 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1266 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1267 gsl_matrix_complex *w = pop_matrix_complex (p, A68G_TRUE);
1268 ASSERT_GSL (gsl_matrix_complex_transpose (w));
1269 gsl_vector_complex *u = pop_vector_complex (p, A68G_TRUE);
1270 size_t len = SIZE (u);
1271 gsl_vector_complex *v = gsl_vector_complex_calloc ((size_t) len);
1272 ASSERT_GSL (gsl_blas_zgemv (CblasNoTrans, one, w, u, zero, v));
1273 push_vector_complex (p, v);
1274 gsl_vector_complex_free (u);
1275 gsl_vector_complex_free (v);
1276 gsl_matrix_complex_free (w);
1277 (void) gsl_set_error_handler (save_handler);
1278 }
1279
1280 //! @brief OP * = ([, ] COMPLEX, [, ] COMPLEX) [, ] COMPLEX
1281
1282 void genie_matrix_complex_times_matrix (NODE_T * p)
1283 {
1284 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
1285 gsl_complex zero, one;
1286 GSL_SET_COMPLEX (&zero, 0.0, 0.0);
1287 GSL_SET_COMPLEX (&one, 1.0, 0.0);
1288 gsl_matrix_complex *v = pop_matrix_complex (p, A68G_TRUE);
1289 gsl_matrix_complex *u = pop_matrix_complex (p, A68G_TRUE);
1290 size_t len1 = SIZE2 (v), len2 = SIZE1 (u);
1291 gsl_matrix_complex *w = gsl_matrix_complex_calloc ((size_t) len1, (size_t) len2);
1292 ASSERT_GSL (gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, one, u, v, zero, w));
1293 gsl_matrix_complex_free (u);
1294 gsl_matrix_complex_free (v);
1295 push_matrix_complex (p, w);
1296 gsl_matrix_complex_free (w);
1297 (void) gsl_set_error_handler (save_handler);
1298 }
1299
1300 #endif
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
