1 /* linalg/test.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2005, 2006, 2007, 2010 Gerard Jungman, Brian Gough
4 * Copyright (C) Huan Wu (test_choleskyc_invert and test_choleskyc_invert_dim)
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /* Author: G. Jungman
22 */
23 #include <config.h>
24 #include <stdlib.h>
25 #include <gsl/gsl_test.h>
26 #include <gsl/gsl_math.h>
27 #include <gsl/gsl_ieee_utils.h>
28 #include <gsl/gsl_permute_vector.h>
29 #include <gsl/gsl_permute_matrix.h>
30 #include <gsl/gsl_blas.h>
31 #include <gsl/gsl_complex_math.h>
32 #include <gsl/gsl_linalg.h>
33 #include <gsl/gsl_rng.h>
34
35 #define TEST_SVD_4X4 1
36
37 int check (double x, double actual, double eps);
38 gsl_matrix * create_hilbert_matrix(unsigned long size);
39 gsl_matrix * create_general_matrix(unsigned long size1, unsigned long size2);
40 gsl_matrix * create_vandermonde_matrix(unsigned long size);
41 gsl_matrix * create_moler_matrix(unsigned long size);
42 gsl_matrix * create_row_matrix(unsigned long size1, unsigned long size2);
43 gsl_matrix * create_2x2_matrix(double a11, double a12, double a21, double a22);
44 gsl_matrix * create_diagonal_matrix(double a[], unsigned long size);
45 gsl_matrix * create_sparse_matrix(unsigned long m, unsigned long n);
46
47 int test_matmult(void);
48 int test_matmult_mod(void);
49 int test_QR_solve_dim(const gsl_matrix * m, const double * actual, double eps);
50 int test_QR_solve(void);
51 int test_QR_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps);
52 int test_QR_QRsolve(void);
53 int test_QR_lssolve_dim(const gsl_matrix * m, const double * actual, double eps);
54 int test_QR_lssolve(void);
55 int test_QRPT_solve_dim(const gsl_matrix * m, const double * actual, double eps);
56 int test_QRPT_solve(void);
57 int test_QRPT_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps);
58 int test_QRPT_QRsolve(void);
59 int test_QRPT_decomp_dim(const gsl_matrix * m, const double expected_rcond, double eps);
60 int test_QRPT_decomp(void);
61 int test_QRPT_lssolve_dim(const gsl_matrix * m, const double * actual, double eps);
62 int test_QRPT_lssolve(void);
63 int test_QRPT_lssolve2_dim(const gsl_matrix * m, const double * actual, double eps);
64 int test_QRPT_lssolve2(void);
65 int test_QR_update_dim(const gsl_matrix * m, double eps);
66 int test_QR_update(void);
67 int test_QRPT_update_dim(const gsl_matrix * m, double eps);
68 int test_QRPT_update(void);
69
70 int test_SV_solve_dim(const gsl_matrix * m, const double * actual, double eps);
71 int test_SV_solve(void);
72 int test_SV_decomp_dim(const gsl_matrix * m, double eps);
73 int test_SV_decomp(void);
74 int test_SV_decomp_mod_dim(const gsl_matrix * m, double eps);
75 int test_SV_decomp_mod(void);
76 int test_SV_decomp_jacobi_dim(const gsl_matrix * m, double eps);
77 int test_SV_decomp_jacobi(void);
78 int test_cholesky_solve_dim(const gsl_matrix * m, const double * actual, double eps);
79 int test_cholesky_solve(void);
80 int test_HH_solve_dim(const gsl_matrix * m, const double * actual, double eps);
81 int test_HH_solve(void);
82 int test_TDS_solve_dim(unsigned long dim, double d, double od, const double * actual, double eps);
83 int test_TDS_solve(void);
84 int test_TDN_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps);
85 int test_TDN_solve(void);
86 int test_TDS_cyc_solve_one(const unsigned long dim, const double * d, const double * od, const double * r,
87 const double * actual, double eps);
88 int test_TDS_cyc_solve(void);
89 int test_TDN_cyc_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps);
90 int test_TDN_cyc_solve(void);
91 int test_bidiag_decomp_dim(const gsl_matrix * m, double eps);
92 int test_bidiag_decomp(void);
93
94 int
check(double x,double actual,double eps)95 check (double x, double actual, double eps)
96 {
97 if (x == actual)
98 {
99 return 0;
100 }
101 else if (actual == 0)
102 {
103 return fabs(x) > eps;
104 }
105 else
106 {
107 return (fabs(x - actual)/fabs(actual)) > eps;
108 }
109 }
110
111
112 gsl_vector *
vector_alloc(size_t n)113 vector_alloc (size_t n)
114 {
115 size_t p[5] = {3, 5, 7, 11, 13};
116 static size_t k = 0;
117
118 size_t stride = p[k];
119 k = (k + 1) % 5;
120
121 {
122 gsl_block * b = gsl_block_alloc (n * stride);
123 gsl_vector * v = gsl_vector_alloc_from_block (b, 0, n, stride);
124 v->owner = 1;
125 return v;
126 }
127 }
128
129 void
vector_free(gsl_vector * v)130 vector_free (gsl_vector * v)
131 {
132 gsl_vector_free (v);
133 }
134
135 gsl_matrix *
create_hilbert_matrix(unsigned long size)136 create_hilbert_matrix(unsigned long size)
137 {
138 unsigned long i, j;
139 gsl_matrix * m = gsl_matrix_alloc(size, size);
140 for(i=0; i<size; i++) {
141 for(j=0; j<size; j++) {
142 gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
143 }
144 }
145 return m;
146 }
147
148 gsl_matrix *
create_general_matrix(unsigned long size1,unsigned long size2)149 create_general_matrix(unsigned long size1, unsigned long size2)
150 {
151 unsigned long i, j;
152 gsl_matrix * m = gsl_matrix_alloc(size1, size2);
153 for(i=0; i<size1; i++) {
154 for(j=0; j<size2; j++) {
155 gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
156 }
157 }
158 return m;
159 }
160
161 gsl_matrix *
create_singular_matrix(unsigned long size1,unsigned long size2)162 create_singular_matrix(unsigned long size1, unsigned long size2)
163 {
164 unsigned long i, j;
165 gsl_matrix * m = gsl_matrix_alloc(size1, size2);
166 for(i=0; i<size1; i++) {
167 for(j=0; j<size2; j++) {
168 gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
169 }
170 }
171
172 /* zero the first column */
173 for(j = 0; j <size2; j++)
174 gsl_matrix_set(m,0,j,0.0);
175
176 return m;
177 }
178
179
180 gsl_matrix *
create_vandermonde_matrix(unsigned long size)181 create_vandermonde_matrix(unsigned long size)
182 {
183 unsigned long i, j;
184 gsl_matrix * m = gsl_matrix_alloc(size, size);
185 for(i=0; i<size; i++) {
186 for(j=0; j<size; j++) {
187 gsl_matrix_set(m, i, j, pow(i + 1.0, size - j - 1.0));
188 }
189 }
190 return m;
191 }
192
193 gsl_matrix *
create_moler_matrix(unsigned long size)194 create_moler_matrix(unsigned long size)
195 {
196 unsigned long i, j;
197 gsl_matrix * m = gsl_matrix_alloc(size, size);
198 for(i=0; i<size; i++) {
199 for(j=0; j<size; j++) {
200 gsl_matrix_set(m, i, j, GSL_MIN(i+1,j+1)-2.0);
201 }
202 }
203 return m;
204 }
205
206 gsl_matrix_complex *
create_complex_matrix(unsigned long size)207 create_complex_matrix(unsigned long size)
208 {
209 unsigned long i, j;
210 gsl_matrix_complex * m = gsl_matrix_complex_alloc(size, size);
211 for(i=0; i<size; i++) {
212 for(j=0; j<size; j++) {
213 gsl_complex z = gsl_complex_rect(1.0/(i+j+1.0), 1/(i*i+j*j+0.5));
214 gsl_matrix_complex_set(m, i, j, z);
215 }
216 }
217 return m;
218 }
219
220 gsl_matrix *
create_row_matrix(unsigned long size1,unsigned long size2)221 create_row_matrix(unsigned long size1, unsigned long size2)
222 {
223 unsigned long i;
224 gsl_matrix * m = gsl_matrix_calloc(size1, size2);
225 for(i=0; i<size1; i++) {
226 gsl_matrix_set(m, i, 0, 1.0/(i+1.0));
227 }
228
229 return m;
230 }
231
232 gsl_matrix *
create_2x2_matrix(double a11,double a12,double a21,double a22)233 create_2x2_matrix(double a11, double a12, double a21, double a22)
234 {
235 gsl_matrix * m = gsl_matrix_alloc(2, 2);
236 gsl_matrix_set(m, 0, 0, a11);
237 gsl_matrix_set(m, 0, 1, a12);
238 gsl_matrix_set(m, 1, 0, a21);
239 gsl_matrix_set(m, 1, 1, a22);
240 return m;
241 }
242
243 gsl_matrix *
create_diagonal_matrix(double a[],unsigned long size)244 create_diagonal_matrix(double a[], unsigned long size)
245 {
246 unsigned long i;
247 gsl_matrix * m = gsl_matrix_calloc(size, size);
248 for(i=0; i<size; i++) {
249 gsl_matrix_set(m, i, i, a[i]);
250 }
251
252 return m;
253 }
254
rand_double()255 double rand_double() {
256 static unsigned int x;
257 x = (69069 * x + 1) & 0xFFFFFFFFUL;
258 return (x / 4294967296.0);
259 }
260
261 gsl_matrix *
create_sparse_matrix(unsigned long m,unsigned long n)262 create_sparse_matrix(unsigned long m, unsigned long n) {
263 gsl_matrix* A = gsl_matrix_calloc(m, n);
264
265 unsigned long int i, j;
266
267 for (i = 0; i < m; i++) {
268 for (j = 0; j < n; j++) {
269 double a = (rand_double() < 0.6 ? 1 : 0);
270 gsl_matrix_set(A, i, j, a);
271 }
272 }
273
274 for (i = 0; i < m; i++) {
275 if (rand_double() < 0.9) {
276 gsl_vector_view row = gsl_matrix_row (A, i);
277 gsl_vector_set_zero (&row.vector);
278 }
279 }
280 for (i = 0; i < n; i++) {
281 if (rand_double() < 0.43) {
282 gsl_vector_view col = gsl_matrix_column (A, i);
283 gsl_vector_set_zero (&col.vector);
284 }
285 }
286
287 return A;
288 }
289
290 static int
create_tri_matrix(CBLAS_UPLO_t Uplo,CBLAS_DIAG_t Diag,gsl_matrix * m,gsl_rng * r)291 create_tri_matrix(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * m, gsl_rng * r)
292 {
293 const size_t N = m->size1;
294 size_t i, j;
295
296 gsl_matrix_set_zero(m);
297
298 for (i = 0; i < N; ++i)
299 {
300 for (j = 0; j <= i; ++j)
301 {
302 double mij = gsl_rng_uniform(r);
303
304 /* ensure diagonally dominant matrix */
305 if (i == j)
306 {
307 if (Diag == CblasUnit)
308 mij = 1.0;
309 else
310 mij += 10.0;
311 }
312 else if (Diag == CblasUnit)
313 mij *= 0.01;
314
315 if (Uplo == CblasLower)
316 gsl_matrix_set(m, i, j, mij);
317 else
318 gsl_matrix_set(m, j, i, mij);
319 }
320 }
321
322 return GSL_SUCCESS;
323 }
324
325
326 gsl_matrix * m11;
327 gsl_matrix * m51;
328
329 gsl_matrix * m35;
330 gsl_matrix * m53;
331 gsl_matrix * m97;
332
333 gsl_matrix * s35;
334 gsl_matrix * s53;
335
336 gsl_matrix * hilb2;
337 gsl_matrix * hilb3;
338 gsl_matrix * hilb4;
339 gsl_matrix * hilb12;
340
341 gsl_matrix * row3;
342 gsl_matrix * row5;
343 gsl_matrix * row12;
344
345 gsl_matrix * A22;
346 gsl_matrix * A33;
347 gsl_matrix * A44;
348 gsl_matrix * A55;
349
350 gsl_matrix_complex * c7;
351
352 gsl_matrix * inf5; double inf5_data[] = {1.0, 0.0, -3.0, 0.0, -5.0};
353 gsl_matrix * nan5;
354 gsl_matrix * dblmin3, * dblmin5, * dblsubnorm5;
355 gsl_matrix * bigsparse;
356
357 double m53_lssolution[] = {52.5992295702070, -337.7263113752073,
358 351.8823436427604};
359 double hilb2_solution[] = {-8.0, 18.0} ;
360 double hilb3_solution[] = {27.0, -192.0, 210.0};
361 double hilb4_solution[] = {-64.0, 900.0, -2520.0, 1820.0};
362 double hilb12_solution[] = {-1728.0, 245388.0, -8528520.0,
363 127026900.0, -1009008000.0, 4768571808.0,
364 -14202796608.0, 27336497760.0, -33921201600.0,
365 26189163000.0, -11437874448.0, 2157916488.0 };
366
367 double c7_solution[] = { 2.40717272023734e+01, -9.84612797621247e+00,
368 -2.69338853034031e+02, 8.75455232472528e+01,
369 2.96661356736296e+03, -1.02624473923993e+03,
370 -1.82073812124749e+04, 5.67384473042410e+03,
371 5.57693879019068e+04, -1.61540963210502e+04,
372 -7.88941207561151e+04, 1.95053812987858e+04,
373 3.95548551241728e+04, -7.76593696255317e+03 };
374
375 gsl_matrix * vander2;
376 gsl_matrix * vander3;
377 gsl_matrix * vander4;
378 gsl_matrix * vander12;
379
380 double vander2_solution[] = {1.0, 0.0};
381 double vander3_solution[] = {0.0, 1.0, 0.0};
382 double vander4_solution[] = {0.0, 0.0, 1.0, 0.0};
383 double vander12_solution[] = {0.0, 0.0, 0.0, 0.0,
384 0.0, 0.0, 0.0, 0.0,
385 0.0, 0.0, 1.0, 0.0};
386
387 gsl_matrix * moler10;
388
389 #include "test_common.c"
390 #include "test_cholesky.c"
391 #include "test_choleskyc.c"
392 #include "test_cod.c"
393 #include "test_ldlt.c"
394 #include "test_lu.c"
395 #include "test_luc.c"
396 #include "test_lu_band.c"
397 #include "test_lq.c"
398 #include "test_tri.c"
399 #include "test_ql.c"
400 #include "test_qr.c"
401 #include "test_qrc.c"
402 #include "test_qr_band.c"
403
404 int
test_QR_solve_dim(const gsl_matrix * m,const double * actual,double eps)405 test_QR_solve_dim(const gsl_matrix * m, const double * actual, double eps)
406 {
407 int s = 0;
408 unsigned long i, dim = m->size1;
409
410 gsl_vector * rhs = gsl_vector_alloc(dim);
411 gsl_matrix * qr = gsl_matrix_alloc(dim,dim);
412 gsl_vector * d = gsl_vector_alloc(dim);
413 gsl_vector * x = gsl_vector_alloc(dim);
414
415 gsl_matrix_memcpy(qr,m);
416 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
417 s += gsl_linalg_QR_decomp(qr, d);
418 s += gsl_linalg_QR_solve(qr, d, rhs, x);
419 for(i=0; i<dim; i++) {
420 int foo = check(gsl_vector_get(x, i), actual[i], eps);
421 if(foo) {
422 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
423 }
424 s += foo;
425 }
426
427 gsl_vector_free(x);
428 gsl_vector_free(d);
429 gsl_matrix_free(qr);
430 gsl_vector_free(rhs);
431
432 return s;
433 }
434
test_QR_solve(void)435 int test_QR_solve(void)
436 {
437 int f;
438 int s = 0;
439
440 f = test_QR_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
441 gsl_test(f, " QR_solve hilbert(2)");
442 s += f;
443
444 f = test_QR_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
445 gsl_test(f, " QR_solve hilbert(3)");
446 s += f;
447
448 f = test_QR_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
449 gsl_test(f, " QR_solve hilbert(4)");
450 s += f;
451
452 f = test_QR_solve_dim(hilb12, hilb12_solution, 0.5);
453 gsl_test(f, " QR_solve hilbert(12)");
454 s += f;
455
456 f = test_QR_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
457 gsl_test(f, " QR_solve vander(2)");
458 s += f;
459
460 f = test_QR_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
461 gsl_test(f, " QR_solve vander(3)");
462 s += f;
463
464 f = test_QR_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
465 gsl_test(f, " QR_solve vander(4)");
466 s += f;
467
468 f = test_QR_solve_dim(vander12, vander12_solution, 0.05);
469 gsl_test(f, " QR_solve vander(12)");
470 s += f;
471
472 return s;
473 }
474
475
476 int
test_QR_QRsolve_dim(const gsl_matrix * m,const double * actual,double eps)477 test_QR_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps)
478 {
479 int s = 0;
480 unsigned long i, dim = m->size1;
481
482 gsl_vector * rhs = gsl_vector_alloc(dim);
483 gsl_matrix * qr = gsl_matrix_alloc(dim,dim);
484 gsl_matrix * q = gsl_matrix_alloc(dim,dim);
485 gsl_matrix * r = gsl_matrix_alloc(dim,dim);
486 gsl_vector * d = gsl_vector_alloc(dim);
487 gsl_vector * x = gsl_vector_alloc(dim);
488
489 gsl_matrix_memcpy(qr,m);
490 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
491 s += gsl_linalg_QR_decomp(qr, d);
492 s += gsl_linalg_QR_unpack(qr, d, q, r);
493 s += gsl_linalg_QR_QRsolve(q, r, rhs, x);
494 for(i=0; i<dim; i++) {
495 int foo = check(gsl_vector_get(x, i), actual[i], eps);
496 if(foo) {
497 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
498 }
499 s += foo;
500 }
501
502 gsl_vector_free(x);
503 gsl_vector_free(d);
504 gsl_matrix_free(qr);
505 gsl_matrix_free(q);
506 gsl_matrix_free(r);
507 gsl_vector_free(rhs);
508 return s;
509 }
510
test_QR_QRsolve(void)511 int test_QR_QRsolve(void)
512 {
513 int f;
514 int s = 0;
515
516 f = test_QR_QRsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
517 gsl_test(f, " QR_QRsolve hilbert(2)");
518 s += f;
519
520 f = test_QR_QRsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
521 gsl_test(f, " QR_QRsolve hilbert(3)");
522 s += f;
523
524 f = test_QR_QRsolve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
525 gsl_test(f, " QR_QRsolve hilbert(4)");
526 s += f;
527
528 f = test_QR_QRsolve_dim(hilb12, hilb12_solution, 0.5);
529 gsl_test(f, " QR_QRsolve hilbert(12)");
530 s += f;
531
532 f = test_QR_QRsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
533 gsl_test(f, " QR_QRsolve vander(2)");
534 s += f;
535
536 f = test_QR_QRsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
537 gsl_test(f, " QR_QRsolve vander(3)");
538 s += f;
539
540 f = test_QR_QRsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
541 gsl_test(f, " QR_QRsolve vander(4)");
542 s += f;
543
544 f = test_QR_QRsolve_dim(vander12, vander12_solution, 0.05);
545 gsl_test(f, " QR_QRsolve vander(12)");
546 s += f;
547
548 return s;
549 }
550
551
552 int
test_QR_lssolve_dim(const gsl_matrix * m,const double * actual,double eps)553 test_QR_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
554 {
555 int s = 0;
556 unsigned long i, M = m->size1, N = m->size2;
557
558 gsl_vector * rhs = gsl_vector_alloc(M);
559 gsl_matrix * qr = gsl_matrix_alloc(M,N);
560 gsl_vector * d = gsl_vector_alloc(N);
561 gsl_vector * x = gsl_vector_alloc(N);
562 gsl_vector * r = gsl_vector_alloc(M);
563 gsl_vector * res = gsl_vector_alloc(M);
564
565 gsl_matrix_memcpy(qr,m);
566 for(i=0; i<M; i++) gsl_vector_set(rhs, i, i+1.0);
567 s += gsl_linalg_QR_decomp(qr, d);
568 s += gsl_linalg_QR_lssolve(qr, d, rhs, x, res);
569
570 for(i=0; i<N; i++) {
571 int foo = check(gsl_vector_get(x, i), actual[i], eps);
572 if(foo) {
573 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
574 }
575 s += foo;
576 }
577
578 /* compute residual r = b - m x */
579 if (M == N) {
580 gsl_vector_set_zero(r);
581 } else {
582 gsl_vector_memcpy(r, rhs);
583 gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
584 };
585
586 for(i=0; i<N; i++) {
587 int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
588 if(foo) {
589 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
590 }
591 s += foo;
592 }
593
594 gsl_vector_free(r);
595 gsl_vector_free(res);
596 gsl_vector_free(x);
597 gsl_vector_free(d);
598 gsl_matrix_free(qr);
599 gsl_vector_free(rhs);
600
601 return s;
602 }
603
test_QR_lssolve(void)604 int test_QR_lssolve(void)
605 {
606 int f;
607 int s = 0;
608
609 f = test_QR_lssolve_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
610 gsl_test(f, " QR_lssolve m(5,3)");
611 s += f;
612
613 f = test_QR_lssolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
614 gsl_test(f, " QR_lssolve hilbert(2)");
615 s += f;
616
617 f = test_QR_lssolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
618 gsl_test(f, " QR_lssolve hilbert(3)");
619 s += f;
620
621 f = test_QR_lssolve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
622 gsl_test(f, " QR_lssolve hilbert(4)");
623 s += f;
624
625 f = test_QR_lssolve_dim(hilb12, hilb12_solution, 0.5);
626 gsl_test(f, " QR_lssolve hilbert(12)");
627 s += f;
628
629 f = test_QR_lssolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
630 gsl_test(f, " QR_lssolve vander(2)");
631 s += f;
632
633 f = test_QR_lssolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
634 gsl_test(f, " QR_lssolve vander(3)");
635 s += f;
636
637 f = test_QR_lssolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
638 gsl_test(f, " QR_lssolve vander(4)");
639 s += f;
640
641 f = test_QR_lssolve_dim(vander12, vander12_solution, 0.05);
642 gsl_test(f, " QR_lssolve vander(12)");
643 s += f;
644
645 return s;
646 }
647
648 int
test_QRPT_lssolve_dim(const gsl_matrix * m,const double * actual,double eps)649 test_QRPT_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
650 {
651 int s = 0;
652 size_t i, M = m->size1, N = m->size2;
653
654 gsl_vector * rhs = gsl_vector_alloc(M);
655 gsl_matrix * QRPT = gsl_matrix_alloc(M,N);
656 gsl_vector * tau = gsl_vector_alloc(N);
657 gsl_vector * x = gsl_vector_alloc(N);
658 gsl_vector * r = gsl_vector_alloc(M);
659 gsl_vector * res = gsl_vector_alloc(M);
660 gsl_permutation * perm = gsl_permutation_alloc(N);
661 gsl_vector * work = gsl_vector_alloc(N);
662 int signum;
663
664 gsl_matrix_memcpy(QRPT, m);
665
666 for(i = 0; i < M; i++)
667 gsl_vector_set(rhs, i, i + 1.0);
668
669 s += gsl_linalg_QRPT_decomp(QRPT, tau, perm, &signum, work);
670 s += gsl_linalg_QRPT_lssolve(QRPT, tau, perm, rhs, x, res);
671
672 for (i = 0; i < N; i++)
673 {
674 int foo = check(gsl_vector_get(x, i), actual[i], eps);
675 if(foo)
676 {
677 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
678 }
679 s += foo;
680 }
681
682 /* compute residual r = b - m x */
683 if (M == N)
684 {
685 gsl_vector_set_zero(r);
686 }
687 else
688 {
689 gsl_vector_memcpy(r, rhs);
690 gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
691 }
692
693 for (i = 0; i < N; i++)
694 {
695 int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
696 if(foo)
697 {
698 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
699 }
700 s += foo;
701 }
702
703 gsl_vector_free(r);
704 gsl_vector_free(res);
705 gsl_vector_free(x);
706 gsl_vector_free(tau);
707 gsl_matrix_free(QRPT);
708 gsl_vector_free(rhs);
709 gsl_permutation_free(perm);
710 gsl_vector_free(work);
711
712 return s;
713 }
714
715 int
test_QRPT_lssolve(void)716 test_QRPT_lssolve(void)
717 {
718 int f;
719 int s = 0;
720
721 f = test_QRPT_lssolve_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
722 gsl_test(f, " QRPT_lssolve m(5,3)");
723 s += f;
724
725 f = test_QRPT_lssolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
726 gsl_test(f, " QRPT_lssolve hilbert(2)");
727 s += f;
728
729 f = test_QRPT_lssolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
730 gsl_test(f, " QRPT_lssolve hilbert(3)");
731 s += f;
732
733 f = test_QRPT_lssolve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
734 gsl_test(f, " QRPT_lssolve hilbert(4)");
735 s += f;
736
737 f = test_QRPT_lssolve_dim(hilb12, hilb12_solution, 0.5);
738 gsl_test(f, " QRPT_lssolve hilbert(12)");
739 s += f;
740
741 f = test_QRPT_lssolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
742 gsl_test(f, " QRPT_lssolve vander(2)");
743 s += f;
744
745 f = test_QRPT_lssolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
746 gsl_test(f, " QRPT_lssolve vander(3)");
747 s += f;
748
749 f = test_QRPT_lssolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
750 gsl_test(f, " QRPT_lssolve vander(4)");
751 s += f;
752
753 f = test_QRPT_lssolve_dim(vander12, vander12_solution, 0.05);
754 gsl_test(f, " QRPT_lssolve vander(12)");
755 s += f;
756
757 return s;
758 }
759
760 int
test_QRPT_lssolve2_dim(const gsl_matrix * m,const double * actual,double eps)761 test_QRPT_lssolve2_dim(const gsl_matrix * m, const double * actual, double eps)
762 {
763 int s = 0;
764 size_t i, M = m->size1, N = m->size2;
765
766 gsl_vector * rhs = gsl_vector_alloc(M);
767 gsl_matrix * QRPT = gsl_matrix_alloc(M,N);
768 gsl_vector * tau = gsl_vector_alloc(N);
769 gsl_vector * x = gsl_vector_alloc(N);
770 gsl_vector * r = gsl_vector_alloc(M);
771 gsl_vector * res = gsl_vector_alloc(M);
772 gsl_permutation * perm = gsl_permutation_alloc(N);
773 gsl_vector * work = gsl_vector_alloc(N);
774 int signum;
775 size_t rank;
776
777 gsl_matrix_memcpy(QRPT, m);
778
779 for(i = 0; i < M; i++)
780 gsl_vector_set(rhs, i, i + 1.0);
781
782 s += gsl_linalg_QRPT_decomp(QRPT, tau, perm, &signum, work);
783
784 rank = gsl_linalg_QRPT_rank(QRPT, -1.0);
785
786 s += gsl_linalg_QRPT_lssolve2(QRPT, tau, perm, rhs, rank, x, res);
787
788 if (M > N)
789 {
790 for (i = 0; i < N; i++)
791 {
792 int foo = check(gsl_vector_get(x, i), actual[i], eps);
793 if(foo)
794 {
795 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
796 }
797 s += foo;
798 }
799 }
800
801 /* compute residual r = b - m x */
802 if (M == N)
803 {
804 gsl_vector_set_zero(r);
805 }
806 else
807 {
808 gsl_vector_memcpy(r, rhs);
809 gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
810 }
811
812 for (i = 0; i < N; i++)
813 {
814 int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
815 if(foo)
816 {
817 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
818 }
819 s += foo;
820 }
821
822 gsl_vector_free(r);
823 gsl_vector_free(res);
824 gsl_vector_free(x);
825 gsl_vector_free(tau);
826 gsl_matrix_free(QRPT);
827 gsl_vector_free(rhs);
828 gsl_permutation_free(perm);
829 gsl_vector_free(work);
830
831 return s;
832 }
833
834 int
test_QRPT_lssolve2(void)835 test_QRPT_lssolve2(void)
836 {
837 int f;
838 int s = 0;
839
840 f = test_QRPT_lssolve2_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
841 gsl_test(f, " QRPT_lssolve2 m(5,3)");
842 s += f;
843
844 f = test_QRPT_lssolve2_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
845 gsl_test(f, " QRPT_lssolve2 hilbert(2)");
846 s += f;
847
848 f = test_QRPT_lssolve2_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
849 gsl_test(f, " QRPT_lssolve2 hilbert(3)");
850 s += f;
851
852 f = test_QRPT_lssolve2_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
853 gsl_test(f, " QRPT_lssolve2 hilbert(4)");
854 s += f;
855
856 f = test_QRPT_lssolve2_dim(hilb12, hilb12_solution, 0.5);
857 gsl_test(f, " QRPT_lssolve2 hilbert(12)");
858 s += f;
859
860 f = test_QRPT_lssolve2_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
861 gsl_test(f, " QRPT_lssolve2 vander(2)");
862 s += f;
863
864 f = test_QRPT_lssolve2_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
865 gsl_test(f, " QRPT_lssolve2 vander(3)");
866 s += f;
867
868 f = test_QRPT_lssolve2_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
869 gsl_test(f, " QRPT_lssolve2 vander(4)");
870 s += f;
871
872 f = test_QRPT_lssolve2_dim(vander12, vander12_solution, 0.05);
873 gsl_test(f, " QRPT_lssolve2 vander(12)");
874 s += f;
875
876 return s;
877 }
878
879 int
test_QRPT_solve_dim(const gsl_matrix * m,const double * actual,double eps)880 test_QRPT_solve_dim(const gsl_matrix * m, const double * actual, double eps)
881 {
882 int s = 0;
883 int signum;
884 unsigned long i, dim = m->size1;
885
886 gsl_permutation * perm = gsl_permutation_alloc(dim);
887 gsl_vector * rhs = gsl_vector_alloc(dim);
888 gsl_matrix * qr = gsl_matrix_alloc(dim,dim);
889 gsl_vector * d = gsl_vector_alloc(dim);
890 gsl_vector * x = gsl_vector_alloc(dim);
891 gsl_vector * norm = gsl_vector_alloc(dim);
892
893 gsl_matrix_memcpy(qr,m);
894 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
895 s += gsl_linalg_QRPT_decomp(qr, d, perm, &signum, norm);
896 s += gsl_linalg_QRPT_solve(qr, d, perm, rhs, x);
897 for(i=0; i<dim; i++) {
898 int foo = check(gsl_vector_get(x, i), actual[i], eps);
899 if(foo) {
900 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
901 }
902 s += foo;
903 }
904
905 gsl_vector_free(norm);
906 gsl_vector_free(x);
907 gsl_vector_free(d);
908 gsl_matrix_free(qr);
909 gsl_vector_free(rhs);
910 gsl_permutation_free(perm);
911
912 return s;
913 }
914
test_QRPT_solve(void)915 int test_QRPT_solve(void)
916 {
917 int f;
918 int s = 0;
919
920 f = test_QRPT_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
921 gsl_test(f, " QRPT_solve hilbert(2)");
922 s += f;
923
924 f = test_QRPT_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
925 gsl_test(f, " QRPT_solve hilbert(3)");
926 s += f;
927
928 f = test_QRPT_solve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
929 gsl_test(f, " QRPT_solve hilbert(4)");
930 s += f;
931
932 f = test_QRPT_solve_dim(hilb12, hilb12_solution, 0.5);
933 gsl_test(f, " QRPT_solve hilbert(12)");
934 s += f;
935
936 f = test_QRPT_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
937 gsl_test(f, " QRPT_solve vander(2)");
938 s += f;
939
940 f = test_QRPT_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
941 gsl_test(f, " QRPT_solve vander(3)");
942 s += f;
943
944 f = test_QRPT_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
945 gsl_test(f, " QRPT_solve vander(4)");
946 s += f;
947
948 f = test_QRPT_solve_dim(vander12, vander12_solution, 0.05);
949 gsl_test(f, " QRPT_solve vander(12)");
950 s += f;
951
952 return s;
953 }
954
955 int
test_QRPT_QRsolve_dim(const gsl_matrix * m,const double * actual,double eps)956 test_QRPT_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps)
957 {
958 int s = 0;
959 int signum;
960 unsigned long i, dim = m->size1;
961
962 gsl_permutation * perm = gsl_permutation_alloc(dim);
963 gsl_vector * rhs = gsl_vector_alloc(dim);
964 gsl_matrix * qr = gsl_matrix_alloc(dim,dim);
965 gsl_matrix * q = gsl_matrix_alloc(dim,dim);
966 gsl_matrix * r = gsl_matrix_alloc(dim,dim);
967 gsl_vector * d = gsl_vector_alloc(dim);
968 gsl_vector * x = gsl_vector_alloc(dim);
969 gsl_vector * norm = gsl_vector_alloc(dim);
970
971 gsl_matrix_memcpy(qr,m);
972 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
973 s += gsl_linalg_QRPT_decomp2(qr, q, r, d, perm, &signum, norm);
974 s += gsl_linalg_QRPT_QRsolve(q, r, perm, rhs, x);
975 for(i=0; i<dim; i++) {
976 int foo = check(gsl_vector_get(x, i), actual[i], eps);
977 if(foo) {
978 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
979 }
980 s += foo;
981 }
982
983 gsl_vector_free(norm);
984 gsl_vector_free(x);
985 gsl_vector_free(d);
986 gsl_matrix_free(qr);
987 gsl_matrix_free(q);
988 gsl_matrix_free(r);
989 gsl_vector_free(rhs);
990 gsl_permutation_free(perm);
991
992 return s;
993 }
994
test_QRPT_QRsolve(void)995 int test_QRPT_QRsolve(void)
996 {
997 int f;
998 int s = 0;
999
1000 f = test_QRPT_QRsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1001 gsl_test(f, " QRPT_QRsolve hilbert(2)");
1002 s += f;
1003
1004 f = test_QRPT_QRsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1005 gsl_test(f, " QRPT_QRsolve hilbert(3)");
1006 s += f;
1007
1008 f = test_QRPT_QRsolve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
1009 gsl_test(f, " QRPT_QRsolve hilbert(4)");
1010 s += f;
1011
1012 f = test_QRPT_QRsolve_dim(hilb12, hilb12_solution, 0.5);
1013 gsl_test(f, " QRPT_QRsolve hilbert(12)");
1014 s += f;
1015
1016 f = test_QRPT_QRsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1017 gsl_test(f, " QRPT_QRsolve vander(2)");
1018 s += f;
1019
1020 f = test_QRPT_QRsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1021 gsl_test(f, " QRPT_QRsolve vander(3)");
1022 s += f;
1023
1024 f = test_QRPT_QRsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1025 gsl_test(f, " QRPT_QRsolve vander(4)");
1026 s += f;
1027
1028 f = test_QRPT_QRsolve_dim(vander12, vander12_solution, 0.05);
1029 gsl_test(f, " QRPT_QRsolve vander(12)");
1030 s += f;
1031
1032 return s;
1033 }
1034
1035 int
test_QRPT_decomp_dim(const gsl_matrix * m,const double expected_rcond,double eps)1036 test_QRPT_decomp_dim(const gsl_matrix * m, const double expected_rcond, double eps)
1037 {
1038 int s = 0, signum;
1039 unsigned long i,j, M = m->size1, N = m->size2;
1040
1041 gsl_matrix * QR = gsl_matrix_alloc(M, N);
1042 gsl_matrix * A = gsl_matrix_alloc(M, N);
1043 gsl_matrix * Q = gsl_matrix_alloc(M, M);
1044 gsl_matrix * R = gsl_matrix_alloc(M, N);
1045 gsl_vector * tau = gsl_vector_alloc(GSL_MIN(M, N));
1046 gsl_vector * norm = gsl_vector_alloc(N);
1047
1048 gsl_permutation * perm = gsl_permutation_alloc(N);
1049
1050 gsl_matrix_memcpy(QR, m);
1051 s += gsl_linalg_QRPT_decomp(QR, tau, perm, &signum, norm);
1052 s += gsl_linalg_QR_unpack(QR, tau, Q, R);
1053
1054 /* compute A = Q R */
1055 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, Q, R, 0.0, A);
1056
1057 /* Compute QR P^T by permuting the elements of the rows of QR */
1058
1059 for (i = 0; i < M; i++)
1060 {
1061 gsl_vector_view row = gsl_matrix_row (A, i);
1062 gsl_permute_vector_inverse (perm, &row.vector);
1063 }
1064
1065 for (i = 0; i < M; i++)
1066 {
1067 for (j = 0; j < N; j++)
1068 {
1069 double aij = gsl_matrix_get(A, i, j);
1070 double mij = gsl_matrix_get(m, i, j);
1071 int foo = check(aij, mij, eps);
1072 if(foo)
1073 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
1074 s += foo;
1075 }
1076 }
1077
1078 if (expected_rcond > 0.0)
1079 {
1080 double rcond;
1081 int foo;
1082 gsl_vector * work = gsl_vector_alloc(3 * N);
1083
1084 gsl_linalg_QRPT_rcond(QR, &rcond, work);
1085 foo = check(rcond, expected_rcond, 1.0e-10);
1086 if (foo)
1087 printf("QRPT_rcond (%3lu,%3lu): %22.18g %22.18g\n", M, N, rcond, expected_rcond);
1088
1089 s += foo;
1090
1091 gsl_vector_free(work);
1092 }
1093
1094 gsl_permutation_free (perm);
1095 gsl_vector_free(norm);
1096 gsl_vector_free(tau);
1097 gsl_matrix_free(QR);
1098 gsl_matrix_free(A);
1099 gsl_matrix_free(Q);
1100 gsl_matrix_free(R);
1101
1102 return s;
1103 }
1104
test_QRPT_decomp(void)1105 int test_QRPT_decomp(void)
1106 {
1107 int f;
1108 int s = 0;
1109
1110 f = test_QRPT_decomp_dim(m35, -1.0, 2 * 8.0 * GSL_DBL_EPSILON);
1111 gsl_test(f, " QRPT_decomp m(3,5)");
1112 s += f;
1113
1114 /* rcond value from LAPACK DTRCON */
1115 f = test_QRPT_decomp_dim(m53, 2.915362697820e-03, 2 * 8.0 * GSL_DBL_EPSILON);
1116 gsl_test(f, " QRPT_decomp m(5,3)");
1117 s += f;
1118
1119 f = test_QRPT_decomp_dim(s35, -1.0, 2 * 8.0 * GSL_DBL_EPSILON);
1120 gsl_test(f, " QRPT_decomp s(3,5)");
1121 s += f;
1122
1123 f = test_QRPT_decomp_dim(s53, 1.002045825443827e-03, 2 * 8.0 * GSL_DBL_EPSILON);
1124 gsl_test(f, " QRPT_decomp s(5,3)");
1125 s += f;
1126
1127 f = test_QRPT_decomp_dim(hilb2, 4.347826086956522e-02, 2 * 8.0 * GSL_DBL_EPSILON);
1128 gsl_test(f, " QRPT_decomp hilbert(2)");
1129 s += f;
1130
1131 f = test_QRPT_decomp_dim(hilb3, 1.505488055305100e-03, 2 * 64.0 * GSL_DBL_EPSILON);
1132 gsl_test(f, " QRPT_decomp hilbert(3)");
1133 s += f;
1134
1135 f = test_QRPT_decomp_dim(hilb4, 4.872100915910022e-05, 2 * 1024.0 * GSL_DBL_EPSILON);
1136 gsl_test(f, " QRPT_decomp hilbert(4)");
1137 s += f;
1138
1139 f = test_QRPT_decomp_dim(hilb12, -1.0, 2 * 1024.0 * GSL_DBL_EPSILON);
1140 gsl_test(f, " QRPT_decomp hilbert(12)");
1141 s += f;
1142
1143 f = test_QRPT_decomp_dim(vander2, 1.249999999999999e-01, 8.0 * GSL_DBL_EPSILON);
1144 gsl_test(f, " QRPT_decomp vander(2)");
1145 s += f;
1146
1147 f = test_QRPT_decomp_dim(vander3, 9.708737864077667e-03, 64.0 * GSL_DBL_EPSILON);
1148 gsl_test(f, " QRPT_decomp vander(3)");
1149 s += f;
1150
1151 f = test_QRPT_decomp_dim(vander4, 5.255631229339451e-04, 1024.0 * GSL_DBL_EPSILON);
1152 gsl_test(f, " QRPT_decomp vander(4)");
1153 s += f;
1154
1155 f = test_QRPT_decomp_dim(vander12, -1.0, 0.0005); /* FIXME: bad accuracy */
1156 gsl_test(f, " QRPT_decomp vander(12)");
1157 s += f;
1158
1159 return s;
1160 }
1161
1162
1163 int
test_QR_update_dim(const gsl_matrix * m,double eps)1164 test_QR_update_dim(const gsl_matrix * m, double eps)
1165 {
1166 int s = 0;
1167 unsigned long i,j,k, M = m->size1, N = m->size2;
1168
1169 gsl_matrix * qr1 = gsl_matrix_alloc(M,N);
1170 gsl_matrix * qr2 = gsl_matrix_alloc(M,N);
1171 gsl_matrix * q1 = gsl_matrix_alloc(M,M);
1172 gsl_matrix * r1 = gsl_matrix_alloc(M,N);
1173 gsl_matrix * q2 = gsl_matrix_alloc(M,M);
1174 gsl_matrix * r2 = gsl_matrix_alloc(M,N);
1175 gsl_vector * d = gsl_vector_alloc(N);
1176 gsl_vector * solution1 = gsl_vector_alloc(N);
1177 gsl_vector * solution2 = gsl_vector_alloc(N);
1178 gsl_vector * u = gsl_vector_alloc(M);
1179 gsl_vector * v = gsl_vector_alloc(N);
1180 gsl_vector * w = gsl_vector_alloc(M);
1181
1182 gsl_matrix_memcpy(qr1,m);
1183 gsl_matrix_memcpy(qr2,m);
1184
1185 for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
1186 for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
1187
1188 for(i=0; i<M; i++)
1189 {
1190 double ui = gsl_vector_get(u, i);
1191 for(j=0; j<N; j++)
1192 {
1193 double vj = gsl_vector_get(v, j);
1194 double qij = gsl_matrix_get(qr1, i, j);
1195 gsl_matrix_set(qr1, i, j, qij + ui * vj);
1196 }
1197 }
1198
1199 s += gsl_linalg_QR_decomp(qr2, d);
1200 s += gsl_linalg_QR_unpack(qr2, d, q2, r2);
1201
1202 /* compute w = Q^T u */
1203
1204 for (j = 0; j < M; j++)
1205 {
1206 double sum = 0;
1207 for (i = 0; i < M; i++)
1208 sum += gsl_matrix_get (q2, i, j) * gsl_vector_get (u, i);
1209 gsl_vector_set (w, j, sum);
1210 }
1211
1212 s += gsl_linalg_QR_update(q2, r2, w, v);
1213
1214 /* compute qr2 = q2 * r2 */
1215
1216 for (i = 0; i < M; i++)
1217 {
1218 for (j = 0; j< N; j++)
1219 {
1220 double sum = 0;
1221 for (k = 0; k <= GSL_MIN(j,M-1); k++)
1222 {
1223 double qik = gsl_matrix_get(q2, i, k);
1224 double rkj = gsl_matrix_get(r2, k, j);
1225 sum += qik * rkj ;
1226 }
1227 gsl_matrix_set (qr2, i, j, sum);
1228 }
1229 }
1230
1231 for(i=0; i<M; i++) {
1232 for(j=0; j<N; j++) {
1233 double s1 = gsl_matrix_get(qr1, i, j);
1234 double s2 = gsl_matrix_get(qr2, i, j);
1235
1236 int foo = check(s1, s2, eps);
1237 if(foo) {
1238 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, s1, s2);
1239 }
1240 s += foo;
1241 }
1242 }
1243
1244 gsl_vector_free(solution1);
1245 gsl_vector_free(solution2);
1246 gsl_vector_free(d);
1247 gsl_vector_free(u);
1248 gsl_vector_free(v);
1249 gsl_vector_free(w);
1250 gsl_matrix_free(qr1);
1251 gsl_matrix_free(qr2);
1252 gsl_matrix_free(q1);
1253 gsl_matrix_free(r1);
1254 gsl_matrix_free(q2);
1255 gsl_matrix_free(r2);
1256
1257 return s;
1258 }
1259
test_QR_update(void)1260 int test_QR_update(void)
1261 {
1262 int f;
1263 int s = 0;
1264
1265 f = test_QR_update_dim(m35, 2 * 512.0 * GSL_DBL_EPSILON);
1266 gsl_test(f, " QR_update m(3,5)");
1267 s += f;
1268
1269 f = test_QR_update_dim(m53, 2 * 512.0 * GSL_DBL_EPSILON);
1270 gsl_test(f, " QR_update m(5,3)");
1271 s += f;
1272
1273 f = test_QR_update_dim(hilb2, 2 * 512.0 * GSL_DBL_EPSILON);
1274 gsl_test(f, " QR_update hilbert(2)");
1275 s += f;
1276
1277 f = test_QR_update_dim(hilb3, 2 * 512.0 * GSL_DBL_EPSILON);
1278 gsl_test(f, " QR_update hilbert(3)");
1279 s += f;
1280
1281 f = test_QR_update_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1282 gsl_test(f, " QR_update hilbert(4)");
1283 s += f;
1284
1285 f = test_QR_update_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1286 gsl_test(f, " QR_update hilbert(12)");
1287 s += f;
1288
1289 f = test_QR_update_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1290 gsl_test(f, " QR_update vander(2)");
1291 s += f;
1292
1293 f = test_QR_update_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1294 gsl_test(f, " QR_update vander(3)");
1295 s += f;
1296
1297 f = test_QR_update_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1298 gsl_test(f, " QR_update vander(4)");
1299 s += f;
1300
1301 f = test_QR_update_dim(vander12, 0.0005); /* FIXME: bad accuracy */
1302 gsl_test(f, " QR_update vander(12)");
1303 s += f;
1304
1305 return s;
1306 }
1307
1308
1309 int
test_QRPT_update_dim(const gsl_matrix * m,double eps)1310 test_QRPT_update_dim(const gsl_matrix * m, double eps)
1311 {
1312 int s = 0, signum;
1313 unsigned long i,j,k, M = m->size1, N = m->size2;
1314
1315 gsl_matrix * qr1 = gsl_matrix_alloc(M,N);
1316 gsl_matrix * qr2 = gsl_matrix_alloc(M,N);
1317 gsl_matrix * q1 = gsl_matrix_alloc(M,M);
1318 gsl_matrix * r1 = gsl_matrix_alloc(M,N);
1319 gsl_matrix * q2 = gsl_matrix_alloc(M,M);
1320 gsl_matrix * r2 = gsl_matrix_alloc(M,N);
1321 gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
1322 gsl_vector * u = gsl_vector_alloc(M);
1323 gsl_vector * v = gsl_vector_alloc(N);
1324 gsl_vector * w = gsl_vector_alloc(M);
1325
1326 gsl_vector * norm = gsl_vector_alloc(N);
1327 gsl_permutation * perm = gsl_permutation_alloc(N);
1328
1329 gsl_matrix_memcpy(qr1,m);
1330 gsl_matrix_memcpy(qr2,m);
1331 for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
1332 for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
1333
1334 for(i=0; i<M; i++)
1335 {
1336 double ui = gsl_vector_get(u, i);
1337 for(j=0; j<N; j++)
1338 {
1339 double vj = gsl_vector_get(v, j);
1340 double qij = gsl_matrix_get(qr1, i, j);
1341 gsl_matrix_set(qr1, i, j, qij + ui * vj);
1342 }
1343 }
1344
1345 s += gsl_linalg_QRPT_decomp(qr2, d, perm, &signum, norm);
1346 s += gsl_linalg_QR_unpack(qr2, d, q2, r2);
1347
1348 /* compute w = Q^T u */
1349
1350 for (j = 0; j < M; j++)
1351 {
1352 double sum = 0;
1353 for (i = 0; i < M; i++)
1354 sum += gsl_matrix_get (q2, i, j) * gsl_vector_get (u, i);
1355 gsl_vector_set (w, j, sum);
1356 }
1357
1358 s += gsl_linalg_QRPT_update(q2, r2, perm, w, v);
1359
1360 /* Now compute qr2 = q2 * r2 * p^T */
1361
1362 /* first multiply q2 * r2 */
1363
1364 for (i = 0; i < M; i++)
1365 {
1366 for (j = 0; j< N; j++)
1367 {
1368 double sum = 0;
1369 for (k = 0; k <= GSL_MIN(j,M-1); k++)
1370 {
1371 double qik = gsl_matrix_get(q2, i, k);
1372 double rkj = gsl_matrix_get(r2, k, j);
1373 sum += qik * rkj ;
1374 }
1375 gsl_matrix_set (qr2, i, j, sum);
1376 }
1377 }
1378
1379 /* now apply permutation to get qr2 = q2 * r2 * p^T */
1380
1381 for (i = 0; i < M ; i++)
1382 {
1383 gsl_vector_view r_i = gsl_matrix_row(qr2, i);
1384 gsl_permute_vector_inverse(perm, &r_i.vector);
1385 }
1386
1387
1388 for(i=0; i<M; i++) {
1389 for(j=0; j<N; j++) {
1390 double s1 = gsl_matrix_get(qr1, i, j);
1391 double s2 = gsl_matrix_get(qr2, i, j);
1392
1393 int foo = check(s1, s2, eps);
1394 if(foo) {
1395 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, s1, s2);
1396 }
1397 s += foo;
1398 }
1399 }
1400
1401 gsl_permutation_free (perm);
1402 gsl_vector_free(norm);
1403 gsl_vector_free(d);
1404 gsl_vector_free(u);
1405 gsl_vector_free(v);
1406 gsl_vector_free(w);
1407 gsl_matrix_free(qr1);
1408 gsl_matrix_free(qr2);
1409 gsl_matrix_free(q1);
1410 gsl_matrix_free(r1);
1411 gsl_matrix_free(q2);
1412 gsl_matrix_free(r2);
1413
1414 return s;
1415 }
1416
test_QRPT_update(void)1417 int test_QRPT_update(void)
1418 {
1419 int f;
1420 int s = 0;
1421
1422 f = test_QRPT_update_dim(m35, 2 * 512.0 * GSL_DBL_EPSILON);
1423 gsl_test(f, " QRPT_update m(3,5)");
1424 s += f;
1425
1426 f = test_QRPT_update_dim(m53, 2 * 1024.0 * GSL_DBL_EPSILON);
1427 gsl_test(f, " QRPT_update m(5,3)");
1428 s += f;
1429
1430 f = test_QRPT_update_dim(hilb2, 2 * 512.0 * GSL_DBL_EPSILON);
1431 gsl_test(f, " QRPT_update hilbert(2)");
1432 s += f;
1433
1434 f = test_QRPT_update_dim(hilb3, 2 * 512.0 * GSL_DBL_EPSILON);
1435 gsl_test(f, " QRPT_update hilbert(3)");
1436 s += f;
1437
1438 f = test_QRPT_update_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1439 gsl_test(f, " QRPT_update hilbert(4)");
1440 s += f;
1441
1442 f = test_QRPT_update_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1443 gsl_test(f, " QRPT_update hilbert(12)");
1444 s += f;
1445
1446 f = test_QRPT_update_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1447 gsl_test(f, " QRPT_update vander(2)");
1448 s += f;
1449
1450 f = test_QRPT_update_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1451 gsl_test(f, " QRPT_update vander(3)");
1452 s += f;
1453
1454 f = test_QRPT_update_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1455 gsl_test(f, " QRPT_update vander(4)");
1456 s += f;
1457
1458 f = test_QRPT_update_dim(vander12, 0.0005); /* FIXME: bad accuracy */
1459 gsl_test(f, " QRPT_update vander(12)");
1460 s += f;
1461
1462 return s;
1463 }
1464
1465 int
test_SV_solve_dim(const gsl_matrix * m,const double * actual,double eps)1466 test_SV_solve_dim(const gsl_matrix * m, const double * actual, double eps)
1467 {
1468 int s = 0;
1469 unsigned long i, dim = m->size1;
1470
1471 gsl_vector * rhs = gsl_vector_alloc(dim);
1472 gsl_matrix * u = gsl_matrix_alloc(dim,dim);
1473 gsl_matrix * q = gsl_matrix_alloc(dim,dim);
1474 gsl_vector * d = gsl_vector_alloc(dim);
1475 gsl_vector * x = gsl_vector_calloc(dim);
1476 gsl_matrix_memcpy(u,m);
1477 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
1478 s += gsl_linalg_SV_decomp(u, q, d, x);
1479 s += gsl_linalg_SV_solve(u, q, d, rhs, x);
1480 for(i=0; i<dim; i++) {
1481 int foo = check(gsl_vector_get(x, i), actual[i], eps);
1482 if(foo) {
1483 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
1484 }
1485 s += foo;
1486 }
1487 gsl_vector_free(x);
1488 gsl_vector_free(d);
1489 gsl_matrix_free(u);
1490 gsl_matrix_free(q);
1491 gsl_vector_free(rhs);
1492
1493 return s;
1494 }
1495
test_SV_solve(void)1496 int test_SV_solve(void)
1497 {
1498 int f;
1499 int s = 0;
1500
1501 f = test_SV_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1502 gsl_test(f, " SV_solve hilbert(2)");
1503 s += f;
1504
1505 f = test_SV_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1506 gsl_test(f, " SV_solve hilbert(3)");
1507 s += f;
1508
1509 f = test_SV_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
1510 gsl_test(f, " SV_solve hilbert(4)");
1511 s += f;
1512
1513 f = test_SV_solve_dim(hilb12, hilb12_solution, 0.5);
1514 gsl_test(f, " SV_solve hilbert(12)");
1515 s += f;
1516
1517 f = test_SV_solve_dim(vander2, vander2_solution, 64.0 * GSL_DBL_EPSILON);
1518 gsl_test(f, " SV_solve vander(2)");
1519 s += f;
1520
1521 f = test_SV_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1522 gsl_test(f, " SV_solve vander(3)");
1523 s += f;
1524
1525 f = test_SV_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1526 gsl_test(f, " SV_solve vander(4)");
1527 s += f;
1528
1529 f = test_SV_solve_dim(vander12, vander12_solution, 0.05);
1530 gsl_test(f, " SV_solve vander(12)");
1531 s += f;
1532
1533 return s;
1534 }
1535
1536 int
test_SV_decomp_dim(const gsl_matrix * m,double eps)1537 test_SV_decomp_dim(const gsl_matrix * m, double eps)
1538 {
1539 int s = 0;
1540 double di1;
1541 unsigned long i,j, M = m->size1, N = m->size2;
1542 unsigned long input_nans = 0;
1543
1544 gsl_matrix * v = gsl_matrix_alloc(M,N);
1545 gsl_matrix * a = gsl_matrix_alloc(M,N);
1546 gsl_matrix * q = gsl_matrix_alloc(N,N);
1547 gsl_matrix * dqt = gsl_matrix_alloc(N,N);
1548 gsl_vector * d = gsl_vector_alloc(N);
1549 gsl_vector * w = gsl_vector_alloc(N);
1550
1551 gsl_matrix_memcpy(v,m);
1552
1553 /* Check for nans in the input */
1554 for (i = 0; i<M; i++) {
1555 for (j = 0; j<N; j++) {
1556 double m_ij = gsl_matrix_get (m, i, j);
1557 if (gsl_isnan (m_ij)) input_nans++;
1558 }
1559 }
1560
1561 s = gsl_linalg_SV_decomp(v, q, d, w);
1562
1563 if (s) printf("returned error code %d = %s\n", s, gsl_strerror(s));
1564
1565 /* Check that singular values are non-negative and in non-decreasing
1566 order */
1567
1568 di1 = 0.0;
1569
1570 for (i = 0; i < N; i++)
1571 {
1572 double di = gsl_vector_get (d, i);
1573
1574 if (gsl_isnan (di))
1575 {
1576 if (input_nans > 0)
1577 continue; /* skip NaNs if present in input */
1578 else
1579 {
1580 s++;
1581 printf("bad singular value %lu = %22.18g\n", i, di);
1582 }
1583 }
1584
1585 if (di < 0) {
1586 s++;
1587 printf("singular value %lu = %22.18g < 0\n", i, di);
1588 }
1589
1590 if(i > 0 && di > di1) {
1591 s++;
1592 printf("singular value %lu = %22.18g vs previous %22.18g\n", i, di, di1);
1593 }
1594
1595 di1 = di;
1596 }
1597
1598 /* Scale dqt = D Q^T */
1599
1600 for (i = 0; i < N ; i++)
1601 {
1602 double di = gsl_vector_get (d, i);
1603
1604 for (j = 0; j < N; j++)
1605 {
1606 double qji = gsl_matrix_get(q, j, i);
1607 gsl_matrix_set (dqt, i, j, qji * di);
1608 }
1609 }
1610
1611 /* compute a = v dqt */
1612 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
1613
1614 for(i=0; i<M; i++) {
1615 for(j=0; j<N; j++) {
1616 double aij = gsl_matrix_get(a, i, j);
1617 double mij = gsl_matrix_get(m, i, j);
1618 int foo = check(aij, mij, eps);
1619 if(foo) {
1620 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
1621 }
1622 s += foo;
1623 }
1624 }
1625 gsl_vector_free(w);
1626 gsl_vector_free(d);
1627 gsl_matrix_free(v);
1628 gsl_matrix_free(a);
1629 gsl_matrix_free(q);
1630 gsl_matrix_free(dqt);
1631
1632 return s;
1633 }
1634
test_SV_decomp(void)1635 int test_SV_decomp(void)
1636 {
1637 int f;
1638 int s = 0;
1639
1640 f = test_SV_decomp_dim(m11, 2 * GSL_DBL_EPSILON);
1641 gsl_test(f, " SV_decomp m(1,1)");
1642 s += f;
1643
1644 f = test_SV_decomp_dim(m51, 2 * 64.0 * GSL_DBL_EPSILON);
1645 gsl_test(f, " SV_decomp m(5,1)");
1646 s += f;
1647
1648 /* M<N not implemented yet */
1649 #if 0
1650 f = test_SV_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
1651 gsl_test(f, " SV_decomp m(3,5)");
1652 s += f;
1653 #endif
1654 f = test_SV_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
1655 gsl_test(f, " SV_decomp m(5,3)");
1656 s += f;
1657
1658 f = test_SV_decomp_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
1659 gsl_test(f, " SV_decomp moler(10)");
1660 s += f;
1661
1662 f = test_SV_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
1663 gsl_test(f, " SV_decomp hilbert(2)");
1664 s += f;
1665
1666 f = test_SV_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
1667 gsl_test(f, " SV_decomp hilbert(3)");
1668 s += f;
1669
1670 f = test_SV_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1671 gsl_test(f, " SV_decomp hilbert(4)");
1672 s += f;
1673
1674 f = test_SV_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1675 gsl_test(f, " SV_decomp hilbert(12)");
1676 s += f;
1677
1678 f = test_SV_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1679 gsl_test(f, " SV_decomp vander(2)");
1680 s += f;
1681
1682 f = test_SV_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1683 gsl_test(f, " SV_decomp vander(3)");
1684 s += f;
1685
1686 f = test_SV_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1687 gsl_test(f, " SV_decomp vander(4)");
1688 s += f;
1689
1690 f = test_SV_decomp_dim(vander12, 1e-4);
1691 gsl_test(f, " SV_decomp vander(12)");
1692 s += f;
1693
1694 f = test_SV_decomp_dim(row3, 10 * GSL_DBL_EPSILON);
1695 gsl_test(f, " SV_decomp row3");
1696 s += f;
1697
1698 f = test_SV_decomp_dim(row5, 128 * GSL_DBL_EPSILON);
1699 gsl_test(f, " SV_decomp row5");
1700 s += f;
1701
1702 f = test_SV_decomp_dim(row12, 1024 * GSL_DBL_EPSILON);
1703 gsl_test(f, " SV_decomp row12");
1704 s += f;
1705
1706 f = test_SV_decomp_dim(inf5, 1024 * GSL_DBL_EPSILON);
1707 gsl_test(f, " SV_decomp inf5");
1708 s += f;
1709
1710 f = test_SV_decomp_dim(nan5, 1024 * GSL_DBL_EPSILON);
1711 gsl_test(f, " SV_decomp nan5");
1712 s += f;
1713
1714 f = test_SV_decomp_dim(dblmin3, 1024 * GSL_DBL_EPSILON);
1715 gsl_test(f, " SV_decomp dblmin3");
1716 s += f;
1717
1718 f = test_SV_decomp_dim(dblmin5, 1024 * GSL_DBL_EPSILON);
1719 gsl_test(f, " SV_decomp dblmin5");
1720 s += f;
1721
1722 f = test_SV_decomp_dim(dblsubnorm5, 100 * 1024 * 1024 * GSL_DBL_EPSILON);
1723 gsl_test(f, " SV_decomp dblsubnorm5");
1724 s += f;
1725
1726 f = test_SV_decomp_dim(bigsparse, 1024 * GSL_DBL_EPSILON);
1727 gsl_test(f, " SV_decomp bigsparse");
1728 s += f;
1729
1730
1731 {
1732 double i1, i2, i3, i4;
1733 double lower = -2, upper = 2;
1734
1735 for (i1 = lower; i1 <= upper; i1++)
1736 {
1737 for (i2 = lower; i2 <= upper; i2++)
1738 {
1739 for (i3 = lower; i3 <= upper; i3++)
1740 {
1741 for (i4 = lower; i4 <= upper; i4++)
1742 {
1743 gsl_matrix_set (A22, 0,0, i1);
1744 gsl_matrix_set (A22, 0,1, i2);
1745 gsl_matrix_set (A22, 1,0, i3);
1746 gsl_matrix_set (A22, 1,1, i4);
1747
1748 f = test_SV_decomp_dim(A22, 16 * GSL_DBL_EPSILON);
1749 gsl_test(f, " SV_decomp (2x2) A=[%g, %g; %g, %g]", i1,i2,i3,i4);
1750 s += f;
1751 }
1752 }
1753 }
1754 }
1755 }
1756
1757 {
1758 int i;
1759 double carry = 0, lower = 0, upper = 1;
1760 double *a = A33->data;
1761
1762 for (i=0; i<9; i++) {
1763 a[i] = lower;
1764 }
1765
1766 while (carry == 0.0) {
1767 f = test_SV_decomp_dim(A33, 64 * GSL_DBL_EPSILON);
1768 gsl_test(f, " SV_decomp (3x3) A=[ %g, %g, %g; %g, %g, %g; %g, %g, %g]",
1769 a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]);
1770
1771 /* increment */
1772 carry=1.0;
1773 for (i=9; carry > 0.0 && i>0 && i--;)
1774 {
1775 double v=a[i]+carry;
1776 carry = (v>upper) ? 1.0 : 0.0;
1777 a[i] = (v>upper) ? lower : v;
1778 }
1779 }
1780 }
1781
1782 #ifdef TEST_SVD_4X4
1783 {
1784 int i;
1785 double carry = 0, lower = 0, upper = 1;
1786 double *a = A44->data;
1787
1788 for (i=0; i<16; i++) {
1789 a[i] = lower;
1790 }
1791
1792 while (carry == 0.0) {
1793 f = test_SV_decomp_dim(A44, 64 * GSL_DBL_EPSILON);
1794 gsl_test(f, " SV_decomp (4x4) A=[ %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g]",
1795 a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9],
1796 a[10], a[11], a[12], a[13], a[14], a[15]);
1797
1798 /* increment */
1799 carry=1.0;
1800 for (i=16; carry > 0.0 && i>0 && i--;)
1801 {
1802 double v=a[i]+carry;
1803 carry = (v>upper) ? 1.0 : 0.0;
1804 a[i] = (v>upper) ? lower : v;
1805 }
1806 }
1807 }
1808 #endif
1809
1810 return s;
1811 }
1812
1813
1814 int
test_SV_decomp_mod_dim(const gsl_matrix * m,double eps)1815 test_SV_decomp_mod_dim(const gsl_matrix * m, double eps)
1816 {
1817 int s = 0;
1818 double di1;
1819 unsigned long i,j, M = m->size1, N = m->size2;
1820
1821 gsl_matrix * v = gsl_matrix_alloc(M,N);
1822 gsl_matrix * a = gsl_matrix_alloc(M,N);
1823 gsl_matrix * q = gsl_matrix_alloc(N,N);
1824 gsl_matrix * x = gsl_matrix_alloc(N,N);
1825 gsl_matrix * dqt = gsl_matrix_alloc(N,N);
1826 gsl_vector * d = gsl_vector_alloc(N);
1827 gsl_vector * w = gsl_vector_alloc(N);
1828
1829 gsl_matrix_memcpy(v,m);
1830
1831 s += gsl_linalg_SV_decomp_mod(v, x, q, d, w);
1832
1833 /* Check that singular values are non-negative and in non-decreasing
1834 order */
1835
1836 di1 = 0.0;
1837
1838 for (i = 0; i < N; i++)
1839 {
1840 double di = gsl_vector_get (d, i);
1841
1842 if (gsl_isnan (di))
1843 {
1844 continue; /* skip NaNs */
1845 }
1846
1847 if (di < 0) {
1848 s++;
1849 printf("singular value %lu = %22.18g < 0\n", i, di);
1850 }
1851
1852 if(i > 0 && di > di1) {
1853 s++;
1854 printf("singular value %lu = %22.18g vs previous %22.18g\n", i, di, di1);
1855 }
1856
1857 di1 = di;
1858 }
1859
1860 /* Scale dqt = D Q^T */
1861
1862 for (i = 0; i < N ; i++)
1863 {
1864 double di = gsl_vector_get (d, i);
1865
1866 for (j = 0; j < N; j++)
1867 {
1868 double qji = gsl_matrix_get(q, j, i);
1869 gsl_matrix_set (dqt, i, j, qji * di);
1870 }
1871 }
1872
1873 /* compute a = v dqt */
1874 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
1875
1876 for(i=0; i<M; i++) {
1877 for(j=0; j<N; j++) {
1878 double aij = gsl_matrix_get(a, i, j);
1879 double mij = gsl_matrix_get(m, i, j);
1880 int foo = check(aij, mij, eps);
1881 if(foo) {
1882 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
1883 }
1884 s += foo;
1885 }
1886 }
1887 gsl_vector_free(w);
1888 gsl_vector_free(d);
1889 gsl_matrix_free(v);
1890 gsl_matrix_free(a);
1891 gsl_matrix_free(q);
1892 gsl_matrix_free(dqt);
1893 gsl_matrix_free (x);
1894
1895 return s;
1896 }
1897
test_SV_decomp_mod(void)1898 int test_SV_decomp_mod(void)
1899 {
1900 int f;
1901 int s = 0;
1902
1903 f = test_SV_decomp_mod_dim(m11, 2 * GSL_DBL_EPSILON);
1904 gsl_test(f, " SV_decomp_mod m(1,1)");
1905 s += f;
1906
1907 f = test_SV_decomp_mod_dim(m51, 2 * 64.0 * GSL_DBL_EPSILON);
1908 gsl_test(f, " SV_decomp_mod m(5,1)");
1909 s += f;
1910
1911 /* M<N not implemented yet */
1912 #if 0
1913 f = test_SV_decomp_mod_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
1914 gsl_test(f, " SV_decomp_mod m(3,5)");
1915 s += f;
1916 #endif
1917 f = test_SV_decomp_mod_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
1918 gsl_test(f, " SV_decomp_mod m(5,3)");
1919 s += f;
1920
1921 f = test_SV_decomp_mod_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
1922 gsl_test(f, " SV_decomp_mod moler(10)");
1923 s += f;
1924
1925 f = test_SV_decomp_mod_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
1926 gsl_test(f, " SV_decomp_mod hilbert(2)");
1927 s += f;
1928
1929 f = test_SV_decomp_mod_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
1930 gsl_test(f, " SV_decomp_mod hilbert(3)");
1931 s += f;
1932
1933 f = test_SV_decomp_mod_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1934 gsl_test(f, " SV_decomp_mod hilbert(4)");
1935 s += f;
1936
1937 f = test_SV_decomp_mod_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1938 gsl_test(f, " SV_decomp_mod hilbert(12)");
1939 s += f;
1940
1941 f = test_SV_decomp_mod_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1942 gsl_test(f, " SV_decomp_mod vander(2)");
1943 s += f;
1944
1945 f = test_SV_decomp_mod_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1946 gsl_test(f, " SV_decomp_mod vander(3)");
1947 s += f;
1948
1949 f = test_SV_decomp_mod_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1950 gsl_test(f, " SV_decomp_mod vander(4)");
1951 s += f;
1952
1953 f = test_SV_decomp_mod_dim(vander12, 1e-4);
1954 gsl_test(f, " SV_decomp_mod vander(12)");
1955 s += f;
1956
1957 f = test_SV_decomp_mod_dim(row3, 10 * GSL_DBL_EPSILON);
1958 gsl_test(f, " SV_decomp_mod row3");
1959 s += f;
1960
1961 f = test_SV_decomp_mod_dim(row5, 128 * GSL_DBL_EPSILON);
1962 gsl_test(f, " SV_decomp_mod row5");
1963 s += f;
1964
1965 f = test_SV_decomp_mod_dim(row12, 1024 * GSL_DBL_EPSILON);
1966 gsl_test(f, " SV_decomp_mod row12");
1967 s += f;
1968
1969 f = test_SV_decomp_mod_dim(inf5, 1024 * GSL_DBL_EPSILON);
1970 gsl_test(f, " SV_decomp_mod inf5");
1971 s += f;
1972
1973 f = test_SV_decomp_mod_dim(nan5, 1024 * GSL_DBL_EPSILON);
1974 gsl_test(f, " SV_decomp_mod nan5");
1975 s += f;
1976
1977
1978 {
1979 double i1, i2, i3, i4;
1980 double lower = -2, upper = 2;
1981
1982 for (i1 = lower; i1 <= upper; i1++)
1983 {
1984 for (i2 = lower; i2 <= upper; i2++)
1985 {
1986 for (i3 = lower; i3 <= upper; i3++)
1987 {
1988 for (i4 = lower; i4 <= upper; i4++)
1989 {
1990 gsl_matrix_set (A22, 0,0, i1);
1991 gsl_matrix_set (A22, 0,1, i2);
1992 gsl_matrix_set (A22, 1,0, i3);
1993 gsl_matrix_set (A22, 1,1, i4);
1994
1995 f = test_SV_decomp_mod_dim(A22, 16 * GSL_DBL_EPSILON);
1996 gsl_test(f, " SV_decomp_mod (2x2) A=[%g, %g; %g, %g]", i1,i2,i3,i4);
1997 s += f;
1998 }
1999 }
2000 }
2001 }
2002 }
2003
2004 {
2005 int i;
2006 double carry = 0, lower = 0, upper = 1;
2007 double *a = A33->data;
2008
2009 for (i=0; i<9; i++) {
2010 a[i] = lower;
2011 }
2012
2013 while (carry == 0.0) {
2014 f = test_SV_decomp_mod_dim(A33, 64 * GSL_DBL_EPSILON);
2015 gsl_test(f, " SV_decomp_mod (3x3) A=[ %g, %g, %g; %g, %g, %g; %g, %g, %g]",
2016 a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]);
2017
2018 /* increment */
2019 carry=1.0;
2020 for (i=9; carry > 0.0 && i>0 && i--;)
2021 {
2022 double v=a[i]+carry;
2023 carry = (v>upper) ? 1.0 : 0.0;
2024 a[i] = (v>upper) ? lower : v;
2025 }
2026 }
2027 }
2028
2029 #ifdef TEST_SVD_4X4
2030 {
2031 int i;
2032 double carry = 0, lower = 0, upper = 1;
2033 double *a = A44->data;
2034
2035 for (i=0; i<16; i++) {
2036 a[i] = lower;
2037 }
2038
2039 while (carry == 0.0) {
2040 f = test_SV_decomp_mod_dim(A44, 64 * GSL_DBL_EPSILON);
2041 gsl_test(f, " SV_decomp_mod (4x4) A=[ %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g]",
2042 a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9],
2043 a[10], a[11], a[12], a[13], a[14], a[15]);
2044
2045 /* increment */
2046 carry=1.0;
2047 for (i=16; carry>0.0 && i>0 && i--;)
2048 {
2049 double v=a[i]+carry;
2050 carry = (v>upper) ? 1.0 : 0.0;
2051 a[i] = (v>upper) ? lower : v;
2052 }
2053 }
2054 }
2055 #endif
2056
2057 return s;
2058 }
2059
2060
2061 int
test_SV_decomp_jacobi_dim(const gsl_matrix * m,double eps)2062 test_SV_decomp_jacobi_dim(const gsl_matrix * m, double eps)
2063 {
2064 int s = 0;
2065 double di1;
2066 unsigned long i,j, M = m->size1, N = m->size2;
2067
2068 gsl_matrix * v = gsl_matrix_alloc(M,N);
2069 gsl_matrix * a = gsl_matrix_alloc(M,N);
2070 gsl_matrix * q = gsl_matrix_alloc(N,N);
2071 gsl_matrix * dqt = gsl_matrix_alloc(N,N);
2072 gsl_vector * d = gsl_vector_alloc(N);
2073
2074 gsl_matrix_memcpy(v,m);
2075
2076 s += gsl_linalg_SV_decomp_jacobi(v, q, d);
2077 if (s)
2078 printf("call returned status = %d\n", s);
2079
2080 /* Check that singular values are non-negative and in non-decreasing
2081 order */
2082
2083 di1 = 0.0;
2084
2085 for (i = 0; i < N; i++)
2086 {
2087 double di = gsl_vector_get (d, i);
2088
2089 if (gsl_isnan (di))
2090 {
2091 continue; /* skip NaNs */
2092 }
2093
2094 if (di < 0) {
2095 s++;
2096 printf("singular value %lu = %22.18g < 0\n", i, di);
2097 }
2098
2099 if(i > 0 && di > di1) {
2100 s++;
2101 printf("singular value %lu = %22.18g vs previous %22.18g\n", i, di, di1);
2102 }
2103
2104 di1 = di;
2105 }
2106
2107 /* Scale dqt = D Q^T */
2108
2109 for (i = 0; i < N ; i++)
2110 {
2111 double di = gsl_vector_get (d, i);
2112
2113 for (j = 0; j < N; j++)
2114 {
2115 double qji = gsl_matrix_get(q, j, i);
2116 gsl_matrix_set (dqt, i, j, qji * di);
2117 }
2118 }
2119
2120 /* compute a = v dqt */
2121 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
2122
2123 for(i=0; i<M; i++) {
2124 for(j=0; j<N; j++) {
2125 double aij = gsl_matrix_get(a, i, j);
2126 double mij = gsl_matrix_get(m, i, j);
2127 int foo = check(aij, mij, eps);
2128 if(foo) {
2129 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
2130 }
2131 s += foo;
2132 }
2133 }
2134 gsl_vector_free(d);
2135 gsl_matrix_free(v);
2136 gsl_matrix_free(a);
2137 gsl_matrix_free(q);
2138 gsl_matrix_free(dqt);
2139
2140 return s;
2141 }
2142
test_SV_decomp_jacobi(void)2143 int test_SV_decomp_jacobi(void)
2144 {
2145 int f;
2146 int s = 0;
2147
2148 f = test_SV_decomp_jacobi_dim(m11, 2 * GSL_DBL_EPSILON);
2149 gsl_test(f, " SV_decomp_jacobi m(1,1)");
2150 s += f;
2151
2152 f = test_SV_decomp_jacobi_dim(m51, 2 * 64.0 * GSL_DBL_EPSILON);
2153 gsl_test(f, " SV_decomp_jacobi m(5,1)");
2154 s += f;
2155
2156 /* M<N not implemented yet */
2157 #if 0
2158 f = test_SV_decomp_jacobi_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
2159 gsl_test(f, " SV_decomp_jacobi m(3,5)");
2160 s += f;
2161 #endif
2162 f = test_SV_decomp_jacobi_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
2163 gsl_test(f, " SV_decomp_jacobi m(5,3)");
2164 s += f;
2165
2166 f = test_SV_decomp_jacobi_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
2167 gsl_test(f, " SV_decomp_jacobi moler(10)");
2168 s += f;
2169
2170 f = test_SV_decomp_jacobi_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
2171 gsl_test(f, " SV_decomp_jacobi hilbert(2)");
2172 s += f;
2173
2174 f = test_SV_decomp_jacobi_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
2175 gsl_test(f, " SV_decomp_jacobi hilbert(3)");
2176 s += f;
2177
2178 f = test_SV_decomp_jacobi_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
2179 gsl_test(f, " SV_decomp_jacobi hilbert(4)");
2180 s += f;
2181
2182 f = test_SV_decomp_jacobi_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
2183 gsl_test(f, " SV_decomp_jacobi hilbert(12)");
2184 s += f;
2185
2186 f = test_SV_decomp_jacobi_dim(vander2, 8.0 * GSL_DBL_EPSILON);
2187 gsl_test(f, " SV_decomp_jacobi vander(2)");
2188 s += f;
2189
2190 f = test_SV_decomp_jacobi_dim(vander3, 64.0 * GSL_DBL_EPSILON);
2191 gsl_test(f, " SV_decomp_jacobi vander(3)");
2192 s += f;
2193
2194 f = test_SV_decomp_jacobi_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
2195 gsl_test(f, " SV_decomp_jacobi vander(4)");
2196 s += f;
2197
2198 f = test_SV_decomp_jacobi_dim(vander12, 1e-4);
2199 gsl_test(f, " SV_decomp_jacobi vander(12)");
2200 s += f;
2201
2202 f = test_SV_decomp_jacobi_dim(row3, 10 * GSL_DBL_EPSILON);
2203 gsl_test(f, " SV_decomp_jacobi row3");
2204 s += f;
2205
2206 f = test_SV_decomp_jacobi_dim(row5, 128 * GSL_DBL_EPSILON);
2207 gsl_test(f, " SV_decomp_jacobi row5");
2208 s += f;
2209
2210 f = test_SV_decomp_jacobi_dim(row12, 1024 * GSL_DBL_EPSILON);
2211 gsl_test(f, " SV_decomp_jacobi row12");
2212 s += f;
2213
2214
2215 #ifdef TEST_JACOBI_INF
2216 f = test_SV_decomp_jacobi_dim(inf5, 1024 * GSL_DBL_EPSILON);
2217 gsl_test(f, " SV_decomp_jacobi inf5");
2218 s += f;
2219
2220 f = test_SV_decomp_jacobi_dim(nan5, 1024 * GSL_DBL_EPSILON);
2221 gsl_test(f, " SV_decomp_jacobi nan5");
2222 s += f;
2223 #endif
2224
2225 {
2226 double i1, i2, i3, i4;
2227 double lower = -2, upper = 2;
2228
2229 for (i1 = lower; i1 <= upper; i1++)
2230 {
2231 for (i2 = lower; i2 <= upper; i2++)
2232 {
2233 for (i3 = lower; i3 <= upper; i3++)
2234 {
2235 for (i4 = lower; i4 <= upper; i4++)
2236 {
2237 gsl_matrix_set (A22, 0,0, i1);
2238 gsl_matrix_set (A22, 0,1, i2);
2239 gsl_matrix_set (A22, 1,0, i3);
2240 gsl_matrix_set (A22, 1,1, i4);
2241
2242 f = test_SV_decomp_jacobi_dim(A22, 16 * GSL_DBL_EPSILON);
2243 gsl_test(f, " SV_decomp_jacobi (2x2) A=[%g, %g; %g, %g]", i1,i2,i3,i4);
2244 s += f;
2245 }
2246 }
2247 }
2248 }
2249 }
2250
2251 {
2252 int i;
2253 double carry = 0, lower = 0, upper = 1;
2254 double *a = A33->data;
2255
2256 for (i=0; i<9; i++) {
2257 a[i] = lower;
2258 }
2259
2260 while (carry == 0.0) {
2261 f = test_SV_decomp_jacobi_dim(A33, 64 * GSL_DBL_EPSILON);
2262 gsl_test(f, " SV_decomp_jacobi (3x3) A=[ %g, %g, %g; %g, %g, %g; %g, %g, %g]",
2263 a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]);
2264
2265 /* increment */
2266 carry=1.0;
2267 for (i=9; carry > 0.0 && i>0 && i--;)
2268 {
2269 double v=a[i]+carry;
2270 carry = (v>upper) ? 1.0 : 0.0;
2271 a[i] = (v>upper) ? lower : v;
2272 }
2273 }
2274 }
2275
2276 #ifdef TEST_SVD_4X4
2277 {
2278 int i;
2279 unsigned long k = 0;
2280 double carry = 0, lower = 0, upper = 1;
2281 double *a = A44->data;
2282
2283 for (i=0; i<16; i++) {
2284 a[i] = lower;
2285 }
2286
2287 while (carry == 0.0) {
2288 k++;
2289 f = test_SV_decomp_jacobi_dim(A44, 64 * GSL_DBL_EPSILON);
2290 gsl_test(f, " SV_decomp_jacobi (4x4) A=[ %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g] %lu",
2291 a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9],
2292 a[10], a[11], a[12], a[13], a[14], a[15], k);
2293 /* increment */
2294 carry=1.0;
2295 for (i=16; carry > 0.0 && i>0 && i--;)
2296 {
2297 double v=a[i]+carry;
2298 carry = (v>upper) ? 1.0 : 0.0;
2299 a[i] = (v>upper) ? lower : v;
2300 }
2301 }
2302 }
2303 #endif
2304
2305 {
2306 int i;
2307 unsigned long k = 0;
2308 double carry = 0, lower = 0, upper = 1;
2309 double *a = A55->data;
2310
2311 for (i=0; i<25; i++) {
2312 a[i] = lower;
2313 }
2314
2315 while (carry == 0.0) {
2316 k++;
2317
2318 if (k % 1001 == 0)
2319 {
2320 f = test_SV_decomp_jacobi_dim(A55, 64 * GSL_DBL_EPSILON);
2321 gsl_test(f, " SV_decomp_jacobi (5x5) case=%lu",k);
2322 }
2323
2324 /* increment */
2325 carry=1.0;
2326 for (i=25; carry >0.0 && i>0 && i--;)
2327 {
2328 double v=a[i]+carry;
2329 carry = (v>upper) ? 1.0 : 0.0;
2330 a[i] = (v>upper) ? lower : v;
2331 }
2332 }
2333 }
2334
2335
2336 return s;
2337 }
2338
2339
2340 int
test_cholesky_solve_dim(const gsl_matrix * m,const double * actual,double eps)2341 test_cholesky_solve_dim(const gsl_matrix * m, const double * actual, double eps)
2342 {
2343 int s = 0;
2344 unsigned long i, dim = m->size1;
2345
2346 gsl_vector * rhs = gsl_vector_alloc(dim);
2347 gsl_matrix * u = gsl_matrix_alloc(dim,dim);
2348 gsl_vector * x = gsl_vector_calloc(dim);
2349 gsl_matrix_memcpy(u,m);
2350 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
2351 s += gsl_linalg_cholesky_decomp1(u);
2352 s += gsl_linalg_cholesky_solve(u, rhs, x);
2353 for(i=0; i<dim; i++) {
2354 int foo = check(gsl_vector_get(x, i), actual[i], eps);
2355 if(foo) {
2356 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2357 }
2358 s += foo;
2359 }
2360 gsl_vector_free(x);
2361 gsl_matrix_free(u);
2362 gsl_vector_free(rhs);
2363
2364 return s;
2365 }
2366
2367 int
test_cholesky_solve2_dim(const gsl_matrix * m,const double * actual,double eps)2368 test_cholesky_solve2_dim(const gsl_matrix * m, const double * actual, double eps)
2369 {
2370 int s = 0;
2371 unsigned long i, dim = m->size1;
2372
2373 gsl_vector * rhs = gsl_vector_alloc(dim);
2374 gsl_matrix * u = gsl_matrix_alloc(dim,dim);
2375 gsl_vector * x = gsl_vector_calloc(dim);
2376 gsl_vector * D = gsl_vector_calloc(dim);
2377 gsl_matrix_memcpy(u,m);
2378 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
2379 s += gsl_linalg_cholesky_decomp2(u, D);
2380 s += gsl_linalg_cholesky_solve2(u, D, rhs, x);
2381 for(i=0; i<dim; i++) {
2382 int foo = check(gsl_vector_get(x, i), actual[i], eps);
2383 if(foo) {
2384 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2385 }
2386 s += foo;
2387 }
2388 gsl_vector_free(x);
2389 gsl_matrix_free(u);
2390 gsl_vector_free(rhs);
2391 gsl_vector_free(D);
2392
2393 return s;
2394 }
2395
2396 int
test_cholesky_solve(void)2397 test_cholesky_solve(void)
2398 {
2399 int f;
2400 int s = 0;
2401
2402 f = test_cholesky_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
2403 gsl_test(f, " cholesky_solve hilbert(2)");
2404 s += f;
2405
2406 f = test_cholesky_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
2407 gsl_test(f, " cholesky_solve hilbert(3)");
2408 s += f;
2409
2410 f = test_cholesky_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
2411 gsl_test(f, " cholesky_solve hilbert(4)");
2412 s += f;
2413
2414 f = test_cholesky_solve_dim(hilb12, hilb12_solution, 0.5);
2415 gsl_test(f, " cholesky_solve hilbert(12)");
2416 s += f;
2417
2418 /* test scaled Cholesky routines */
2419
2420 f = test_cholesky_solve2_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
2421 gsl_test(f, " cholesky_solve2 hilbert(2)");
2422 s += f;
2423
2424 f = test_cholesky_solve2_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
2425 gsl_test(f, " cholesky_solve2 hilbert(3)");
2426 s += f;
2427
2428 f = test_cholesky_solve2_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
2429 gsl_test(f, " cholesky_solve2 hilbert(4)");
2430 s += f;
2431
2432 f = test_cholesky_solve2_dim(hilb12, hilb12_solution, 0.5);
2433 gsl_test(f, " cholesky_solve2 hilbert(12)");
2434 s += f;
2435
2436 return s;
2437 }
2438
2439 int
test_cholesky_decomp_unit_dim(const gsl_matrix * m,double eps)2440 test_cholesky_decomp_unit_dim(const gsl_matrix * m, double eps)
2441 {
2442 int s = 0;
2443 const unsigned long M = m->size1;
2444 const unsigned long N = m->size2;
2445 unsigned long i,j;
2446
2447 gsl_matrix * v = gsl_matrix_alloc(M,N);
2448 gsl_matrix * a = gsl_matrix_alloc(M,N);
2449 gsl_matrix * l = gsl_matrix_alloc(M,N);
2450 gsl_matrix * lt = gsl_matrix_alloc(N,N);
2451 gsl_matrix * dm = gsl_matrix_alloc(M,N);
2452 gsl_vector * dv = gsl_vector_alloc(M);
2453
2454 gsl_matrix_memcpy(v,m);
2455
2456 s += gsl_linalg_cholesky_decomp_unit(v, dv);
2457
2458 /*
2459 for(i = 0; i < M; i++)
2460 {
2461 for(j = 0; j < N; j++)
2462 {
2463 printf("v[%lu,%lu]: %22.18e\n", i,j, gsl_matrix_get(v, i, j));
2464 }
2465 }
2466
2467
2468 for(i = 0; i < M; i++)
2469 {
2470 printf("d[%lu]: %22.18e\n", i, gsl_vector_get(dv, i));
2471 }
2472 */
2473
2474 /* put L and transpose(L) into separate matrices */
2475
2476 for(i = 0; i < N ; i++)
2477 {
2478 for(j = 0; j < N; j++)
2479 {
2480 const double vij = gsl_matrix_get(v, i, j);
2481 gsl_matrix_set (l, i, j, i>=j ? vij : 0);
2482 gsl_matrix_set (lt, i, j, i<=j ? vij : 0);
2483 }
2484 }
2485
2486 /* put D into its own matrix */
2487
2488 gsl_matrix_set_zero(dm);
2489 for(i = 0; i < M; ++i) gsl_matrix_set(dm, i, i, gsl_vector_get(dv, i));
2490
2491 /* compute a = L * D * transpose(L); uses v for temp space */
2492
2493 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, dm, lt, 0.0, v);
2494 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, v, 0.0, a);
2495
2496 for(i = 0; i < M; i++)
2497 {
2498 for(j = 0; j < N; j++)
2499 {
2500 const double aij = gsl_matrix_get(a, i, j);
2501 const double mij = gsl_matrix_get(m, i, j);
2502 int foo = check(aij, mij, eps);
2503 if(foo)
2504 {
2505 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
2506 }
2507 s += foo;
2508 }
2509 }
2510
2511 gsl_vector_free(dv);
2512 gsl_matrix_free(dm);
2513 gsl_matrix_free(lt);
2514 gsl_matrix_free(l);
2515 gsl_matrix_free(v);
2516 gsl_matrix_free(a);
2517
2518 return s;
2519 }
2520
test_cholesky_decomp_unit(void)2521 int test_cholesky_decomp_unit(void)
2522 {
2523 int f;
2524 int s = 0;
2525
2526 f = test_cholesky_decomp_unit_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
2527 gsl_test(f, " cholesky_decomp_unit hilbert(2)");
2528 s += f;
2529
2530 f = test_cholesky_decomp_unit_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
2531 gsl_test(f, " cholesky_decomp_unit hilbert(3)");
2532 s += f;
2533
2534 f = test_cholesky_decomp_unit_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
2535 gsl_test(f, " cholesky_decomp_unit hilbert(4)");
2536 s += f;
2537
2538 f = test_cholesky_decomp_unit_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
2539 gsl_test(f, " cholesky_decomp_unit hilbert(12)");
2540 s += f;
2541
2542 return s;
2543 }
2544
2545 int
test_choleskyc_solve_dim(const gsl_matrix_complex * m,const gsl_vector_complex * actual,double eps)2546 test_choleskyc_solve_dim(const gsl_matrix_complex * m, const gsl_vector_complex * actual, double eps)
2547 {
2548 int s = 0;
2549 unsigned long i, dim = m->size1;
2550 gsl_complex z;
2551 gsl_vector_complex * rhs = gsl_vector_complex_alloc(dim);
2552 gsl_matrix_complex * u = gsl_matrix_complex_alloc(dim,dim);
2553 gsl_vector_complex * x = gsl_vector_complex_calloc(dim);
2554 GSL_SET_IMAG(&z, 0.0);
2555 gsl_matrix_complex_memcpy(u,m);
2556 for(i=0; i<dim; i++)
2557 {
2558 GSL_SET_REAL(&z, i + 1.0);
2559 gsl_vector_complex_set(rhs, i, z);
2560 }
2561 s += gsl_linalg_complex_cholesky_decomp(u);
2562 s += gsl_linalg_complex_cholesky_solve(u, rhs, x);
2563 for(i=0; i<dim; i++) {
2564 gsl_complex y = gsl_vector_complex_get(x, i);
2565 gsl_complex a = gsl_vector_complex_get(actual, i);
2566 int foo = check(GSL_REAL(y), GSL_REAL(a), eps);
2567 int foo2 = check(GSL_IMAG(y), GSL_IMAG(a), eps);
2568 if(foo || foo2) {
2569 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, GSL_REAL(y), GSL_REAL(a));
2570 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, GSL_IMAG(y), GSL_IMAG(a));
2571 }
2572 s += foo + foo2;
2573 }
2574 gsl_vector_complex_free(x);
2575 gsl_matrix_complex_free(u);
2576 gsl_vector_complex_free(rhs);
2577
2578 return s;
2579 } /* test_choleskyc_solve_dim() */
2580
2581 int
test_choleskyc_solve(void)2582 test_choleskyc_solve(void)
2583 {
2584 double data7[] = { 66,0, 0,64, 126,63, 124,-62, 61,-61, 60,60, 0,-59,
2585 0,-64, 65,0, 62,-124, -61,-122, -60,-60, 59,-59, -58,0,
2586 126,-63, 62,124, 308,0, 180,-240, 59,-177, 174,58, -57,-114,
2587 124,62, -61,122, 180,240, 299,0, 174,-58, 57,171, 56,-112,
2588 61,61, -60,60, 59,177, 174,58, 119,0, 0,112, 55,-55,
2589 60,-60, 59,59, 174,-58, 57,-171, 0,-112, 116,0, -54,-54,
2590 0,59, -58,0, -57,114, 56,112, 55,55, -54,54, 60,0 };
2591 double data7_sol[] = { -0.524944196428570,0.209123883928571,
2592 1.052873883928572,0.712444196428571,
2593 0.117568824404762,0.443191964285714,
2594 0.412862723214286,-0.356696428571429,
2595 0.815931919642858,-0.265820312500000,
2596 0.777929687500000,0.119484747023810,
2597 1.058733258928571,-0.132087053571429 };
2598 gsl_matrix_complex_view cp7 = gsl_matrix_complex_view_array(data7, 7, 7);
2599 gsl_vector_complex_view cp7_sol = gsl_vector_complex_view_array(data7_sol, 7);
2600 int f;
2601 int s = 0;
2602
2603 f = test_choleskyc_solve_dim(&cp7.matrix, &cp7_sol.vector, 1024.0 * 1024.0 * GSL_DBL_EPSILON);
2604 gsl_test(f, " complex_cholesky_solve complex(7)");
2605 s += f;
2606
2607 return s;
2608 }
2609
2610 int
test_HH_solve_dim(const gsl_matrix * m,const double * actual,double eps)2611 test_HH_solve_dim(const gsl_matrix * m, const double * actual, double eps)
2612 {
2613 int s = 0;
2614 unsigned long i, dim = m->size1;
2615
2616 gsl_permutation * perm = gsl_permutation_alloc(dim);
2617 gsl_matrix * hh = gsl_matrix_alloc(dim,dim);
2618 gsl_vector * d = gsl_vector_alloc(dim);
2619 gsl_vector * x = gsl_vector_alloc(dim);
2620 gsl_matrix_memcpy(hh,m);
2621 for(i=0; i<dim; i++) gsl_vector_set(x, i, i+1.0);
2622 s += gsl_linalg_HH_svx(hh, x);
2623 for(i=0; i<dim; i++) {
2624 int foo = check(gsl_vector_get(x, i),actual[i],eps);
2625 if( foo) {
2626 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2627 }
2628 s += foo;
2629 }
2630 gsl_vector_free(x);
2631 gsl_vector_free(d);
2632 gsl_matrix_free(hh);
2633 gsl_permutation_free(perm);
2634
2635 return s;
2636 }
2637
test_HH_solve(void)2638 int test_HH_solve(void)
2639 {
2640 int f;
2641 int s = 0;
2642
2643 f = test_HH_solve_dim(hilb2, hilb2_solution, 8.0 * GSL_DBL_EPSILON);
2644 gsl_test(f, " HH_solve hilbert(2)");
2645 s += f;
2646
2647 f = test_HH_solve_dim(hilb3, hilb3_solution, 128.0 * GSL_DBL_EPSILON);
2648 gsl_test(f, " HH_solve hilbert(3)");
2649 s += f;
2650
2651 f = test_HH_solve_dim(hilb4, hilb4_solution, 2.0 * 1024.0 * GSL_DBL_EPSILON);
2652 gsl_test(f, " HH_solve hilbert(4)");
2653 s += f;
2654
2655 f = test_HH_solve_dim(hilb12, hilb12_solution, 0.5);
2656 gsl_test(f, " HH_solve hilbert(12)");
2657 s += f;
2658
2659 f = test_HH_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
2660 gsl_test(f, " HH_solve vander(2)");
2661 s += f;
2662
2663 f = test_HH_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
2664 gsl_test(f, " HH_solve vander(3)");
2665 s += f;
2666
2667 f = test_HH_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
2668 gsl_test(f, " HH_solve vander(4)");
2669 s += f;
2670
2671 f = test_HH_solve_dim(vander12, vander12_solution, 0.05);
2672 gsl_test(f, " HH_solve vander(12)");
2673 s += f;
2674
2675 return s;
2676 }
2677
2678
2679 int
test_TDS_solve_dim(unsigned long dim,double d,double od,const double * actual,double eps)2680 test_TDS_solve_dim(unsigned long dim, double d, double od, const double * actual, double eps)
2681 {
2682 int s = 0;
2683 unsigned long i;
2684
2685 gsl_vector * offdiag = vector_alloc(dim-1);
2686 gsl_vector * diag = vector_alloc(dim);
2687 gsl_vector * rhs = vector_alloc(dim);
2688 gsl_vector * x = vector_alloc(dim);
2689
2690 for(i=0; i<dim; i++) {
2691 gsl_vector_set(diag, i, d);
2692 gsl_vector_set(rhs, i, i + 1.0);
2693 }
2694 for(i=0; i<dim-1; i++) {
2695 gsl_vector_set(offdiag, i, od);
2696 }
2697
2698 s += gsl_linalg_solve_symm_tridiag(diag, offdiag, rhs, x);
2699
2700 for(i=0; i<dim; i++) {
2701 double si = gsl_vector_get(x, i);
2702 double ai = actual[i];
2703 int foo = check(si, ai, eps);
2704 if(foo) {
2705 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2706 }
2707 s += foo;
2708 }
2709
2710 vector_free(x);
2711 vector_free(rhs);
2712 vector_free(diag);
2713 vector_free(offdiag);
2714
2715 return s;
2716 }
2717
2718
test_TDS_solve(void)2719 int test_TDS_solve(void)
2720 {
2721 int f;
2722 int s = 0;
2723
2724 {
2725 double actual[] = {0.0, 2.0};
2726 f = test_TDS_solve_dim(2, 1.0, 0.5, actual, 8.0 * GSL_DBL_EPSILON);
2727 gsl_test(f, " solve_TDS dim=2 A");
2728 s += f;
2729 }
2730
2731 {
2732 double actual[] = {3.0/8.0, 15.0/8.0};
2733 f = test_TDS_solve_dim(2, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
2734 gsl_test(f, " solve_TDS dim=2 B");
2735 s += f;
2736 }
2737
2738 {
2739 double actual[] = {5.0/8.0, 9.0/8.0, 2.0, 15.0/8.0, 35.0/8.0};
2740 f = test_TDS_solve_dim(5, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
2741 gsl_test(f, " solve_TDS dim=5");
2742 s += f;
2743 }
2744
2745 return s;
2746 }
2747
2748 int
test_TDS_cyc_solve_one(const unsigned long dim,const double * d,const double * od,const double * r,const double * actual,double eps)2749 test_TDS_cyc_solve_one(const unsigned long dim,
2750 const double * d, const double * od,
2751 const double * r, const double * actual, double eps)
2752 {
2753 int s = 0;
2754 unsigned long i;
2755
2756 gsl_vector * offdiag = vector_alloc(dim);
2757 gsl_vector * diag = vector_alloc(dim);
2758 gsl_vector * rhs = vector_alloc(dim);
2759 gsl_vector * x = vector_alloc(dim);
2760
2761 for(i=0; i<dim; i++) {
2762 gsl_vector_set(diag, i, d[i]);
2763 gsl_vector_set(rhs, i, r[i]);
2764 gsl_vector_set(offdiag, i, od[i]);
2765 }
2766
2767 s += gsl_linalg_solve_symm_cyc_tridiag(diag, offdiag, rhs, x);
2768
2769 for(i=0; i<dim; i++) {
2770 double si = gsl_vector_get(x, i);
2771 double ai = actual[i];
2772 int foo = check(si, ai, eps);
2773 if(foo) {
2774 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2775 }
2776 s += foo;
2777 }
2778
2779 vector_free(x);
2780 vector_free(rhs);
2781 vector_free(diag);
2782 vector_free(offdiag);
2783
2784 return s;
2785 }
2786
test_TDS_cyc_solve(void)2787 int test_TDS_cyc_solve(void)
2788 {
2789 int f;
2790 int s = 0;
2791
2792 #ifdef SUPPORT_UNDERSIZE_CYC
2793 {
2794 unsigned long dim = 1;
2795 double diag[] = { 2 };
2796 double offdiag[] = { 3 };
2797 double rhs[] = { 7 };
2798 double actual[] = { 3.5 };
2799
2800 f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 28.0 * GSL_DBL_EPSILON);
2801 gsl_test(f, " solve_TDS_cyc dim=%lu A", dim);
2802 s += f;
2803 }
2804
2805 {
2806 unsigned long dim = 2;
2807 double diag[] = { 1, 2 };
2808 double offdiag[] = { 3, 4 };
2809 double rhs[] = { 7, -7 };
2810 double actual[] = { -5, 4 };
2811
2812 f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 28.0 * GSL_DBL_EPSILON);
2813 gsl_test(f, " solve_TDS_cyc dim=%lu A", dim);
2814 s += f;
2815 }
2816 #endif
2817
2818 {
2819 unsigned long dim = 3;
2820 double diag[] = { 1, 1, 1 };
2821 double offdiag[] = { 3, 3, 3 };
2822 double rhs[] = { 7, -7, 7 };
2823 double actual[] = { -2, 5, -2 };
2824
2825 f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 28.0 * GSL_DBL_EPSILON);
2826 gsl_test(f, " solve_TDS_cyc dim=%lu A", dim);
2827 s += f;
2828 }
2829
2830 {
2831 unsigned long dim = 5;
2832 double diag[] = { 4, 2, 1, 2, 4 };
2833 double offdiag[] = { 1, 1, 1, 1, 1 };
2834 double rhs[] = { 30, -24, 3, 21, -30 };
2835 double actual[] = { 12, 3, -42, 42, -21 };
2836
2837 /* f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 7.0 * GSL_DBL_EPSILON);
2838 FIXME: bad accuracy */
2839 f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 40.0 * GSL_DBL_EPSILON);
2840 gsl_test(f, " solve_TDS_cyc dim=%lu B", dim);
2841 s += f;
2842 }
2843
2844 return s;
2845 }
2846
2847 int
test_TDN_solve_dim(unsigned long dim,double d,double a,double b,const double * actual,double eps)2848 test_TDN_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps)
2849 {
2850 int s = 0;
2851 unsigned long i;
2852
2853 gsl_vector * abovediag = vector_alloc(dim-1);
2854 gsl_vector * belowdiag = vector_alloc(dim-1);
2855 gsl_vector * diag = vector_alloc(dim);
2856 gsl_vector * rhs = vector_alloc(dim);
2857 gsl_vector * x = vector_alloc(dim);
2858
2859 for(i=0; i<dim; i++) {
2860 gsl_vector_set(diag, i, d);
2861 gsl_vector_set(rhs, i, i + 1.0);
2862 }
2863 for(i=0; i<dim-1; i++) {
2864 gsl_vector_set(abovediag, i, a);
2865 gsl_vector_set(belowdiag, i, b);
2866 }
2867
2868 s += gsl_linalg_solve_tridiag(diag, abovediag, belowdiag, rhs, x);
2869
2870 for(i=0; i<dim; i++) {
2871 double si = gsl_vector_get(x, i);
2872 double ai = actual[i];
2873 int foo = check(si, ai, eps);
2874 if(foo) {
2875 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2876 }
2877 s += foo;
2878 }
2879
2880 vector_free(x);
2881 vector_free(rhs);
2882 vector_free(diag);
2883 vector_free(abovediag);
2884 vector_free(belowdiag);
2885
2886 return s;
2887 }
2888
2889
test_TDN_solve(void)2890 int test_TDN_solve(void)
2891 {
2892 int f;
2893 int s = 0;
2894 double actual[16];
2895
2896 actual[0] = -7.0/3.0;
2897 actual[1] = 5.0/3.0;
2898 actual[2] = 4.0/3.0;
2899 f = test_TDN_solve_dim(3, 1.0, 2.0, 1.0, actual, 2.0 * GSL_DBL_EPSILON);
2900 gsl_test(f, " solve_TDN dim=2 A");
2901 s += f;
2902
2903 actual[0] = 0.75;
2904 actual[1] = 0.75;
2905 actual[2] = 2.625;
2906
2907 f = test_TDN_solve_dim(3, 1.0, 1.0/3.0, 1.0/2.0, actual, 2.0 * GSL_DBL_EPSILON);
2908 gsl_test(f, " solve_TDN dim=2 B");
2909 s += f;
2910
2911 actual[0] = 99.0/140.0;
2912 actual[1] = 41.0/35.0;
2913 actual[2] = 19.0/10.0;
2914 actual[3] = 72.0/35.0;
2915 actual[4] = 139.0/35.0;
2916 f = test_TDN_solve_dim(5, 1.0, 1.0/4.0, 1.0/2.0, actual, 35.0/8.0 * GSL_DBL_EPSILON);
2917 gsl_test(f, " solve_TDN dim=5");
2918 s += f;
2919
2920 return s;
2921 }
2922
2923 int
test_TDN_cyc_solve_dim(unsigned long dim,double d,double a,double b,const double * actual,double eps)2924 test_TDN_cyc_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps)
2925 {
2926 int s = 0;
2927 unsigned long i;
2928
2929 gsl_vector * abovediag = vector_alloc(dim);
2930 gsl_vector * belowdiag = vector_alloc(dim);
2931 gsl_vector * diag = vector_alloc(dim);
2932 gsl_vector * rhs = vector_alloc(dim);
2933 gsl_vector * x = vector_alloc(dim);
2934
2935 for(i=0; i<dim; i++) {
2936 gsl_vector_set(diag, i, d);
2937 gsl_vector_set(rhs, i, i + 1.0);
2938 }
2939 for(i=0; i<dim; i++) {
2940 gsl_vector_set(abovediag, i, a);
2941 gsl_vector_set(belowdiag, i, b);
2942 }
2943
2944 s += gsl_linalg_solve_cyc_tridiag(diag, abovediag, belowdiag, rhs, x);
2945
2946 for(i=0; i<dim; i++) {
2947 double si = gsl_vector_get(x, i);
2948 double ai = actual[i];
2949 int foo = check(si, ai, eps);
2950 if(foo) {
2951 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2952 }
2953 s += foo;
2954 }
2955
2956 vector_free(x);
2957 vector_free(rhs);
2958 vector_free(diag);
2959 vector_free(abovediag);
2960 vector_free(belowdiag);
2961
2962 return s;
2963 }
2964
2965
test_TDN_cyc_solve(void)2966 int test_TDN_cyc_solve(void)
2967 {
2968 int f;
2969 int s = 0;
2970 double actual[16];
2971
2972 actual[0] = 3.0/2.0;
2973 actual[1] = -1.0/2.0;
2974 actual[2] = 1.0/2.0;
2975 f = test_TDN_cyc_solve_dim(3, 1.0, 2.0, 1.0, actual, 32.0 * GSL_DBL_EPSILON);
2976 gsl_test(f, " solve_TDN_cyc dim=2 A");
2977 s += f;
2978
2979 actual[0] = -5.0/22.0;
2980 actual[1] = -3.0/22.0;
2981 actual[2] = 29.0/22.0;
2982 actual[3] = -9.0/22.0;
2983 actual[4] = 43.0/22.0;
2984
2985 f = test_TDN_cyc_solve_dim(5, 3.0, 2.0, 1.0, actual, 66.0 * GSL_DBL_EPSILON);
2986 gsl_test(f, " solve_TDN_cyc dim=5");
2987 s += f;
2988
2989 return s;
2990 }
2991
2992 int
test_bidiag_decomp_dim(const gsl_matrix * m,double eps)2993 test_bidiag_decomp_dim(const gsl_matrix * m, double eps)
2994 {
2995 int s = 0;
2996 unsigned long i,j,k,r, M = m->size1, N = m->size2;
2997
2998 gsl_matrix * A = gsl_matrix_alloc(M,N);
2999 gsl_matrix * a = gsl_matrix_alloc(M,N);
3000 gsl_matrix * b = gsl_matrix_alloc(N,N);
3001
3002 gsl_matrix * u = gsl_matrix_alloc(M,N);
3003 gsl_matrix * v = gsl_matrix_alloc(N,N);
3004
3005 gsl_vector * tau1 = gsl_vector_alloc(N);
3006 gsl_vector * tau2 = gsl_vector_alloc(N-1);
3007 gsl_vector * d = gsl_vector_alloc(N);
3008 gsl_vector * sd = gsl_vector_alloc(N-1);
3009
3010 gsl_matrix_memcpy(A,m);
3011
3012 s += gsl_linalg_bidiag_decomp(A, tau1, tau2);
3013 s += gsl_linalg_bidiag_unpack(A, tau1, u, tau2, v, d, sd);
3014
3015 gsl_matrix_set_zero(b);
3016 for (i = 0; i < N; i++) gsl_matrix_set(b, i,i, gsl_vector_get(d,i));
3017 for (i = 0; i < N-1; i++) gsl_matrix_set(b, i,i+1, gsl_vector_get(sd,i));
3018
3019 /* Compute A = U B V^T */
3020
3021 for (i = 0; i < M ; i++)
3022 {
3023 for (j = 0; j < N; j++)
3024 {
3025 double sum = 0;
3026
3027 for (k = 0; k < N; k++)
3028 {
3029 for (r = 0; r < N; r++)
3030 {
3031 sum += gsl_matrix_get(u, i, k) * gsl_matrix_get (b, k, r)
3032 * gsl_matrix_get(v, j, r);
3033 }
3034 }
3035 gsl_matrix_set (a, i, j, sum);
3036 }
3037 }
3038
3039 for(i=0; i<M; i++) {
3040 for(j=0; j<N; j++) {
3041 double aij = gsl_matrix_get(a, i, j);
3042 double mij = gsl_matrix_get(m, i, j);
3043 int foo = check(aij, mij, eps);
3044 if(foo) {
3045 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
3046 }
3047 s += foo;
3048 }
3049 }
3050
3051 gsl_matrix_free(A);
3052 gsl_matrix_free(a);
3053 gsl_matrix_free(u);
3054 gsl_matrix_free(v);
3055 gsl_matrix_free(b);
3056 gsl_vector_free(tau1);
3057 gsl_vector_free(tau2);
3058 gsl_vector_free(d);
3059 gsl_vector_free(sd);
3060
3061 return s;
3062 }
3063
test_bidiag_decomp(void)3064 int test_bidiag_decomp(void)
3065 {
3066 int f;
3067 int s = 0;
3068
3069 f = test_bidiag_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
3070 gsl_test(f, " bidiag_decomp m(5,3)");
3071 s += f;
3072
3073 f = test_bidiag_decomp_dim(m97, 2 * 64.0 * GSL_DBL_EPSILON);
3074 gsl_test(f, " bidiag_decomp m(9,7)");
3075 s += f;
3076
3077 f = test_bidiag_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
3078 gsl_test(f, " bidiag_decomp hilbert(2)");
3079 s += f;
3080
3081 f = test_bidiag_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
3082 gsl_test(f, " bidiag_decomp hilbert(3)");
3083 s += f;
3084
3085 f = test_bidiag_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
3086 gsl_test(f, " bidiag_decomp hilbert(4)");
3087 s += f;
3088
3089 f = test_bidiag_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
3090 gsl_test(f, " bidiag_decomp hilbert(12)");
3091 s += f;
3092
3093 return s;
3094 }
3095
3096 int
test_tri_invert2(CBLAS_UPLO_t Uplo,CBLAS_DIAG_t Diag,gsl_rng * r,const double tol)3097 test_tri_invert2(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_rng * r, const double tol)
3098 {
3099 const size_t N_max = 200;
3100 int s = 0;
3101 size_t n, i, j;
3102
3103 for (n = 1; n <= N_max; ++n)
3104 {
3105 gsl_matrix *T = gsl_matrix_alloc(n, n);
3106 gsl_matrix *B = gsl_matrix_alloc(n, n);
3107
3108 /* generate random triangular matrix */
3109 create_tri_matrix(Uplo, Diag, T, r);
3110
3111 /* compute B = T^{-1} */
3112 gsl_matrix_memcpy(B, T);
3113 gsl_linalg_tri_invert(Uplo, Diag, B);
3114
3115 /* compute B = T * T^{-1} */
3116 gsl_blas_dtrmm(CblasLeft, Uplo, CblasNoTrans, Diag, 1.0, T, B);
3117
3118 /* test B = I */
3119 if (Uplo == CblasUpper)
3120 {
3121 for (i = 0; i < n; ++i)
3122 {
3123 double Bii = gsl_matrix_get(B, i, i);
3124
3125 gsl_test_abs(Bii, 1.0, tol, "tri_invert[%zu,%zu] N=%zu upper %s",
3126 i, i, n,
3127 (Diag == CblasNonUnit) ? "NonUnit" : "Unit");
3128
3129 for (j = i + 1; j < n; ++j)
3130 {
3131 double Bij = gsl_matrix_get(B, i, j);
3132
3133 gsl_test_abs(Bij, 0.0, tol, "tri_invert[%zu,%zu] N=%zu upper %s",
3134 i, j, n,
3135 (Diag == CblasNonUnit) ? "NonUnit" : "Unit");
3136 }
3137 }
3138 }
3139 else
3140 {
3141 for (i = 0; i < n; ++i)
3142 {
3143 double Bii = gsl_matrix_get(B, i, i);
3144
3145 gsl_test_abs(Bii, 1.0, tol, "tri_invert[%zu,%zu] N=%zu lower %s",
3146 i, i, n,
3147 (Diag == CblasNonUnit) ? "NonUnit" : "Unit");
3148
3149 for (j = 0; j < i; ++j)
3150 {
3151 double Bij = gsl_matrix_get(B, i, j);
3152
3153 gsl_test_abs(Bij, 0.0, tol, "tri_invert[%zu,%zu] N=%zu lower %s",
3154 i, j, n,
3155 (Diag == CblasNonUnit) ? "NonUnit" : "Unit");
3156 }
3157 }
3158 }
3159
3160 gsl_matrix_free(T);
3161 gsl_matrix_free(B);
3162 }
3163
3164 return s;
3165 }
3166
3167 int
test_tri_invert(gsl_rng * r)3168 test_tri_invert(gsl_rng * r)
3169 {
3170 int s = 0;
3171
3172 s += test_tri_invert2(CblasLower, CblasNonUnit, r, 1.0e-10);
3173 s += test_tri_invert2(CblasLower, CblasUnit, r, 1.0e-10);
3174
3175 s += test_tri_invert2(CblasUpper, CblasNonUnit, r, 1.0e-10);
3176 s += test_tri_invert2(CblasUpper, CblasUnit, r, 1.0e-10);
3177
3178 return s;
3179 }
3180
3181 void
my_error_handler(const char * reason,const char * file,int line,int err)3182 my_error_handler (const char *reason, const char *file, int line, int err)
3183 {
3184 if (1) printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err) ;
3185 }
3186
3187 int
main(void)3188 main(void)
3189 {
3190 gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
3191
3192 gsl_ieee_env_setup ();
3193 gsl_set_error_handler (&my_error_handler);
3194
3195 m11 = create_general_matrix(1,1);
3196 m51 = create_general_matrix(5,1);
3197
3198 m35 = create_general_matrix(3,5);
3199 m53 = create_general_matrix(5,3);
3200 m97 = create_general_matrix(9,7);
3201
3202 s35 = create_singular_matrix(3,5);
3203 s53 = create_singular_matrix(5,3);
3204
3205 hilb2 = create_hilbert_matrix(2);
3206 hilb3 = create_hilbert_matrix(3);
3207 hilb4 = create_hilbert_matrix(4);
3208 hilb12 = create_hilbert_matrix(12);
3209
3210 vander2 = create_vandermonde_matrix(2);
3211 vander3 = create_vandermonde_matrix(3);
3212 vander4 = create_vandermonde_matrix(4);
3213 vander12 = create_vandermonde_matrix(12);
3214
3215 moler10 = create_moler_matrix(10);
3216
3217 c7 = create_complex_matrix(7);
3218
3219 row3 = create_row_matrix(3,3);
3220 row5 = create_row_matrix(5,5);
3221 row12 = create_row_matrix(12,12);
3222
3223 A22 = create_2x2_matrix (0.0, 0.0, 0.0, 0.0);
3224 A33 = gsl_matrix_alloc(3,3);
3225 A44 = gsl_matrix_alloc(4,4);
3226 A55 = gsl_matrix_alloc(5,5);
3227
3228 inf5 = create_diagonal_matrix (inf5_data, 5);
3229 gsl_matrix_set(inf5, 3, 3, GSL_POSINF);
3230
3231 nan5 = create_diagonal_matrix (inf5_data, 5);
3232 gsl_matrix_set(nan5, 3, 3, GSL_NAN);
3233
3234 dblmin3 = create_general_matrix (3, 3);
3235 gsl_matrix_scale(dblmin3, GSL_DBL_MIN);
3236
3237 dblmin5 = create_general_matrix (5, 5);
3238 gsl_matrix_scale(dblmin5, GSL_DBL_MIN);
3239
3240 dblsubnorm5 = create_general_matrix (5, 5);
3241 gsl_matrix_scale(dblsubnorm5, GSL_DBL_MIN/1024);
3242 gsl_matrix_set(dblsubnorm5, 0, 0, 0.0);
3243
3244 bigsparse = create_sparse_matrix(100, 100);
3245
3246 /* Matmult now obsolete */
3247 #ifdef MATMULT
3248 gsl_test(test_matmult(), "Matrix Multiply");
3249 gsl_test(test_matmult_mod(), "Matrix Multiply with Modification");
3250 #endif
3251
3252 gsl_test(test_tri_invert(r), "Triangular Inverse");
3253
3254 gsl_test(test_bidiag_decomp(), "Bidiagonal Decomposition");
3255 gsl_test(test_LU_decomp(r), "LU Decomposition");
3256 gsl_test(test_LU_solve(r), "LU Solve");
3257 gsl_test(test_LU_invert(r), "LU Inverse");
3258 gsl_test(test_LUc_decomp(r), "Complex LU Decomposition");
3259 gsl_test(test_LUc_solve(r), "Complex LU Solve");
3260 gsl_test(test_LUc_invert(r), "Complex LU Inverse");
3261 gsl_test(test_QR_decomp(), "QR Decomposition");
3262 gsl_test(test_QR_solve(), "QR Solve");
3263 gsl_test(test_QRc_decomp(r), "Complex QR Decomposition");
3264 gsl_test(test_QRc_solve(r), "Complex QR Solve");
3265 gsl_test(test_QRc_lssolve(), "Complex QR LS Solve");
3266 gsl_test(test_LQ_solve(), "LQ Solve");
3267 gsl_test(test_PTLQ_solve(), "PTLQ Solve");
3268
3269 gsl_test(test_QR_decomp_r(r), "QR Decomposition (recursive)");
3270 gsl_test(test_QR_QTmat_r(r), "QR QTmat (recursive)");
3271 gsl_test(test_QR_solve_r(r), "QR Solve (recursive)");
3272 gsl_test(test_QR_lssolve_r(r), "QR LS Solve (recursive)");
3273
3274 gsl_test(test_QRc_decomp_r(r), "Complex QR Decomposition (recursive)");
3275 gsl_test(test_QRc_solve_r(r), "Complex QR Solve (recursive)");
3276 gsl_test(test_QRc_lssolve_r(), "Complex QR LS Solve (recursive)");
3277
3278 gsl_test(test_LU_band_decomp(r), "Banded LU Decomposition");
3279 gsl_test(test_LU_band_solve(r), "Banded LU Solve");
3280
3281 gsl_test(test_QR_band_decomp(r), "Banded QR Decomposition");
3282
3283 gsl_test(test_QR_UR_decomp(r), "QR_UR Decomposition");
3284 gsl_test(test_QR_UZ_decomp(r), "QR_UZ Decomposition");
3285 gsl_test(test_QR_UU_decomp(r), "QR_UU Decomposition");
3286 gsl_test(test_QR_UD_decomp(r), "QR_UD Decomposition");
3287
3288 gsl_test(test_QR_UU_lssolve(r), "QR_UU LS Solve");
3289 gsl_test(test_QR_UD_lssolve(r), "QR_UD LS Solve");
3290
3291 gsl_test(test_QL_decomp(r), "QL Decomposition");
3292
3293 gsl_test(test_LQ_decomp(), "LQ Decomposition");
3294 gsl_test(test_LQ_LQsolve(), "LQ LQ Solve");
3295 gsl_test(test_LQ_lssolve_T(), "LQ LS Solve_T");
3296 gsl_test(test_LQ_lssolve(), "LQ LS Solve");
3297 gsl_test(test_LQ_update(), "LQ Rank-1 Update");
3298 gsl_test(test_QRPT_decomp(), "PTLQ Decomposition");
3299 gsl_test(test_PTLQ_solve(), "PTLQ Solve");
3300
3301 gsl_test(test_QR_QRsolve(), "QR QR Solve");
3302 gsl_test(test_QR_lssolve(), "QR LS Solve");
3303 gsl_test(test_QR_update(), "QR Rank-1 Update");
3304 gsl_test(test_QRPT_decomp(), "QRPT Decomposition");
3305 gsl_test(test_QRPT_lssolve(), "QRPT LS Solve");
3306 gsl_test(test_QRPT_lssolve2(), "QRPT LS Solve 2");
3307 gsl_test(test_QRPT_solve(), "QRPT Solve");
3308 gsl_test(test_QRPT_QRsolve(), "QRPT QR Solve");
3309 gsl_test(test_QRPT_update(), "QRPT Rank-1 Update");
3310
3311 gsl_test(test_COD_decomp(r), "COD Decomposition");
3312 gsl_test(test_COD_lssolve(), "COD LS Solve");
3313 gsl_test(test_COD_lssolve2(r), "COD LS Solve 2");
3314
3315 gsl_test(test_SV_decomp(), "Singular Value Decomposition");
3316 gsl_test(test_SV_decomp_jacobi(), "Singular Value Decomposition (Jacobi)");
3317 gsl_test(test_SV_decomp_mod(), "Singular Value Decomposition (Mod)");
3318 gsl_test(test_SV_solve(), "SVD Solve");
3319
3320 gsl_test(test_cholesky_decomp_unit(), "Cholesky Decomposition [unit triangular]");
3321 gsl_test(test_cholesky_solve(), "Cholesky Solve");
3322 gsl_test(test_cholesky_decomp(r), "Cholesky Decomposition");
3323 gsl_test(test_cholesky_invert(r), "Cholesky Inverse");
3324
3325 gsl_test(test_pcholesky_decomp(r), "Pivoted Cholesky Decomposition");
3326 gsl_test(test_pcholesky_solve(r), "Pivoted Cholesky Solve");
3327 gsl_test(test_pcholesky_invert(r), "Pivoted Cholesky Inverse");
3328
3329 gsl_test(test_mcholesky_decomp(r), "Modified Cholesky Decomposition");
3330 gsl_test(test_mcholesky_solve(r), "Modified Cholesky Solve");
3331 gsl_test(test_mcholesky_invert(r), "Modified Cholesky Inverse");
3332
3333 gsl_test(test_choleskyc_decomp(r), "Complex Cholesky Decomposition");
3334 gsl_test(test_choleskyc_solve(), "Complex Cholesky Solve");
3335 gsl_test(test_choleskyc_invert(r), "Complex Cholesky Inverse");
3336
3337 gsl_test(test_cholesky_band_decomp(r), "Banded Cholesky Decomposition");
3338 gsl_test(test_cholesky_band_solve(r), "Banded Cholesky Solve");
3339 gsl_test(test_cholesky_band_invert(r), "Banded Cholesky Inverse");
3340
3341 gsl_test(test_ldlt_decomp(r), "LDLT Decomposition");
3342 gsl_test(test_ldlt_solve(r), "LDLT Solve");
3343
3344 gsl_test(test_ldlt_band_decomp(r), "Banded LDLT Decomposition");
3345 gsl_test(test_ldlt_band_solve(r), "Banded LDLT Solve");
3346
3347 gsl_test(test_symmtd_decomp(r), "Symmetric Tridiagonal Decomposition");
3348 gsl_test(test_hermtd_decomp(r), "Hermitian Tridiagonal Decomposition");
3349
3350 gsl_test(test_HH_solve(), "Householder solve");
3351 gsl_test(test_TDS_solve(), "Tridiagonal symmetric solve");
3352 gsl_test(test_TDS_cyc_solve(), "Tridiagonal symmetric cyclic solve");
3353 gsl_test(test_TDN_solve(), "Tridiagonal nonsymmetric solve");
3354 gsl_test(test_TDN_cyc_solve(), "Tridiagonal nonsymmetric cyclic solve");
3355
3356 gsl_matrix_free(m11);
3357 gsl_matrix_free(m35);
3358 gsl_matrix_free(m51);
3359 gsl_matrix_free(m53);
3360 gsl_matrix_free(m97);
3361 gsl_matrix_free(s35);
3362 gsl_matrix_free(s53);
3363
3364 gsl_matrix_free(hilb2);
3365 gsl_matrix_free(hilb3);
3366 gsl_matrix_free(hilb4);
3367 gsl_matrix_free(hilb12);
3368
3369 gsl_matrix_free(vander2);
3370 gsl_matrix_free(vander3);
3371 gsl_matrix_free(vander4);
3372 gsl_matrix_free(vander12);
3373
3374 gsl_matrix_free(moler10);
3375
3376 gsl_matrix_complex_free(c7);
3377 gsl_matrix_free(row3);
3378 gsl_matrix_free(row5);
3379 gsl_matrix_free(row12);
3380
3381 gsl_matrix_free(A22);
3382 gsl_matrix_free(A33);
3383 gsl_matrix_free(A44);
3384 gsl_matrix_free(A55);
3385
3386 gsl_matrix_free (inf5);
3387 gsl_matrix_free (nan5);
3388
3389 gsl_matrix_free (dblmin3);
3390 gsl_matrix_free (dblmin5);
3391 gsl_matrix_free (dblsubnorm5);
3392
3393 gsl_matrix_free (bigsparse);
3394
3395 gsl_rng_free(r);
3396
3397 exit (gsl_test_summary());
3398 }
3399