single-torrix-gsl.c

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