1 /* linalg/test_qr.c
2 *
3 * Copyright (C) 2019-2020 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <stdlib.h>
22 #include <gsl/gsl_test.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_ieee_utils.h>
25 #include <gsl/gsl_permute_vector.h>
26 #include <gsl/gsl_blas.h>
27 #include <gsl/gsl_linalg.h>
28 #include <gsl/gsl_rng.h>
29 #include <gsl/gsl_permutation.h>
30
31 static int
test_QR_decomp_dim(const gsl_matrix * m,double eps)32 test_QR_decomp_dim(const gsl_matrix * m, double eps)
33 {
34 int s = 0;
35 const size_t M = m->size1;
36 const size_t N = m->size2;
37 size_t i, j;
38
39 gsl_matrix * qr = gsl_matrix_alloc(M, N);
40 gsl_matrix * a = gsl_matrix_alloc(M, N);
41 gsl_matrix * q = gsl_matrix_alloc(M, M);
42 gsl_matrix * r = gsl_matrix_alloc(M, N);
43 gsl_vector * d = gsl_vector_alloc(N);
44
45 gsl_matrix_memcpy(qr, m);
46
47 s += gsl_linalg_QR_decomp(qr, d);
48 s += gsl_linalg_QR_unpack(qr, d, q, r);
49
50 /* compute a = q r */
51 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, q, r, 0.0, a);
52
53 for(i=0; i<M; i++) {
54 for(j=0; j<N; j++) {
55 double aij = gsl_matrix_get(a, i, j);
56 double mij = gsl_matrix_get(m, i, j);
57 int foo = check(aij, mij, eps);
58 if(foo) {
59 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
60 }
61 s += foo;
62 }
63 }
64
65 gsl_vector_free(d);
66 gsl_matrix_free(qr);
67 gsl_matrix_free(a);
68 gsl_matrix_free(q);
69 gsl_matrix_free(r);
70
71 return s;
72 }
73
74 static int
test_QR_decomp(void)75 test_QR_decomp(void)
76 {
77 int f;
78 int s = 0;
79
80 f = test_QR_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
81 gsl_test(f, " QR_decomp m(3,5)");
82 s += f;
83
84 f = test_QR_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
85 gsl_test(f, " QR_decomp m(5,3)");
86 s += f;
87
88 f = test_QR_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
89 gsl_test(f, " QR_decomp hilbert(2)");
90 s += f;
91
92 f = test_QR_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
93 gsl_test(f, " QR_decomp hilbert(3)");
94 s += f;
95
96 f = test_QR_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
97 gsl_test(f, " QR_decomp hilbert(4)");
98 s += f;
99
100 f = test_QR_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
101 gsl_test(f, " QR_decomp hilbert(12)");
102 s += f;
103
104 f = test_QR_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
105 gsl_test(f, " QR_decomp vander(2)");
106 s += f;
107
108 f = test_QR_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
109 gsl_test(f, " QR_decomp vander(3)");
110 s += f;
111
112 f = test_QR_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
113 gsl_test(f, " QR_decomp vander(4)");
114 s += f;
115
116 f = test_QR_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
117 gsl_test(f, " QR_decomp vander(12)");
118 s += f;
119
120 return s;
121 }
122
123 static int
test_QR_decomp_r_eps(const gsl_matrix * m,const double eps,const char * desc)124 test_QR_decomp_r_eps(const gsl_matrix * m, const double eps, const char * desc)
125 {
126 int s = 0;
127 const size_t M = m->size1;
128 const size_t N = m->size2;
129 size_t i, j;
130
131 gsl_matrix * QR = gsl_matrix_alloc(M, N);
132 gsl_matrix * T = gsl_matrix_alloc(N, N);
133 gsl_matrix * A = gsl_matrix_alloc(M, N);
134 gsl_matrix * R = gsl_matrix_alloc(N, N);
135 gsl_matrix * Q = gsl_matrix_alloc(M, M);
136 gsl_matrix_view Q1 = gsl_matrix_submatrix(Q, 0, 0, M, N);
137
138 gsl_matrix_memcpy(QR, m);
139
140 s += gsl_linalg_QR_decomp_r(QR, T);
141 s += gsl_linalg_QR_unpack_r(QR, T, Q, R);
142
143 /* compute A = Q R */
144 gsl_matrix_memcpy(A, &Q1.matrix);
145 gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, A);
146
147 for (i = 0; i < M; i++)
148 {
149 for (j = 0; j < N; j++)
150 {
151 double aij = gsl_matrix_get(A, i, j);
152 double mij = gsl_matrix_get(m, i, j);
153
154 gsl_test_rel(aij, mij, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
155 desc, M, N, i,j, aij, mij);
156 }
157 }
158
159 if (M > N)
160 {
161 gsl_matrix * R_alt = gsl_matrix_alloc(M, N);
162 gsl_matrix * Q_alt = gsl_matrix_alloc(M, M);
163 gsl_vector_view tau = gsl_matrix_diagonal(T);
164
165 /* test that Q2 was computed correctly by comparing with Level 2 algorithm */
166 gsl_linalg_QR_unpack(QR, &tau.vector, Q_alt, R_alt);
167
168 for (i = 0; i < M; i++)
169 {
170 for (j = 0; j < M; j++)
171 {
172 double aij = gsl_matrix_get(Q, i, j);
173 double bij = gsl_matrix_get(Q_alt, i, j);
174
175 gsl_test_rel(aij, bij, eps, "%s Q (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
176 desc, M, N, i, j, aij, bij);
177 }
178 }
179
180 gsl_matrix_free(R_alt);
181 gsl_matrix_free(Q_alt);
182 }
183
184 gsl_matrix_free(QR);
185 gsl_matrix_free(T);
186 gsl_matrix_free(A);
187 gsl_matrix_free(Q);
188 gsl_matrix_free(R);
189
190 return s;
191 }
192
193 static int
test_QR_decomp_r(gsl_rng * r)194 test_QR_decomp_r(gsl_rng * r)
195 {
196 int s = 0;
197 size_t M, N;
198
199 for (M = 1; M <= 50; ++M)
200 {
201 for (N = 1; N <= M; ++N)
202 {
203 gsl_matrix * A = gsl_matrix_alloc(M, N);
204
205 create_random_matrix(A, r);
206 s += test_QR_decomp_r_eps(A, 1.0e6 * M * GSL_DBL_EPSILON, "QR_decomp_r random");
207
208 gsl_matrix_free(A);
209 }
210 }
211
212 s += test_QR_decomp_r_eps(m53, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r m(5,3)");
213 s += test_QR_decomp_r_eps(hilb2, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(2)");
214 s += test_QR_decomp_r_eps(hilb3, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(3)");
215 s += test_QR_decomp_r_eps(hilb4, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(4)");
216 s += test_QR_decomp_r_eps(hilb12, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(12)");
217 s += test_QR_decomp_r_eps(vander2, 1.0e1 * GSL_DBL_EPSILON, "QR_decomp_r vander(2)");
218 s += test_QR_decomp_r_eps(vander3, 1.0e1 * GSL_DBL_EPSILON, "QR_decomp_r vander(3)");
219 s += test_QR_decomp_r_eps(vander4, 1.0e1 * GSL_DBL_EPSILON, "QR_decomp_r vander(4)");
220
221 return s;
222 }
223
224 static int
test_QR_QTmat_r_eps(const gsl_matrix * A,const gsl_matrix * B,const double eps,const char * desc)225 test_QR_QTmat_r_eps(const gsl_matrix * A, const gsl_matrix * B, const double eps, const char * desc)
226 {
227 int s = 0;
228 const size_t M = A->size1;
229 const size_t N = A->size2;
230 const size_t K = B->size2;
231 size_t i, j;
232
233 gsl_matrix * QR = gsl_matrix_alloc(M, N);
234 gsl_matrix * T = gsl_matrix_alloc(N, N);
235 gsl_matrix * work = gsl_matrix_alloc(N, K);
236 gsl_matrix * B1 = gsl_matrix_alloc(M, K);
237 gsl_matrix * B2 = gsl_matrix_alloc(M, K);
238 gsl_vector_view tau = gsl_matrix_diagonal(T);
239
240 gsl_matrix_memcpy(QR, A);
241 gsl_matrix_memcpy(B1, B);
242 gsl_matrix_memcpy(B2, B);
243
244 s += gsl_linalg_QR_decomp_r(QR, T);
245
246 /* compute Q^T B with both recursive and non-recursive methods and compare */
247 s += gsl_linalg_QR_QTmat_r(QR, T, B1, work);
248 s += gsl_linalg_QR_QTmat(QR, &tau.vector, B2);
249
250 for (i = 0; i < M; i++)
251 {
252 for (j = 0; j < K; j++)
253 {
254 double aij = gsl_matrix_get(B1, i, j);
255 double bij = gsl_matrix_get(B2, i, j);
256
257 gsl_test_rel(aij, bij, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
258 desc, M, K, i,j, aij, bij);
259 }
260 }
261
262 gsl_matrix_free(QR);
263 gsl_matrix_free(T);
264 gsl_matrix_free(B1);
265 gsl_matrix_free(B2);
266 gsl_matrix_free(work);
267
268 return s;
269 }
270
271 static int
test_QR_QTmat_r(gsl_rng * r)272 test_QR_QTmat_r(gsl_rng * r)
273 {
274 int s = 0;
275 size_t M, N;
276
277 for (M = 1; M <= 50; ++M)
278 {
279 for (N = 1; N <= M; ++N)
280 {
281 const size_t K = GSL_MAX(N / 2, 1);
282 gsl_matrix * A = gsl_matrix_alloc(M, N);
283 gsl_matrix * B = gsl_matrix_alloc(M, K);
284
285 create_random_matrix(A, r);
286 create_random_matrix(B, r);
287 s += test_QR_QTmat_r_eps(A, B, 1.0e6 * M * GSL_DBL_EPSILON, "QR_QTmat_r random");
288
289 gsl_matrix_free(A);
290 gsl_matrix_free(B);
291 }
292 }
293
294 return s;
295 }
296
297 static int
test_QR_solve_r_eps(const gsl_matrix * A,const gsl_vector * rhs,const gsl_vector * sol,const double eps,const char * desc)298 test_QR_solve_r_eps(const gsl_matrix * A, const gsl_vector * rhs, const gsl_vector * sol,
299 const double eps, const char * desc)
300 {
301 int s = 0;
302 const size_t M = A->size1;
303 const size_t N = A->size2;
304 size_t i;
305
306 gsl_matrix * QR = gsl_matrix_alloc(M, N);
307 gsl_matrix * T = gsl_matrix_alloc(N, N);
308 gsl_vector * x = gsl_vector_alloc(N);
309
310 gsl_matrix_memcpy(QR, A);
311
312 s += gsl_linalg_QR_decomp_r(QR, T);
313 s += gsl_linalg_QR_solve_r(QR, T, rhs, x);
314
315 for (i = 0; i < M; i++)
316 {
317 double xi = gsl_vector_get(x, i);
318 double yi = gsl_vector_get(sol, i);
319
320 gsl_test_rel(xi, yi, eps, "%s (%3lu)[%lu]: %22.18g %22.18g\n",
321 desc, N, i, xi, yi);
322 }
323
324 gsl_matrix_free(QR);
325 gsl_matrix_free(T);
326 gsl_vector_free(x);
327
328 return s;
329 }
330
331 static int
test_QR_solve_r(gsl_rng * r)332 test_QR_solve_r(gsl_rng * r)
333 {
334 int s = 0;
335 size_t N;
336
337 for (N = 1; N <= 50; ++N)
338 {
339 gsl_matrix * A = gsl_matrix_alloc(N, N);
340 gsl_vector * sol = gsl_vector_alloc(N);
341 gsl_vector * rhs = gsl_vector_alloc(N);
342
343 create_random_matrix(A, r);
344 create_random_vector(sol, r);
345 gsl_blas_dgemv(CblasNoTrans, 1.0, A, sol, 0.0, rhs);
346
347 s += test_QR_solve_r_eps(A, rhs, sol, 1.0e5 * N * GSL_DBL_EPSILON, "QR_solve_r random");
348
349 gsl_matrix_free(A);
350 gsl_vector_free(sol);
351 gsl_vector_free(rhs);
352 }
353
354 return s;
355 }
356
357 static int
test_QR_lssolve_r_eps(const gsl_matrix * A,const gsl_vector * b,const double eps,const char * desc)358 test_QR_lssolve_r_eps(const gsl_matrix * A, const gsl_vector * b, const double eps, const char * desc)
359 {
360 int s = 0;
361 const size_t M = A->size1;
362 const size_t N = A->size2;
363 size_t i;
364
365 gsl_matrix * QR = gsl_matrix_alloc(M, N);
366 gsl_matrix * T = gsl_matrix_alloc(N, N);
367 gsl_vector * x = gsl_vector_alloc(M);
368 gsl_vector * work = gsl_vector_alloc(N);
369
370 gsl_matrix * U = gsl_matrix_alloc(M, N);
371 gsl_matrix * V = gsl_matrix_alloc(N, N);
372 gsl_vector * S = gsl_vector_alloc(N);
373 gsl_vector * x_svd = gsl_vector_alloc(N);
374
375 gsl_vector * residual = gsl_vector_alloc(M);
376
377 gsl_matrix_memcpy(QR, A);
378 s += gsl_linalg_QR_decomp_r(QR, T);
379 s += gsl_linalg_QR_lssolve_r(QR, T, b, x, work);
380
381 gsl_matrix_memcpy(U, A);
382 gsl_linalg_SV_decomp(U, V, S, work);
383 gsl_linalg_SV_solve(U, V, S, b, x_svd);
384
385 /* compare QR with SVD solution */
386 for (i = 0; i < N; i++)
387 {
388 double xi = gsl_vector_get(x, i);
389 double yi = gsl_vector_get(x_svd, i);
390
391 gsl_test_rel(xi, yi, eps, "%s (%3lu,%3lu)[%lu]: %22.18g %22.18g\n",
392 desc, M, N, i, xi, yi);
393 }
394
395 if (M > N)
396 {
397 gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);
398 gsl_vector_view x2 = gsl_vector_subvector(x, N, M - N);
399 double norm = gsl_blas_dnrm2(&x2.vector);
400 double norm_expected;
401
402 /* compute residual and check || x(N+1:end) || = || b - A x || */
403 gsl_vector_memcpy(residual, b);
404 gsl_blas_dgemv(CblasNoTrans, -1.0, A, &x1.vector, 1.0, residual);
405 norm_expected = gsl_blas_dnrm2(residual);
406
407 gsl_test_rel(norm, norm_expected, eps, "%s rnorm (%3lu,%3lu): %22.18g %22.18g\n",
408 desc, M, N, norm, norm_expected);
409 }
410
411 gsl_matrix_free(QR);
412 gsl_matrix_free(T);
413 gsl_vector_free(x);
414 gsl_matrix_free(U);
415 gsl_matrix_free(V);
416 gsl_vector_free(x_svd);
417 gsl_vector_free(work);
418 gsl_vector_free(S);
419 gsl_vector_free(residual);
420
421 return s;
422 }
423
424 static int
test_QR_lssolve_r(gsl_rng * r)425 test_QR_lssolve_r(gsl_rng * r)
426 {
427 int s = 0;
428 size_t M, N;
429
430 for (M = 1; M <= 30; ++M)
431 {
432 for (N = 1; N <= M; ++N)
433 {
434 gsl_matrix * A = gsl_matrix_alloc(M, N);
435 gsl_vector * b = gsl_vector_alloc(M);
436
437 create_random_matrix(A, r);
438 create_random_vector(b, r);
439
440 s += test_QR_lssolve_r_eps(A, b, 1.0e5 * M * GSL_DBL_EPSILON, "QR_lssolve_r random");
441
442 gsl_matrix_free(A);
443 gsl_vector_free(b);
444 }
445 }
446
447 return s;
448 }
449
450 static int
test_QR_UR_decomp_eps(const gsl_matrix * S,const gsl_matrix * A,const double eps,const char * desc)451 test_QR_UR_decomp_eps(const gsl_matrix * S, const gsl_matrix * A, const double eps, const char * desc)
452 {
453 int s = 0;
454 const size_t M = A->size1;
455 const size_t N = A->size2;
456 size_t i, j;
457
458 gsl_matrix * R = gsl_matrix_calloc(N, N);
459 gsl_matrix * V = gsl_matrix_alloc(M, N);
460 gsl_matrix * T = gsl_matrix_alloc(N, N);
461 gsl_matrix * B1 = gsl_matrix_calloc(N, N);
462 gsl_matrix * B2 = gsl_matrix_alloc(M, N);
463
464 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, S);
465 gsl_matrix_memcpy(V, A);
466
467 s += gsl_linalg_QR_UR_decomp(R, V, T);
468
469 /*
470 * compute B = Q R = [ R - T R ]
471 * [ -V~ T R ]
472 */
473
474 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);
475 gsl_matrix_memcpy(B2, V);
476
477 gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1); /* B1 = T R */
478 gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, B1, B2); /* B2 = -V~ T R */
479 gsl_matrix_sub (B1, R); /* B1 = T R - R */
480 gsl_matrix_scale (B1, -1.0); /* B1 = R - T R */
481
482 /* test S = B1 */
483 for (i = 0; i < N; i++)
484 {
485 for (j = i; j < N; j++)
486 {
487 double aij = gsl_matrix_get(B1, i, j);
488 double mij = gsl_matrix_get(S, i, j);
489
490 gsl_test_rel(aij, mij, eps, "%s B1 (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
491 desc, M, N, i,j, aij, mij);
492 }
493 }
494
495 /* test A = B2 */
496 for (i = 0; i < M; i++)
497 {
498 for (j = 0; j < N; j++)
499 {
500 double aij = gsl_matrix_get(B2, i, j);
501 double mij = gsl_matrix_get(A, i, j);
502
503 gsl_test_rel(aij, mij, eps, "%s B2 (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
504 desc, M, N, i,j, aij, mij);
505 }
506 }
507
508 gsl_matrix_free(R);
509 gsl_matrix_free(V);
510 gsl_matrix_free(T);
511 gsl_matrix_free(B1);
512 gsl_matrix_free(B2);
513
514 return s;
515 }
516
517 static int
test_QR_UR_decomp(gsl_rng * r)518 test_QR_UR_decomp(gsl_rng * r)
519 {
520 int s = 0;
521 size_t M, N;
522
523 for (M = 1; M <= 50; ++M)
524 {
525 for (N = 1; N <= M; ++N)
526 {
527 gsl_matrix * S = gsl_matrix_alloc(N, N);
528 gsl_matrix * A = gsl_matrix_alloc(M, N);
529
530 create_random_matrix(A, r);
531 create_random_matrix(S, r);
532 s += test_QR_UR_decomp_eps(S, A, 1.0e6 * M * GSL_DBL_EPSILON, "QR_UR_decomp random");
533
534 gsl_matrix_free(S);
535 gsl_matrix_free(A);
536 }
537 }
538
539 return s;
540 }
541
542 static int
test_QR_UZ_decomp_eps(const gsl_matrix * S,const gsl_matrix * A,const double eps,const char * desc)543 test_QR_UZ_decomp_eps(const gsl_matrix * S, const gsl_matrix * A, const double eps, const char * desc)
544 {
545 int s = 0;
546 const size_t M = A->size1;
547 const size_t N = A->size2;
548 size_t i, j;
549
550 gsl_matrix * R = gsl_matrix_calloc(N, N);
551 gsl_matrix * V = gsl_matrix_alloc(M, N);
552 gsl_matrix * T = gsl_matrix_alloc(N, N);
553 gsl_matrix * B1 = gsl_matrix_calloc(N, N);
554 gsl_matrix * B2 = gsl_matrix_alloc(M, N);
555
556 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, S);
557 gsl_matrix_memcpy(V, A);
558
559 s += gsl_linalg_QR_UZ_decomp(R, V, T);
560
561 /*
562 * compute B = Q R = [ R - T R ]
563 * [ -V~ T R ]
564 */
565
566 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);
567
568 gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1); /* B1 = T R */
569 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, -1.0, V, B1, 0.0, B2); /* B2 = -V~ T R */
570 gsl_matrix_sub (B1, R); /* B1 = T R - R */
571 gsl_matrix_scale (B1, -1.0); /* B1 = R - T R */
572
573 /* test S = B1 */
574 for (i = 0; i < N; i++)
575 {
576 for (j = i; j < N; j++)
577 {
578 double aij = gsl_matrix_get(B1, i, j);
579 double mij = gsl_matrix_get(S, i, j);
580
581 gsl_test_rel(aij, mij, eps, "%s B1 (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
582 desc, M, N, i,j, aij, mij);
583 }
584 }
585
586 /* test A = B2 */
587 for (i = 0; i < M; i++)
588 {
589 for (j = 0; j < N; j++)
590 {
591 double aij = gsl_matrix_get(B2, i, j);
592 double mij = gsl_matrix_get(A, i, j);
593
594 gsl_test_rel(aij, mij, eps, "%s B2 (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
595 desc, M, N, i,j, aij, mij);
596 }
597 }
598
599 gsl_matrix_free(R);
600 gsl_matrix_free(V);
601 gsl_matrix_free(T);
602 gsl_matrix_free(B1);
603 gsl_matrix_free(B2);
604
605 return s;
606 }
607
608 static int
test_QR_UZ_decomp(gsl_rng * r)609 test_QR_UZ_decomp(gsl_rng * r)
610 {
611 int s = 0;
612 const size_t N_max = 20;
613 const size_t M_max = 2*N_max;
614 gsl_matrix * U1 = gsl_matrix_alloc(N_max, N_max);
615 gsl_matrix * U2 = gsl_matrix_alloc(N_max, N_max);
616 gsl_matrix * A = gsl_matrix_alloc(M_max, N_max);
617 size_t M, N;
618
619 for (N = 1; N <= N_max; ++N)
620 {
621 gsl_matrix_view S = gsl_matrix_submatrix(U1, 0, 0, N, N);
622 gsl_matrix_view T = gsl_matrix_submatrix(U2, 0, 0, N, N);
623
624 for (M = N; M <= 2*N_max; ++M)
625 {
626 gsl_matrix_view B = gsl_matrix_submatrix(A, 0, 0, M, N);
627 gsl_matrix_view Bu = gsl_matrix_submatrix(&B.matrix, M - N, 0, N, N);
628
629 gsl_matrix_set_zero(&B.matrix);
630
631 if (M > N)
632 {
633 gsl_matrix_view Bd = gsl_matrix_submatrix(&B.matrix, 0, 0, M - N, N);
634 create_random_matrix(&Bd.matrix, r);
635 }
636
637 create_random_matrix(&S.matrix, r);
638 create_random_matrix(&T.matrix, r);
639 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &Bu.matrix, &T.matrix);
640
641 s += test_QR_UZ_decomp_eps(&S.matrix, &B.matrix, 1.0e6 * M * GSL_DBL_EPSILON,
642 "QR_UZ_decomp random");
643 }
644 }
645
646 gsl_matrix_free(U1);
647 gsl_matrix_free(U2);
648 gsl_matrix_free(A);
649
650 return s;
651 }
652
653 static int
test_QR_UU_decomp_eps(const gsl_matrix * U,const gsl_matrix * S,const double eps,const char * desc)654 test_QR_UU_decomp_eps(const gsl_matrix * U, const gsl_matrix * S, const double eps, const char * desc)
655 {
656 int s = 0;
657 const size_t N = U->size2;
658 size_t i, j;
659
660 gsl_matrix * R = gsl_matrix_calloc(N, N);
661 gsl_matrix * V = gsl_matrix_calloc(N, N);
662 gsl_matrix * T = gsl_matrix_alloc(N, N);
663 gsl_matrix * B1 = gsl_matrix_calloc(N, N);
664 gsl_matrix * B2 = gsl_matrix_calloc(N, N);
665
666 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);
667 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, V, S);
668
669 s += gsl_linalg_QR_UU_decomp(R, V, T);
670
671 /*
672 * compute B = Q R = [ R - T R ]
673 * [ -V~ T R ]
674 */
675
676 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);
677
678 gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1); /* B1 = T R */
679 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B2, B1); /* B2 := T R */
680 gsl_blas_dtrmm (CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, V, B2); /* B2 = -V~ T R */
681 gsl_matrix_sub (B1, R); /* B1 = T R - R */
682 gsl_matrix_scale (B1, -1.0); /* B1 = R - T R */
683
684 /* test U = B1 */
685 for (i = 0; i < N; i++)
686 {
687 for (j = i; j < N; j++)
688 {
689 double aij = gsl_matrix_get(B1, i, j);
690 double mij = gsl_matrix_get(U, i, j);
691
692 gsl_test_rel(aij, mij, eps, "%s B1 (%3lu)[%lu,%lu]: %22.18g %22.18g\n",
693 desc, N, i, j, aij, mij);
694 }
695 }
696
697 /* test S = B2 */
698 for (i = 0; i < N; i++)
699 {
700 for (j = 0; j < N; j++)
701 {
702 double aij = gsl_matrix_get(B2, i, j);
703 double mij = gsl_matrix_get(S, i, j);
704
705 gsl_test_rel(aij, mij, eps, "%s B2 (%3lu)[%lu,%lu]: %22.18g %22.18g\n",
706 desc, N, i,j, aij, mij);
707 }
708 }
709
710 gsl_matrix_free(R);
711 gsl_matrix_free(V);
712 gsl_matrix_free(T);
713 gsl_matrix_free(B1);
714 gsl_matrix_free(B2);
715
716 return s;
717 }
718
719 static int
test_QR_UU_decomp(gsl_rng * r)720 test_QR_UU_decomp(gsl_rng * r)
721 {
722 int s = 0;
723 const size_t N_max = 50;
724 gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
725 gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
726 gsl_matrix * T = gsl_matrix_alloc(N_max, N_max);
727 size_t N;
728
729 for (N = 1; N <= N_max; ++N)
730 {
731 gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
732 gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
733 gsl_matrix_view c = gsl_matrix_submatrix(T, 0, 0, N, N);
734
735 create_random_matrix(&c.matrix, r);
736 gsl_matrix_set_zero(&a.matrix);
737 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &c.matrix);
738
739 create_random_matrix(&c.matrix, r);
740 gsl_matrix_set_zero(&b.matrix);
741 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &b.matrix, &c.matrix);
742
743 s += test_QR_UU_decomp_eps(&a.matrix, &b.matrix, 1.0e6 * N * GSL_DBL_EPSILON,
744 "QR_UU_decomp random");
745 }
746
747 gsl_matrix_free(U);
748 gsl_matrix_free(S);
749 gsl_matrix_free(T);
750
751 return s;
752 }
753
754 static int
test_QR_UU_lssolve_eps(const gsl_matrix * U,const gsl_matrix * S,const gsl_vector * b,const double eps,const char * desc)755 test_QR_UU_lssolve_eps(const gsl_matrix * U, const gsl_matrix * S, const gsl_vector * b,
756 const double eps, const char * desc)
757 {
758 int s = 0;
759 const size_t N = U->size1;
760 const size_t M = 2 * N;
761 size_t i;
762
763 gsl_matrix * R = gsl_matrix_calloc(N, N);
764 gsl_matrix * Y = gsl_matrix_calloc(N, N);
765 gsl_matrix * T = gsl_matrix_alloc(N, N);
766 gsl_vector * x = gsl_vector_alloc(M);
767 gsl_vector * work = gsl_vector_alloc(N);
768
769 gsl_matrix * A = gsl_matrix_calloc(M, N);
770 gsl_matrix * U_svd = gsl_matrix_alloc(M, N);
771 gsl_matrix * V_svd = gsl_matrix_alloc(N, N);
772 gsl_vector * S_svd = gsl_vector_alloc(N);
773 gsl_vector * x_svd = gsl_vector_alloc(N);
774
775 gsl_vector * residual = gsl_vector_alloc(M);
776
777 gsl_matrix_view tmp;
778
779 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);
780 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, Y, S);
781
782 s += gsl_linalg_QR_UU_decomp(R, Y, T);
783 s += gsl_linalg_QR_UU_lssolve(R, Y, T, b, x, work);
784
785 tmp = gsl_matrix_submatrix(A, 0, 0, N, N);
786 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &tmp.matrix, U);
787
788 tmp = gsl_matrix_submatrix(A, N, 0, N, N);
789 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &tmp.matrix, S);
790
791 gsl_matrix_memcpy(U_svd, A);
792 gsl_linalg_SV_decomp(U_svd, V_svd, S_svd, work);
793 gsl_linalg_SV_solve(U_svd, V_svd, S_svd, b, x_svd);
794
795 /* compare QR with SVD solution */
796 for (i = 0; i < N; i++)
797 {
798 double xi = gsl_vector_get(x, i);
799 double yi = gsl_vector_get(x_svd, i);
800
801 gsl_test_rel(xi, yi, eps, "%s (%3lu,%3lu)[%lu]: %22.18g %22.18g\n",
802 desc, M, N, i, xi, yi);
803 }
804
805 if (M > N)
806 {
807 gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);
808 gsl_vector_view x2 = gsl_vector_subvector(x, N, M - N);
809 double norm = gsl_blas_dnrm2(&x2.vector);
810 double norm_expected;
811
812 /* compute residual and check || x(N+1:end) || = || b - A x || */
813 gsl_vector_memcpy(residual, b);
814 gsl_blas_dgemv(CblasNoTrans, -1.0, A, &x1.vector, 1.0, residual);
815 norm_expected = gsl_blas_dnrm2(residual);
816
817 gsl_test_rel(norm, norm_expected, eps, "%s rnorm (%3lu,%3lu): %22.18g %22.18g\n",
818 desc, M, N, norm, norm_expected);
819 }
820
821 gsl_matrix_free(A);
822 gsl_matrix_free(R);
823 gsl_matrix_free(Y);
824 gsl_matrix_free(T);
825 gsl_vector_free(x);
826 gsl_matrix_free(U_svd);
827 gsl_matrix_free(V_svd);
828 gsl_vector_free(x_svd);
829 gsl_vector_free(work);
830 gsl_vector_free(S_svd);
831 gsl_vector_free(residual);
832
833 return s;
834 }
835
836 static int
test_QR_UU_lssolve(gsl_rng * r)837 test_QR_UU_lssolve(gsl_rng * r)
838 {
839 int s = 0;
840 const size_t N_max = 30;
841 gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
842 gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
843 gsl_matrix * T = gsl_matrix_alloc(N_max, N_max);
844 gsl_vector * rhs = gsl_vector_alloc(2*N_max);
845 size_t N;
846
847 for (N = 1; N <= N_max; ++N)
848 {
849 gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
850 gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
851 gsl_matrix_view c = gsl_matrix_submatrix(T, 0, 0, N, N);
852 gsl_vector_view rhsv = gsl_vector_subvector(rhs, 0, 2*N);
853
854 create_random_vector(&rhsv.vector, r);
855
856 create_random_matrix(&c.matrix, r);
857 gsl_matrix_set_zero(&a.matrix);
858 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &c.matrix);
859
860 create_random_matrix(&c.matrix, r);
861 gsl_matrix_set_zero(&b.matrix);
862 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &b.matrix, &c.matrix);
863
864 s += test_QR_UU_lssolve_eps(&a.matrix, &b.matrix, &rhsv.vector, 1.0e4 * N * GSL_DBL_EPSILON,
865 "QR_UU_lssolve random");
866 }
867
868 gsl_matrix_free(U);
869 gsl_matrix_free(S);
870 gsl_matrix_free(T);
871 gsl_vector_free(rhs);
872
873 return s;
874 }
875
876 static int
test_QR_UD_decomp_eps(const gsl_matrix * U,const gsl_vector * D,const double eps,const char * desc)877 test_QR_UD_decomp_eps(const gsl_matrix * U, const gsl_vector * D, const double eps, const char * desc)
878 {
879 int s = 0;
880 const size_t N = U->size2;
881 size_t i, j;
882
883 gsl_matrix * R = gsl_matrix_calloc(N, N);
884 gsl_matrix * V = gsl_matrix_calloc(N, N);
885 gsl_matrix * T = gsl_matrix_alloc(N, N);
886 gsl_matrix * B1 = gsl_matrix_calloc(N, N);
887 gsl_matrix * B2 = gsl_matrix_calloc(N, N);
888
889 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);
890
891 s += gsl_linalg_QR_UD_decomp(R, D, V, T);
892
893 /*
894 * compute B = Q R = [ R - T R ]
895 * [ -V~ T R ]
896 */
897
898 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);
899
900 gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1); /* B1 = T R */
901 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B2, B1); /* B2 := T R */
902 gsl_blas_dtrmm (CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, V, B2); /* B2 = -V~ T R */
903 gsl_matrix_sub (B1, R); /* B1 = T R - R */
904 gsl_matrix_scale (B1, -1.0); /* B1 = R - T R */
905
906 /* test U = B1 */
907 for (i = 0; i < N; i++)
908 {
909 for (j = i; j < N; j++)
910 {
911 double aij = gsl_matrix_get(B1, i, j);
912 double mij = gsl_matrix_get(U, i, j);
913
914 gsl_test_rel(aij, mij, eps, "%s B1 (%3lu)[%lu,%lu]: %22.18g %22.18g\n",
915 desc, N, i, j, aij, mij);
916 }
917 }
918
919 /* test S = B2 */
920 for (i = 0; i < N; i++)
921 {
922 for (j = 0; j < N; j++)
923 {
924 double aij = gsl_matrix_get(B2, i, j);
925 double mij = (i == j) ? gsl_vector_get(D, i) : 0.0;
926
927 gsl_test_rel(aij, mij, eps, "%s B2 (%3lu)[%lu,%lu]: %22.18g %22.18g\n",
928 desc, N, i,j, aij, mij);
929 }
930 }
931
932 gsl_matrix_free(R);
933 gsl_matrix_free(V);
934 gsl_matrix_free(T);
935 gsl_matrix_free(B1);
936 gsl_matrix_free(B2);
937
938 return s;
939 }
940
941 static int
test_QR_UD_decomp(gsl_rng * r)942 test_QR_UD_decomp(gsl_rng * r)
943 {
944 int s = 0;
945 const size_t N_max = 30;
946 gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
947 gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
948 gsl_vector * D = gsl_vector_alloc(N_max);
949 size_t N;
950
951 for (N = 1; N <= N_max; ++N)
952 {
953 gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
954 gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
955 gsl_vector_view diag = gsl_vector_subvector(D, 0, N);
956
957 create_random_matrix(&b.matrix, r);
958 gsl_matrix_set_zero(&a.matrix);
959 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &b.matrix);
960
961 create_random_vector(&diag.vector, r);
962
963 s += test_QR_UD_decomp_eps(&a.matrix, &diag.vector, 1.0e4 * N * GSL_DBL_EPSILON,
964 "QR_UD_decomp random");
965 }
966
967 gsl_matrix_free(U);
968 gsl_matrix_free(S);
969 gsl_vector_free(D);
970
971 return s;
972 }
973
974 static int
test_QR_UD_lssolve_eps(const gsl_matrix * U,const gsl_vector * D,const gsl_vector * b,const double eps,const char * desc)975 test_QR_UD_lssolve_eps(const gsl_matrix * U, const gsl_vector * D, const gsl_vector * b,
976 const double eps, const char * desc)
977 {
978 int s = 0;
979 const size_t N = U->size1;
980 const size_t M = 2 * N;
981 size_t i;
982
983 gsl_matrix * R = gsl_matrix_calloc(N, N);
984 gsl_matrix * Y = gsl_matrix_calloc(N, N);
985 gsl_matrix * T = gsl_matrix_alloc(N, N);
986 gsl_vector * x = gsl_vector_alloc(M);
987 gsl_vector * work = gsl_vector_alloc(N);
988
989 gsl_matrix * A = gsl_matrix_calloc(M, N);
990 gsl_matrix * U_svd = gsl_matrix_alloc(M, N);
991 gsl_matrix * V_svd = gsl_matrix_alloc(N, N);
992 gsl_vector * S_svd = gsl_vector_alloc(N);
993 gsl_vector * x_svd = gsl_vector_alloc(N);
994
995 gsl_vector * residual = gsl_vector_alloc(M);
996
997 gsl_matrix_view tmp;
998 gsl_vector_view diag;
999
1000 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);
1001
1002 s += gsl_linalg_QR_UD_decomp(R, D, Y, T);
1003 s += gsl_linalg_QR_UD_lssolve(R, Y, T, b, x, work);
1004
1005 tmp = gsl_matrix_submatrix(A, 0, 0, N, N);
1006 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &tmp.matrix, U);
1007
1008 tmp = gsl_matrix_submatrix(A, N, 0, N, N);
1009 diag = gsl_matrix_diagonal(&tmp.matrix);
1010 gsl_vector_memcpy(&diag.vector, D);
1011
1012 gsl_matrix_memcpy(U_svd, A);
1013 gsl_linalg_SV_decomp(U_svd, V_svd, S_svd, work);
1014 gsl_linalg_SV_solve(U_svd, V_svd, S_svd, b, x_svd);
1015
1016 /* compare QR with SVD solution */
1017 for (i = 0; i < N; i++)
1018 {
1019 double xi = gsl_vector_get(x, i);
1020 double yi = gsl_vector_get(x_svd, i);
1021
1022 gsl_test_rel(xi, yi, eps, "%s (%3lu,%3lu)[%lu]: %22.18g %22.18g\n",
1023 desc, M, N, i, xi, yi);
1024 }
1025
1026 if (M > N)
1027 {
1028 gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);
1029 gsl_vector_view x2 = gsl_vector_subvector(x, N, M - N);
1030 double norm = gsl_blas_dnrm2(&x2.vector);
1031 double norm_expected;
1032
1033 /* compute residual and check || x(N+1:end) || = || b - A x || */
1034 gsl_vector_memcpy(residual, b);
1035 gsl_blas_dgemv(CblasNoTrans, -1.0, A, &x1.vector, 1.0, residual);
1036 norm_expected = gsl_blas_dnrm2(residual);
1037
1038 gsl_test_rel(norm, norm_expected, eps, "%s rnorm (%3lu,%3lu): %22.18g %22.18g\n",
1039 desc, M, N, norm, norm_expected);
1040 }
1041
1042 gsl_matrix_free(A);
1043 gsl_matrix_free(R);
1044 gsl_matrix_free(Y);
1045 gsl_matrix_free(T);
1046 gsl_vector_free(x);
1047 gsl_matrix_free(U_svd);
1048 gsl_matrix_free(V_svd);
1049 gsl_vector_free(x_svd);
1050 gsl_vector_free(work);
1051 gsl_vector_free(S_svd);
1052 gsl_vector_free(residual);
1053
1054 return s;
1055 }
1056
1057 static int
test_QR_UD_lssolve(gsl_rng * r)1058 test_QR_UD_lssolve(gsl_rng * r)
1059 {
1060 int s = 0;
1061 const size_t N_max = 30;
1062 gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
1063 gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
1064 gsl_vector * D = gsl_vector_alloc(N_max);
1065 gsl_vector * rhs = gsl_vector_alloc(2*N_max);
1066 size_t N;
1067
1068 for (N = 1; N <= N_max; ++N)
1069 {
1070 gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
1071 gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
1072 gsl_vector_view rhsv = gsl_vector_subvector(rhs, 0, 2*N);
1073 gsl_vector_view diag = gsl_vector_subvector(D, 0, N);
1074
1075 create_random_vector(&rhsv.vector, r);
1076 create_random_vector(&diag.vector, r);
1077
1078 create_random_matrix(&b.matrix, r);
1079 gsl_matrix_set_zero(&a.matrix);
1080 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &b.matrix);
1081
1082 s += test_QR_UD_lssolve_eps(&a.matrix, &diag.vector, &rhsv.vector, 1.0e4 * N * GSL_DBL_EPSILON,
1083 "QR_UD_lssolve random");
1084 }
1085
1086 gsl_matrix_free(U);
1087 gsl_matrix_free(S);
1088 gsl_vector_free(D);
1089 gsl_vector_free(rhs);
1090
1091 return s;
1092 }
1093