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