single-decomposition.c

     
   1  //! @file single-decomposition.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 LU, QR and Choleski decomposition.
  25  
  26  #include "a68g.h"
  27  #include "a68g-torrix.h"
  28  
  29  #if defined (HAVE_GSL)
  30  
  31  //! @brief PROC lu decomp = ([, ] REAL, REF [] INT, REF INT) [, ] REAL
  32  
  33  void genie_matrix_lu (NODE_T * p)
  34  {
  35    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  36    torrix_error_node = p;
  37    A68_REF ref_signum;
  38    POP_REF (p, &ref_signum);
  39    CHECK_REF (p, ref_signum, M_REF_INT);
  40    A68_REF ref_q;
  41    POP_REF (p, &ref_q);
  42    CHECK_REF (p, ref_q, M_REF_ROW_INT);
  43    PUSH_REF (p, *DEREF (A68_ROW, &ref_q));
  44    gsl_permutation *q = pop_permutation (p, A68_FALSE);
  45    gsl_matrix *u = pop_matrix (p, A68_TRUE);
  46    int sign;
  47    ASSERT_GSL (gsl_linalg_LU_decomp (u, q, &sign));
  48    A68_INT signum;
  49    VALUE (&signum) = sign;
  50    STATUS (&signum) = INIT_MASK;
  51    *DEREF (A68_INT, &ref_signum) = signum;
  52    push_permutation (p, q);
  53    POP_REF (p, DEREF (A68_ROW, &ref_q));
  54    push_matrix (p, u);
  55    gsl_matrix_free (u);
  56    gsl_permutation_free (q);
  57    (void) gsl_set_error_handler (save_handler);
  58  }
  59  
  60  //! @brief PROC lu det = ([, ] REAL, INT) REAL
  61  
  62  void genie_matrix_lu_det (NODE_T * p)
  63  {
  64    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  65    torrix_error_node = p;
  66    A68_INT signum;
  67    POP_OBJECT (p, &signum, A68_INT);
  68    gsl_matrix *lu = pop_matrix (p, A68_TRUE);
  69    PUSH_VALUE (p, gsl_linalg_LU_det (lu, VALUE (&signum)), A68_REAL);
  70    gsl_matrix_free (lu);
  71    (void) gsl_set_error_handler (save_handler);
  72  }
  73  
  74  //! @brief PROC lu inv = ([, ] REAL, [] INT) [, ] REAL
  75  
  76  void genie_matrix_lu_inv (NODE_T * p)
  77  {
  78    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  79    torrix_error_node = p;
  80    gsl_permutation *q = pop_permutation (p, A68_TRUE);
  81    gsl_matrix *lu = pop_matrix (p, A68_TRUE);
  82    gsl_matrix *inv = gsl_matrix_calloc (SIZE1 (lu), SIZE2 (lu));
  83    ASSERT_GSL (gsl_linalg_LU_invert (lu, q, inv));
  84    push_matrix (p, inv);
  85    gsl_matrix_free (lu);
  86    gsl_matrix_free (inv);
  87    gsl_permutation_free (q);
  88    (void) gsl_set_error_handler (save_handler);
  89  }
  90  
  91  //! @brief PROC lu solve ([, ] REAL, [, ] REAL, [] INT, [] REAL) [] REAL
  92  
  93  void genie_matrix_lu_solve (NODE_T * p)
  94  {
  95    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  96    torrix_error_node = p;
  97    gsl_vector *b = pop_vector (p, A68_TRUE);
  98    gsl_permutation *q = pop_permutation (p, A68_TRUE);
  99    gsl_matrix *lu = pop_matrix (p, A68_TRUE);
 100    gsl_matrix *a = pop_matrix (p, A68_TRUE);
 101    gsl_vector *x = gsl_vector_calloc (SIZE (b));
 102    gsl_vector *r = gsl_vector_calloc (SIZE (b));
 103    ASSERT_GSL (gsl_linalg_LU_solve (lu, q, b, x));
 104    ASSERT_GSL (gsl_linalg_LU_refine (a, lu, q, b, x, r));
 105    push_vector (p, x);
 106    gsl_matrix_free (a);
 107    gsl_matrix_free (lu);
 108    gsl_vector_free (b);
 109    gsl_vector_free (r);
 110    gsl_vector_free (x);
 111    gsl_permutation_free (q);
 112    (void) gsl_set_error_handler (save_handler);
 113  }
 114  
 115  //! @brief PROC complex lu decomp = ([, ] COMPLEX, REF [] INT, REF INT) [, ] COMPLEX
 116  
 117  void genie_matrix_complex_lu (NODE_T * p)
 118  {
 119    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 120    A68_REF ref_signum, ref_q;
 121    gsl_permutation *q;
 122    gsl_matrix_complex *u;
 123    int sign;
 124    A68_INT signum;
 125    torrix_error_node = p;
 126    POP_REF (p, &ref_signum);
 127    CHECK_REF (p, ref_signum, M_REF_INT);
 128    POP_REF (p, &ref_q);
 129    CHECK_REF (p, ref_q, M_REF_ROW_INT);
 130    PUSH_REF (p, *DEREF (A68_ROW, &ref_q));
 131    q = pop_permutation (p, A68_FALSE);
 132    u = pop_matrix_complex (p, A68_TRUE);
 133    ASSERT_GSL (gsl_linalg_complex_LU_decomp (u, q, &sign));
 134    VALUE (&signum) = sign;
 135    STATUS (&signum) = INIT_MASK;
 136    *DEREF (A68_INT, &ref_signum) = signum;
 137    push_permutation (p, q);
 138    POP_REF (p, DEREF (A68_ROW, &ref_q));
 139    push_matrix_complex (p, u);
 140    gsl_matrix_complex_free (u);
 141    gsl_permutation_free (q);
 142    (void) gsl_set_error_handler (save_handler);
 143  }
 144  
 145  //! @brief PROC complex lu det = ([, ] COMPLEX, INT) COMPLEX
 146  
 147  void genie_matrix_complex_lu_det (NODE_T * p)
 148  {
 149    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 150    gsl_matrix_complex *lu;
 151    A68_INT signum;
 152    gsl_complex det;
 153    torrix_error_node = p;
 154    POP_OBJECT (p, &signum, A68_INT);
 155    lu = pop_matrix_complex (p, A68_TRUE);
 156    det = gsl_linalg_complex_LU_det (lu, VALUE (&signum));
 157    PUSH_VALUE (p, GSL_REAL (det), A68_REAL);
 158    PUSH_VALUE (p, GSL_IMAG (det), A68_REAL);
 159    gsl_matrix_complex_free (lu);
 160    (void) gsl_set_error_handler (save_handler);
 161  }
 162  
 163  //! @brief PROC complex lu inv = ([, ] COMPLEX, [] INT) [, ] COMPLEX
 164  
 165  void genie_matrix_complex_lu_inv (NODE_T * p)
 166  {
 167    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 168    gsl_permutation *q;
 169    gsl_matrix_complex *lu, *inv;
 170    torrix_error_node = p;
 171    q = pop_permutation (p, A68_TRUE);
 172    lu = pop_matrix_complex (p, A68_TRUE);
 173    inv = gsl_matrix_complex_calloc (SIZE1 (lu), SIZE2 (lu));
 174    ASSERT_GSL (gsl_linalg_complex_LU_invert (lu, q, inv));
 175    push_matrix_complex (p, inv);
 176    gsl_matrix_complex_free (lu);
 177    gsl_matrix_complex_free (inv);
 178    gsl_permutation_free (q);
 179    (void) gsl_set_error_handler (save_handler);
 180  }
 181  
 182  //! @brief PROC complex lu solve ([, ] COMPLEX, [, ] COMPLEX, [] INT, [] COMPLEX) [] COMPLEX
 183  
 184  void genie_matrix_complex_lu_solve (NODE_T * p)
 185  {
 186    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 187    torrix_error_node = p;
 188    gsl_vector_complex *b = pop_vector_complex (p, A68_TRUE);
 189    gsl_permutation *q = pop_permutation (p, A68_TRUE);
 190    gsl_matrix_complex *lu = pop_matrix_complex (p, A68_TRUE);
 191    gsl_matrix_complex *a = pop_matrix_complex (p, A68_TRUE);
 192    gsl_vector_complex *x = gsl_vector_complex_calloc (SIZE (b));
 193    gsl_vector_complex *r = gsl_vector_complex_calloc (SIZE (b));
 194    ASSERT_GSL (gsl_linalg_complex_LU_solve (lu, q, b, x));
 195    ASSERT_GSL (gsl_linalg_complex_LU_refine (a, lu, q, b, x, r));
 196    gsl_matrix_complex_free (a);
 197    gsl_matrix_complex_free (lu);
 198    gsl_permutation_free (q);
 199    gsl_vector_complex_free (b);
 200    gsl_vector_complex_free (r);
 201    push_vector_complex (p, x);
 202    gsl_vector_complex_free (x);
 203    (void) gsl_set_error_handler (save_handler);
 204  }
 205  
 206  //! @brief PROC qr decomp = ([, ] REAL, [] REAL) [, ] REAL
 207  
 208  void genie_matrix_qr (NODE_T * p)
 209  {
 210    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 211    torrix_error_node = p;
 212    A68_REF ref_t;
 213    POP_REF (p, &ref_t);
 214    CHECK_REF (p, ref_t, M_REF_ROW_REAL);
 215    PUSH_REF (p, *DEREF (A68_ROW, &ref_t));
 216    gsl_vector *t = pop_vector (p, A68_FALSE);
 217    gsl_matrix *a = pop_matrix (p, A68_TRUE);
 218    ASSERT_GSL (gsl_linalg_QR_decomp (a, t));
 219    push_vector (p, t);
 220    gsl_vector_free (t);
 221    POP_REF (p, DEREF (A68_ROW, &ref_t));
 222    push_matrix (p, a);
 223    gsl_matrix_free (a);
 224    (void) gsl_set_error_handler (save_handler);
 225  }
 226  
 227  //! @brief PROC qr solve = ([, ] REAL, [] REAL, [] REAL) [] REAL
 228  
 229  void genie_matrix_qr_solve (NODE_T * p)
 230  {
 231    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 232    torrix_error_node = p;
 233    gsl_vector *b = pop_vector (p, A68_TRUE);
 234    gsl_vector *t = pop_vector (p, A68_TRUE);
 235    gsl_matrix *q = pop_matrix (p, A68_TRUE);
 236    gsl_vector *x = gsl_vector_calloc (SIZE (b));
 237    ASSERT_GSL (gsl_linalg_QR_solve (q, t, b, x));
 238    push_vector (p, x);
 239    gsl_vector_free (x);
 240    gsl_vector_free (b);
 241    gsl_vector_free (t);
 242    gsl_matrix_free (q);
 243    (void) gsl_set_error_handler (save_handler);
 244  }
 245  
 246  //! @brief PROC qr ls solve = ([, ] REAL, [] REAL, [] REAL) [] REAL
 247  
 248  void genie_matrix_qr_ls_solve (NODE_T * p)
 249  {
 250    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 251    torrix_error_node = p;
 252    gsl_vector *b = pop_vector (p, A68_TRUE);
 253    gsl_vector *t = pop_vector (p, A68_TRUE);
 254    gsl_matrix *q = pop_matrix (p, A68_TRUE);
 255    gsl_vector *r = gsl_vector_calloc (SIZE (b));
 256    gsl_vector *x = gsl_vector_calloc (SIZE (b));
 257    ASSERT_GSL (gsl_linalg_QR_lssolve (q, t, b, x, r));
 258    push_vector (p, x);
 259    gsl_vector_free (x);
 260    gsl_vector_free (r);
 261    gsl_vector_free (b);
 262    gsl_vector_free (t);
 263    gsl_matrix_free (q);
 264    (void) gsl_set_error_handler (save_handler);
 265  }
 266  
 267  //! @brief PROC cholesky decomp = ([, ] REAL) [, ] REAL
 268  
 269  void genie_matrix_ch (NODE_T * p)
 270  {
 271    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 272    torrix_error_node = p;
 273    gsl_matrix *a = pop_matrix (p, A68_TRUE);
 274    ASSERT_GSL (gsl_linalg_cholesky_decomp (a));
 275    push_matrix (p, a);
 276    gsl_matrix_free (a);
 277    (void) gsl_set_error_handler (save_handler);
 278  }
 279  
 280  //! @brief PROC cholesky solve = ([, ] REAL, [] REAL) [] REAL
 281  
 282  void genie_matrix_ch_solve (NODE_T * p)
 283  {
 284    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 285    gsl_matrix *c;
 286    gsl_vector *b, *x;
 287    torrix_error_node = p;
 288    b = pop_vector (p, A68_TRUE);
 289    c = pop_matrix (p, A68_TRUE);
 290    x = gsl_vector_calloc (SIZE (b));
 291    ASSERT_GSL (gsl_linalg_cholesky_solve (c, b, x));
 292    push_vector (p, x);
 293    gsl_vector_free (x);
 294    gsl_vector_free (b);
 295    gsl_matrix_free (c);
 296    (void) gsl_set_error_handler (save_handler);
 297  }
 298  
 299  #endif