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-2024 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 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 A68_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, A68_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 (A68_REF, &ref_u) = matrix_to_row (p, u);
53 }
54 if (!IS_NIL (ref_s)) {
55 *DEREF (A68_REF, &ref_s) = vector_to_row (p, s);
56 }
57 if (!IS_NIL (ref_v)) {
58 *DEREF (A68_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 (A68_REF, &ref_u) = matrix_to_row (p, v);
70 }
71 if (!IS_NIL (ref_s)) {
72 *DEREF (A68_REF, &ref_s) = vector_to_row (p, s);
73 }
74 if (!IS_NIL (ref_v)) {
75 *DEREF (A68_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, A68_TRUE);
92 gsl_matrix *v = pop_matrix (p, A68_TRUE);
93 gsl_vector *s = pop_vector (p, A68_TRUE);
94 gsl_matrix *u = pop_matrix (p, A68_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
