Algol 68 Genie

single-svd.c

     
   1  //! @file single-svd.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-2026 J. Marcel van der Veer [algol68g@algol68genie.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 matrix SVD decomposition.
  25  
  26  #include "a68g.h"
  27  #include "a68g-torrix.h"
  28  
  29  #if defined (HAVE_GSL)
  30  
  31  //! @brief PROC svd decomp = ([, ] REAL, REF [, ] REAL, REF [] REAL, REF [, ] REAL) VOID
  32  
  33  void genie_matrix_svd (NODE_T * p)
  34  {
  35    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  36    A68G_REF ref_u, ref_s, ref_v;
  37    POP_REF (p, &ref_v);
  38    POP_REF (p, &ref_s);
  39    POP_REF (p, &ref_u);
  40    gsl_matrix *a = pop_matrix (p, A68G_TRUE), *u, *v; gsl_vector *s, *w;
  41    unt M = SIZE1 (a), N = SIZE2 (a);
  42  // GSL computes thin SVD, only handles M >= N. GSL returns V, not Vᵀ.
  43    if (M >= N) {
  44  // X = USVᵀ
  45      s = gsl_vector_calloc (N);
  46      u = gsl_matrix_calloc (M, N);
  47      v = gsl_matrix_calloc (N, N);
  48      w = gsl_vector_calloc (N);
  49      gsl_matrix_memcpy (u, a);
  50      ASSERT_GSL (gsl_linalg_SV_decomp (u, v, s, w));
  51      if (!IS_NIL (ref_u)) {
  52        *DEREF (A68G_REF, &ref_u) = matrix_to_row (p, u);
  53      }
  54      if (!IS_NIL (ref_s)) {
  55        *DEREF (A68G_REF, &ref_s) = vector_to_row (p, s);
  56      }
  57      if (!IS_NIL (ref_v)) {
  58        *DEREF (A68G_REF, &ref_v) = matrix_to_row (p, v);
  59      }
  60    } else {
  61  // Xᵀ = VSUᵀ
  62      s = gsl_vector_calloc (M);
  63      u = gsl_matrix_calloc (N, M);
  64      v = gsl_matrix_calloc (M, M);
  65      w = gsl_vector_calloc (M);
  66      gsl_matrix_transpose_memcpy (u, a);
  67      ASSERT_GSL (gsl_linalg_SV_decomp (u, v, s, w));
  68      if (!IS_NIL (ref_u)) {
  69        *DEREF (A68G_REF, &ref_u) = matrix_to_row (p, v);
  70      }
  71      if (!IS_NIL (ref_s)) {
  72        *DEREF (A68G_REF, &ref_s) = vector_to_row (p, s);
  73      }
  74      if (!IS_NIL (ref_v)) {
  75        *DEREF (A68G_REF, &ref_v) = matrix_to_row (p, u);
  76      }
  77    }
  78    gsl_matrix_free (a);
  79    gsl_matrix_free (u);
  80    gsl_matrix_free (v);
  81    gsl_vector_free (s);
  82    gsl_vector_free (w);
  83    (void) gsl_set_error_handler (save_handler);
  84  }
  85  
  86  //! @brief PROC svd solve = ([, ] REAL, [] REAL, [,] REAL, [] REAL) [] REAL
  87  
  88  void genie_matrix_svd_solve (NODE_T * p)
  89  {
  90    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  91    gsl_vector *b = pop_vector (p, A68G_TRUE);
  92    gsl_matrix *v = pop_matrix (p, A68G_TRUE);
  93    gsl_vector *s = pop_vector (p, A68G_TRUE);
  94    gsl_matrix *u = pop_matrix (p, A68G_TRUE);
  95    gsl_vector *x = gsl_vector_calloc (SIZE (b));
  96    ASSERT_GSL (gsl_linalg_SV_solve (u, v, s, b, x));
  97    push_vector (p, x);
  98    gsl_matrix_free (u);
  99    gsl_matrix_free (v);
 100    gsl_vector_free (b);
 101    gsl_vector_free (s);
 102    gsl_vector_free (x);
 103    (void) gsl_set_error_handler (save_handler);
 104  }
 105  
 106  #endif