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-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 matrix SVD decomposition using GSL.
  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    torrix_error_node = p;
  37  //
  38    A68_REF ref_u, ref_s, ref_v;
  39    POP_REF (p, &ref_v);
  40    POP_REF (p, &ref_s);
  41    POP_REF (p, &ref_u);
  42  //
  43    gsl_matrix *a = pop_matrix (p, A68_TRUE), *u, *v; gsl_vector *s, *w;
  44    unt M = SIZE1 (a), N = SIZE2 (a);
  45  // GSL computes thin SVD, only handles M >= N. GSL returns V, not Vᵀ.
  46    if (M >= N) {
  47  // X = USVᵀ
  48      s = gsl_vector_calloc (N);
  49      u = gsl_matrix_calloc (M, N);
  50      v = gsl_matrix_calloc (N, N);
  51      w = gsl_vector_calloc (N);
  52      gsl_matrix_memcpy (u, a);
  53      ASSERT_GSL (gsl_linalg_SV_decomp (u, v, s, w));
  54      if (!IS_NIL (ref_u)) {
  55        *DEREF (A68_REF, &ref_u) = matrix_to_row (p, u);
  56      }
  57      if (!IS_NIL (ref_s)) {
  58        *DEREF (A68_REF, &ref_s) = vector_to_row (p, s);
  59      }
  60      if (!IS_NIL (ref_v)) {
  61        *DEREF (A68_REF, &ref_v) = matrix_to_row (p, v);
  62      }
  63    } else {
  64  // Xᵀ = VSUᵀ
  65      s = gsl_vector_calloc (M);
  66      u = gsl_matrix_calloc (N, M);
  67      v = gsl_matrix_calloc (M, M);
  68      w = gsl_vector_calloc (M);
  69      gsl_matrix_transpose_memcpy (u, a);
  70      ASSERT_GSL (gsl_linalg_SV_decomp (u, v, s, w));
  71      if (!IS_NIL (ref_u)) {
  72        *DEREF (A68_REF, &ref_u) = matrix_to_row (p, v);
  73      }
  74      if (!IS_NIL (ref_s)) {
  75        *DEREF (A68_REF, &ref_s) = vector_to_row (p, s);
  76      }
  77      if (!IS_NIL (ref_v)) {
  78        *DEREF (A68_REF, &ref_v) = matrix_to_row (p, u);
  79      }
  80    }
  81    gsl_matrix_free (a);
  82    gsl_matrix_free (u);
  83    gsl_matrix_free (v);
  84    gsl_vector_free (s);
  85    gsl_vector_free (w);
  86    (void) gsl_set_error_handler (save_handler);
  87  }
  88  
  89  //! @brief PROC svd solve = ([, ] REAL, [] REAL, [,] REAL, [] REAL) [] REAL
  90  
  91  void genie_matrix_svd_solve (NODE_T * p)
  92  {
  93    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
  94    torrix_error_node = p;
  95    gsl_vector *b = pop_vector (p, A68_TRUE);
  96    gsl_matrix *v = pop_matrix (p, A68_TRUE);
  97    gsl_vector *s = pop_vector (p, A68_TRUE);
  98    gsl_matrix *u = pop_matrix (p, A68_TRUE);
  99    gsl_vector *x = gsl_vector_calloc (SIZE (b));
 100    ASSERT_GSL (gsl_linalg_SV_solve (u, v, s, b, x));
 101    push_vector (p, x);
 102    gsl_matrix_free (u);
 103    gsl_matrix_free (v);
 104    gsl_vector_free (b);
 105    gsl_vector_free (s);
 106    gsl_vector_free (x);
 107    (void) gsl_set_error_handler (save_handler);
 108  }
 109  
 110  #endif