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
