1 /* linalg/test_lq.c
2 *
3 * Copyright (C) 2018 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <stdlib.h>
22 #include <gsl/gsl_test.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_ieee_utils.h>
25 #include <gsl/gsl_permute_vector.h>
26 #include <gsl/gsl_blas.h>
27 #include <gsl/gsl_linalg.h>
28 #include <gsl/gsl_rng.h>
29 #include <gsl/gsl_permutation.h>
30
31 int test_LQ_solve_dim(const gsl_matrix * m, const double * actual, double eps);
32 int test_LQ_solve(void);
33 int test_LQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps);
34 int test_LQ_LQsolve(void);
35 int test_LQ_lssolve_T_dim(const gsl_matrix * m, const double * actual, double eps);
36 int test_LQ_lssolve_T(void);
37 int test_LQ_lssolve_dim(const gsl_matrix * m, const double * actual, double eps);
38 int test_LQ_lssolve(void);
39 int test_LQ_decomp_dim(const gsl_matrix * m, double eps);
40 int test_LQ_decomp(void);
41 int test_PTLQ_solve_dim(const gsl_matrix * m, const double * actual, double eps);
42 int test_PTLQ_solve(void);
43 int test_PTLQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps);
44 int test_PTLQ_LQsolve(void);
45 int test_PTLQ_decomp_dim(const gsl_matrix * m, double eps);
46 int test_PTLQ_decomp(void);
47 int test_LQ_update_dim(const gsl_matrix * m, double eps);
48 int test_LQ_update(void);
49
50 int
test_LQ_solve_dim(const gsl_matrix * m,const double * actual,double eps)51 test_LQ_solve_dim(const gsl_matrix * m, const double * actual, double eps)
52 {
53 int s = 0;
54 unsigned long i, dim = m->size1;
55
56 gsl_vector * rhs = gsl_vector_alloc(dim);
57 gsl_matrix * lq = gsl_matrix_alloc(dim,dim);
58 gsl_vector * d = gsl_vector_alloc(dim);
59 gsl_vector * x = gsl_vector_alloc(dim);
60
61 gsl_matrix_transpose_memcpy(lq,m);
62 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
63 s += gsl_linalg_LQ_decomp(lq, d);
64 s += gsl_linalg_LQ_solve_T(lq, d, rhs, x);
65 for(i=0; i<dim; i++) {
66 int foo = check(gsl_vector_get(x, i), actual[i], eps);
67 if(foo) {
68 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
69 }
70 s += foo;
71 }
72
73 gsl_vector_free(x);
74 gsl_vector_free(d);
75 gsl_matrix_free(lq);
76 gsl_vector_free(rhs);
77
78 return s;
79 }
80
81 int
test_LQ_solve(void)82 test_LQ_solve(void)
83 {
84 int f;
85 int s = 0;
86
87 f = test_LQ_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
88 gsl_test(f, " LQ_solve hilbert(2)");
89 s += f;
90
91 f = test_LQ_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
92 gsl_test(f, " LQ_solve hilbert(3)");
93 s += f;
94
95 f = test_LQ_solve_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
96 gsl_test(f, " LQ_solve hilbert(4)");
97 s += f;
98
99 f = test_LQ_solve_dim(hilb12, hilb12_solution, 0.5);
100 gsl_test(f, " LQ_solve hilbert(12)");
101 s += f;
102
103 f = test_LQ_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
104 gsl_test(f, " LQ_solve vander(2)");
105 s += f;
106
107 f = test_LQ_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
108 gsl_test(f, " LQ_solve vander(3)");
109 s += f;
110
111 f = test_LQ_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
112 gsl_test(f, " LQ_solve vander(4)");
113 s += f;
114
115 f = test_LQ_solve_dim(vander12, vander12_solution, 0.05);
116 gsl_test(f, " LQ_solve vander(12)");
117 s += f;
118
119 return s;
120 }
121
122 int
test_LQ_LQsolve_dim(const gsl_matrix * m,const double * actual,double eps)123 test_LQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps)
124 {
125 int s = 0;
126 unsigned long i, dim = m->size1;
127
128 gsl_vector * rhs = gsl_vector_alloc(dim);
129 gsl_matrix * lq = gsl_matrix_alloc(dim,dim);
130 gsl_matrix * q = gsl_matrix_alloc(dim,dim);
131 gsl_matrix * l = gsl_matrix_alloc(dim,dim);
132 gsl_vector * d = gsl_vector_alloc(dim);
133 gsl_vector * x = gsl_vector_alloc(dim);
134
135 gsl_matrix_transpose_memcpy(lq,m);
136 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
137 s += gsl_linalg_LQ_decomp(lq, d);
138 s += gsl_linalg_LQ_unpack(lq, d, q, l);
139 s += gsl_linalg_LQ_LQsolve(q, l, rhs, x);
140 for(i=0; i<dim; i++) {
141 int foo = check(gsl_vector_get(x, i), actual[i], eps);
142 if(foo) {
143 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
144 }
145 s += foo;
146 }
147
148 gsl_vector_free(x);
149 gsl_vector_free(d);
150 gsl_matrix_free(lq);
151 gsl_matrix_free(q);
152 gsl_matrix_free(l);
153 gsl_vector_free(rhs);
154
155 return s;
156 }
157
158 int
test_LQ_LQsolve(void)159 test_LQ_LQsolve(void)
160 {
161 int f;
162 int s = 0;
163
164 f = test_LQ_LQsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
165 gsl_test(f, " LQ_LQsolve hilbert(2)");
166 s += f;
167
168 f = test_LQ_LQsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
169 gsl_test(f, " LQ_LQsolve hilbert(3)");
170 s += f;
171
172 f = test_LQ_LQsolve_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
173 gsl_test(f, " LQ_LQsolve hilbert(4)");
174 s += f;
175
176 f = test_LQ_LQsolve_dim(hilb12, hilb12_solution, 0.5);
177 gsl_test(f, " LQ_LQsolve hilbert(12)");
178 s += f;
179
180 f = test_LQ_LQsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
181 gsl_test(f, " LQ_LQsolve vander(2)");
182 s += f;
183
184 f = test_LQ_LQsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
185 gsl_test(f, " LQ_LQsolve vander(3)");
186 s += f;
187
188 f = test_LQ_LQsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
189 gsl_test(f, " LQ_LQsolve vander(4)");
190 s += f;
191
192 f = test_LQ_LQsolve_dim(vander12, vander12_solution, 0.05);
193 gsl_test(f, " LQ_LQsolve vander(12)");
194 s += f;
195
196 return s;
197 }
198
199 int
test_LQ_lssolve_T_dim(const gsl_matrix * m,const double * actual,double eps)200 test_LQ_lssolve_T_dim(const gsl_matrix * m, const double * actual, double eps)
201 {
202 int s = 0;
203 unsigned long i, M = m->size1, N = m->size2;
204
205 gsl_vector * rhs = gsl_vector_alloc(M);
206 gsl_matrix * lq = gsl_matrix_alloc(N,M);
207 gsl_vector * d = gsl_vector_alloc(N);
208 gsl_vector * x = gsl_vector_alloc(N);
209 gsl_vector * r = gsl_vector_alloc(M);
210 gsl_vector * res = gsl_vector_alloc(M);
211
212 gsl_matrix_transpose_memcpy(lq,m);
213 for(i=0; i<M; i++) gsl_vector_set(rhs, i, i+1.0);
214 s += gsl_linalg_LQ_decomp(lq, d);
215 s += gsl_linalg_LQ_lssolve_T(lq, d, rhs, x, res);
216
217 for(i=0; i<N; i++) {
218 int foo = check(gsl_vector_get(x, i), actual[i], eps);
219 if(foo) {
220 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
221 }
222 s += foo;
223 }
224
225
226 /* compute residual r = b - m x */
227 if (M == N) {
228 gsl_vector_set_zero(r);
229 } else {
230 gsl_vector_memcpy(r, rhs);
231 gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
232 };
233
234 for(i=0; i<N; i++) {
235 int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
236 if(foo) {
237 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
238 }
239 s += foo;
240 }
241
242 gsl_vector_free(r);
243 gsl_vector_free(res);
244 gsl_vector_free(x);
245 gsl_vector_free(d);
246 gsl_matrix_free(lq);
247 gsl_vector_free(rhs);
248
249 return s;
250 }
251
252 int
test_LQ_lssolve_dim(const gsl_matrix * m,const double * actual,double eps)253 test_LQ_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
254 {
255 int s = 0;
256 unsigned long i, M = m->size1, N = m->size2;
257
258 gsl_vector * rhs = gsl_vector_alloc(M);
259 gsl_matrix * lq = gsl_matrix_alloc(M, N);
260 gsl_vector * tau = gsl_vector_alloc(N);
261 gsl_vector * x = gsl_vector_alloc(N);
262 gsl_vector * r = gsl_vector_alloc(M);
263 gsl_vector * res = gsl_vector_alloc(M);
264
265 for (i = 0; i < M; i++)
266 gsl_vector_set(rhs, i, i + 1.0);
267
268 gsl_matrix_memcpy(lq, m);
269 gsl_linalg_LQ_decomp(lq, tau);
270 gsl_linalg_LQ_lssolve(lq, tau, rhs, x, res);
271
272 for (i = 0; i < N; i++)
273 {
274 int foo = check(gsl_vector_get(x, i), actual[i], eps);
275 if(foo) {
276 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
277 }
278 s += foo;
279 }
280
281 /* compute residual r = b - m x */
282 if (M == N) {
283 gsl_vector_set_zero(r);
284 } else {
285 gsl_vector_memcpy(r, rhs);
286 gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
287 };
288
289 for (i = 0; i < N; i++)
290 {
291 int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
292 if(foo) {
293 printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
294 }
295 s += foo;
296 }
297
298 gsl_vector_free(r);
299 gsl_vector_free(res);
300 gsl_vector_free(x);
301 gsl_vector_free(tau);
302 gsl_matrix_free(lq);
303 gsl_vector_free(rhs);
304
305 return s;
306 }
307
308 int
test_LQ_lssolve_T(void)309 test_LQ_lssolve_T(void)
310 {
311 int f;
312 int s = 0;
313
314 f = test_LQ_lssolve_T_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
315 gsl_test(f, " LQ_lssolve_T m(5,3)");
316 s += f;
317
318 f = test_LQ_lssolve_T_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
319 gsl_test(f, " LQ_lssolve_T hilbert(2)");
320 s += f;
321
322 f = test_LQ_lssolve_T_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
323 gsl_test(f, " LQ_lssolve_T hilbert(3)");
324 s += f;
325
326 f = test_LQ_lssolve_T_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
327 gsl_test(f, " LQ_lssolve_T hilbert(4)");
328 s += f;
329
330 f = test_LQ_lssolve_T_dim(hilb12, hilb12_solution, 0.5);
331 gsl_test(f, " LQ_lssolve_T hilbert(12)");
332 s += f;
333
334 f = test_LQ_lssolve_T_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
335 gsl_test(f, " LQ_lssolve_T vander(2)");
336 s += f;
337
338 f = test_LQ_lssolve_T_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
339 gsl_test(f, " LQ_lssolve_T vander(3)");
340 s += f;
341
342 f = test_LQ_lssolve_T_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
343 gsl_test(f, " LQ_lssolve_T vander(4)");
344 s += f;
345
346 f = test_LQ_lssolve_T_dim(vander12, vander12_solution, 0.05);
347 gsl_test(f, " LQ_lssolve_T vander(12)");
348 s += f;
349
350 return s;
351 }
352
353 int
test_LQ_lssolve(void)354 test_LQ_lssolve(void)
355 {
356 int f;
357 int s = 0;
358
359 f = test_LQ_lssolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
360 gsl_test(f, " LQ_lssolve hilbert(2)");
361 s += f;
362
363 f = test_LQ_lssolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
364 gsl_test(f, " LQ_lssolve hilbert(3)");
365 s += f;
366
367 f = test_LQ_lssolve_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
368 gsl_test(f, " LQ_lssolve hilbert(4)");
369 s += f;
370
371 f = test_LQ_lssolve_dim(hilb12, hilb12_solution, 0.5);
372 gsl_test(f, " LQ_lssolve hilbert(12)");
373 s += f;
374
375 f = test_LQ_lssolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
376 gsl_test(f, " LQ_lssolve vander(2)");
377 s += f;
378
379 f = test_LQ_lssolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
380 gsl_test(f, " LQ_lssolve vander(3)");
381 s += f;
382
383 f = test_LQ_lssolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
384 gsl_test(f, " LQ_lssolve vander(4)");
385 s += f;
386
387 f = test_LQ_lssolve_dim(vander12, vander12_solution, 0.05);
388 gsl_test(f, " LQ_lssolve vander(12)");
389 s += f;
390
391 return s;
392 }
393
394 int
test_LQ_decomp_dim(const gsl_matrix * m,double eps)395 test_LQ_decomp_dim(const gsl_matrix * m, double eps)
396 {
397 int s = 0;
398 unsigned long i,j, M = m->size1, N = m->size2;
399
400 gsl_matrix * lq = gsl_matrix_alloc(M,N);
401 gsl_matrix * a = gsl_matrix_alloc(M,N);
402 gsl_matrix * q = gsl_matrix_alloc(N,N);
403 gsl_matrix * l = gsl_matrix_alloc(M,N);
404 gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
405
406 gsl_matrix_memcpy(lq,m);
407
408 s += gsl_linalg_LQ_decomp(lq, d);
409 s += gsl_linalg_LQ_unpack(lq, d, q, l);
410
411 /* compute a = q r */
412 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, q, 0.0, a);
413
414 for(i=0; i<M; i++) {
415 for(j=0; j<N; j++) {
416 double aij = gsl_matrix_get(a, i, j);
417 double mij = gsl_matrix_get(m, i, j);
418 int foo = check(aij, mij, eps);
419 if(foo) {
420 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
421 }
422 s += foo;
423 }
424 }
425
426 gsl_vector_free(d);
427 gsl_matrix_free(lq);
428 gsl_matrix_free(a);
429 gsl_matrix_free(q);
430 gsl_matrix_free(l);
431
432 return s;
433 }
434
435 int
test_LQ_decomp(void)436 test_LQ_decomp(void)
437 {
438 int f;
439 int s = 0;
440
441 f = test_LQ_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
442 gsl_test(f, " LQ_decomp m(3,5)");
443 s += f;
444
445 f = test_LQ_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
446 gsl_test(f, " LQ_decomp m(5,3)");
447 s += f;
448
449 f = test_LQ_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
450 gsl_test(f, " LQ_decomp hilbert(2)");
451 s += f;
452
453 f = test_LQ_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
454 gsl_test(f, " LQ_decomp hilbert(3)");
455 s += f;
456
457 f = test_LQ_decomp_dim(hilb4, 4 * 1024.0 * GSL_DBL_EPSILON);
458 gsl_test(f, " LQ_decomp hilbert(4)");
459 s += f;
460
461 f = test_LQ_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
462 gsl_test(f, " LQ_decomp hilbert(12)");
463 s += f;
464
465 f = test_LQ_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
466 gsl_test(f, " LQ_decomp vander(2)");
467 s += f;
468
469 f = test_LQ_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
470 gsl_test(f, " LQ_decomp vander(3)");
471 s += f;
472
473 f = test_LQ_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
474 gsl_test(f, " LQ_decomp vander(4)");
475 s += f;
476
477 f = test_LQ_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
478 gsl_test(f, " LQ_decomp vander(12)");
479 s += f;
480
481 return s;
482 }
483
484 int
test_PTLQ_solve_dim(const gsl_matrix * m,const double * actual,double eps)485 test_PTLQ_solve_dim(const gsl_matrix * m, const double * actual, double eps)
486 {
487 int s = 0;
488 int signum;
489 unsigned long i, dim = m->size1;
490
491 gsl_permutation * perm = gsl_permutation_alloc(dim);
492 gsl_vector * rhs = gsl_vector_alloc(dim);
493 gsl_matrix * lq = gsl_matrix_alloc(dim,dim);
494 gsl_vector * d = gsl_vector_alloc(dim);
495 gsl_vector * x = gsl_vector_alloc(dim);
496 gsl_vector * norm = gsl_vector_alloc(dim);
497
498 gsl_matrix_transpose_memcpy(lq,m);
499 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
500 s += gsl_linalg_PTLQ_decomp(lq, d, perm, &signum, norm);
501 s += gsl_linalg_PTLQ_solve_T(lq, d, perm, rhs, x);
502 for(i=0; i<dim; i++) {
503 int foo = check(gsl_vector_get(x, i), actual[i], eps);
504 if(foo) {
505 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
506 }
507 s += foo;
508 }
509
510 gsl_vector_free(norm);
511 gsl_vector_free(x);
512 gsl_vector_free(d);
513 gsl_matrix_free(lq);
514 gsl_vector_free(rhs);
515 gsl_permutation_free(perm);
516
517 return s;
518 }
519
test_PTLQ_solve(void)520 int test_PTLQ_solve(void)
521 {
522 int f;
523 int s = 0;
524
525 f = test_PTLQ_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
526 gsl_test(f, " PTLQ_solve hilbert(2)");
527 s += f;
528
529 f = test_PTLQ_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
530 gsl_test(f, " PTLQ_solve hilbert(3)");
531 s += f;
532
533 f = test_PTLQ_solve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
534 gsl_test(f, " PTLQ_solve hilbert(4)");
535 s += f;
536
537 f = test_PTLQ_solve_dim(hilb12, hilb12_solution, 0.5);
538 gsl_test(f, " PTLQ_solve hilbert(12)");
539 s += f;
540
541 f = test_PTLQ_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
542 gsl_test(f, " PTLQ_solve vander(2)");
543 s += f;
544
545 f = test_PTLQ_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
546 gsl_test(f, " PTLQ_solve vander(3)");
547 s += f;
548
549 f = test_PTLQ_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
550 gsl_test(f, " PTLQ_solve vander(4)");
551 s += f;
552
553 f = test_PTLQ_solve_dim(vander12, vander12_solution, 0.05);
554 gsl_test(f, " PTLQ_solve vander(12)");
555 s += f;
556
557 return s;
558 }
559
560 int
test_PTLQ_LQsolve_dim(const gsl_matrix * m,const double * actual,double eps)561 test_PTLQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps)
562 {
563 int s = 0;
564 int signum;
565 unsigned long i, dim = m->size1;
566
567 gsl_permutation * perm = gsl_permutation_alloc(dim);
568 gsl_vector * rhs = gsl_vector_alloc(dim);
569 gsl_matrix * lq = gsl_matrix_alloc(dim,dim);
570 gsl_matrix * q = gsl_matrix_alloc(dim,dim);
571 gsl_matrix * l = gsl_matrix_alloc(dim,dim);
572 gsl_vector * d = gsl_vector_alloc(dim);
573 gsl_vector * x = gsl_vector_alloc(dim);
574 gsl_vector * norm = gsl_vector_alloc(dim);
575
576 gsl_matrix_transpose_memcpy(lq,m);
577 for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
578 s += gsl_linalg_PTLQ_decomp2(lq, q, l, d, perm, &signum, norm);
579 s += gsl_linalg_PTLQ_LQsolve_T(q, l, perm, rhs, x);
580 for(i=0; i<dim; i++) {
581 int foo = check(gsl_vector_get(x, i), actual[i], eps);
582 if(foo) {
583 printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
584 }
585 s += foo;
586 }
587
588 gsl_vector_free(norm);
589 gsl_vector_free(x);
590 gsl_vector_free(d);
591 gsl_matrix_free(lq);
592 gsl_vector_free(rhs);
593 gsl_permutation_free(perm);
594
595 return s;
596 }
597
test_PTLQ_LQsolve(void)598 int test_PTLQ_LQsolve(void)
599 {
600 int f;
601 int s = 0;
602
603 f = test_PTLQ_LQsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
604 gsl_test(f, " PTLQ_LQsolve hilbert(2)");
605 s += f;
606
607 f = test_PTLQ_LQsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
608 gsl_test(f, " PTLQ_LQsolve hilbert(3)");
609 s += f;
610
611 f = test_PTLQ_LQsolve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
612 gsl_test(f, " PTLQ_LQsolve hilbert(4)");
613 s += f;
614
615 f = test_PTLQ_LQsolve_dim(hilb12, hilb12_solution, 0.5);
616 gsl_test(f, " PTLQ_LQsolve hilbert(12)");
617 s += f;
618
619 f = test_PTLQ_LQsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
620 gsl_test(f, " PTLQ_LQsolve vander(2)");
621 s += f;
622
623 f = test_PTLQ_LQsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
624 gsl_test(f, " PTLQ_LQsolve vander(3)");
625 s += f;
626
627 f = test_PTLQ_LQsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
628 gsl_test(f, " PTLQ_LQsolve vander(4)");
629 s += f;
630
631 f = test_PTLQ_LQsolve_dim(vander12, vander12_solution, 0.05);
632 gsl_test(f, " PTLQ_LQsolve vander(12)");
633 s += f;
634
635 return s;
636 }
637
638 int
test_PTLQ_decomp_dim(const gsl_matrix * m,double eps)639 test_PTLQ_decomp_dim(const gsl_matrix * m, double eps)
640 {
641 int s = 0, signum;
642 unsigned long i,j, M = m->size1, N = m->size2;
643
644 gsl_matrix * lq = gsl_matrix_alloc(N,M);
645 gsl_matrix * a = gsl_matrix_alloc(N,M);
646 gsl_matrix * q = gsl_matrix_alloc(M,M);
647 gsl_matrix * l = gsl_matrix_alloc(N,M);
648 gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
649 gsl_vector * norm = gsl_vector_alloc(N);
650
651 gsl_permutation * perm = gsl_permutation_alloc(N);
652
653 gsl_matrix_transpose_memcpy(lq,m);
654
655 s += gsl_linalg_PTLQ_decomp(lq, d, perm, &signum, norm);
656 s += gsl_linalg_LQ_unpack(lq, d, q, l);
657
658 /* compute a = l q */
659 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, q, 0.0, a);
660
661
662 /* Compute P LQ by permuting the rows of LQ */
663
664 for (i = 0; i < M; i++) {
665 gsl_vector_view col = gsl_matrix_column (a, i);
666 gsl_permute_vector_inverse (perm, &col.vector);
667 }
668
669 for(i=0; i<M; i++) {
670 for(j=0; j<N; j++) {
671 double aij = gsl_matrix_get(a, j, i);
672 double mij = gsl_matrix_get(m, i, j);
673 int foo = check(aij, mij, eps);
674 if(foo) {
675 printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, aij, mij);
676 }
677 s += foo;
678 }
679 }
680
681 gsl_permutation_free (perm);
682 gsl_vector_free(norm);
683 gsl_vector_free(d);
684 gsl_matrix_free(lq);
685 gsl_matrix_free(a);
686 gsl_matrix_free(q);
687 gsl_matrix_free(l);
688
689 return s;
690 }
691
test_PTLQ_decomp(void)692 int test_PTLQ_decomp(void)
693 {
694 int f;
695 int s = 0;
696
697 f = test_PTLQ_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
698 gsl_test(f, " PTLQ_decomp m(3,5)");
699 s += f;
700
701 f = test_PTLQ_decomp_dim(m53, 2 * 8.0 * GSL_DBL_EPSILON);
702 gsl_test(f, " PTLQ_decomp m(5,3)");
703 s += f;
704
705 f = test_PTLQ_decomp_dim(s35, 2 * 8.0 * GSL_DBL_EPSILON);
706 gsl_test(f, " PTLQ_decomp s(3,5)");
707 s += f;
708
709 f = test_PTLQ_decomp_dim(s53, 2 * 8.0 * GSL_DBL_EPSILON);
710 gsl_test(f, " PTLQ_decomp s(5,3)");
711 s += f;
712
713 f = test_PTLQ_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
714 gsl_test(f, " PTLQ_decomp hilbert(2)");
715 s += f;
716
717 f = test_PTLQ_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
718 gsl_test(f, " PTLQ_decomp hilbert(3)");
719 s += f;
720
721 f = test_PTLQ_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
722 gsl_test(f, " PTLQ_decomp hilbert(4)");
723 s += f;
724
725 f = test_PTLQ_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
726 gsl_test(f, " PTLQ_decomp hilbert(12)");
727 s += f;
728
729 f = test_PTLQ_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
730 gsl_test(f, " PTLQ_decomp vander(2)");
731 s += f;
732
733 f = test_PTLQ_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
734 gsl_test(f, " PTLQ_decomp vander(3)");
735 s += f;
736
737 f = test_PTLQ_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
738 gsl_test(f, " PTLQ_decomp vander(4)");
739 s += f;
740
741 f = test_PTLQ_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
742 gsl_test(f, " PTLQ_decomp vander(12)");
743 s += f;
744
745 return s;
746 }
747
748 int
test_LQ_update_dim(const gsl_matrix * m,double eps)749 test_LQ_update_dim(const gsl_matrix * m, double eps)
750 {
751 int s = 0;
752 unsigned long i,j, M = m->size1, N = m->size2;
753
754 gsl_matrix * lq1 = gsl_matrix_alloc(N,M);
755 gsl_matrix * lq2 = gsl_matrix_alloc(N,M);
756 gsl_matrix * q1 = gsl_matrix_alloc(M,M);
757 gsl_matrix * l1 = gsl_matrix_alloc(N,M);
758 gsl_matrix * q2 = gsl_matrix_alloc(M,M);
759 gsl_matrix * l2 = gsl_matrix_alloc(N,M);
760 gsl_vector * d2 = gsl_vector_alloc(GSL_MIN(M,N));
761 gsl_vector * u = gsl_vector_alloc(M);
762 gsl_vector * v = gsl_vector_alloc(N);
763 gsl_vector * w = gsl_vector_alloc(M);
764
765 gsl_matrix_transpose_memcpy(lq1,m);
766 gsl_matrix_transpose_memcpy(lq2,m);
767 for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
768 for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
769
770 /* lq1 is updated */
771
772 gsl_blas_dger(1.0, v, u, lq1);
773
774 /* lq2 is first decomposed, updated later */
775
776 s += gsl_linalg_LQ_decomp(lq2, d2);
777 s += gsl_linalg_LQ_unpack(lq2, d2, q2, l2);
778
779 /* compute w = Q^T u */
780
781 gsl_blas_dgemv(CblasNoTrans, 1.0, q2, u, 0.0, w);
782
783 /* now lq2 is updated */
784
785 s += gsl_linalg_LQ_update(q2, l2, v, w);
786
787 /* multiply q2*l2 */
788
789 gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,l2,q2,0.0,lq2);
790
791 /* check lq1==lq2 */
792
793 for(i=0; i<N; i++) {
794 for(j=0; j<M; j++) {
795 double s1 = gsl_matrix_get(lq1, i, j);
796 double s2 = gsl_matrix_get(lq2, i, j);
797
798 int foo = check(s1, s2, eps);
799 #if 0
800 if(foo) {
801 printf("LQ:(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j, s1, s2);
802 }
803 #endif
804 s += foo;
805 }
806 }
807
808 gsl_vector_free(d2);
809 gsl_vector_free(u);
810 gsl_vector_free(v);
811 gsl_vector_free(w);
812 gsl_matrix_free(lq1);
813 gsl_matrix_free(lq2);
814 gsl_matrix_free(q1);
815 gsl_matrix_free(l1);
816 gsl_matrix_free(q2);
817 gsl_matrix_free(l2);
818
819 return s;
820 }
821
822 int
test_LQ_update(void)823 test_LQ_update(void)
824 {
825 int f;
826 int s = 0;
827
828 f = test_LQ_update_dim(m35, 2 * 512.0 * GSL_DBL_EPSILON);
829 gsl_test(f, " LQ_update m(3,5)");
830 s += f;
831
832 f = test_LQ_update_dim(m53, 2 * 2048.0 * GSL_DBL_EPSILON);
833 gsl_test(f, " LQ_update m(5,3)");
834 s += f;
835
836 f = test_LQ_update_dim(hilb2, 2 * 512.0 * GSL_DBL_EPSILON);
837 gsl_test(f, " LQ_update hilbert(2)");
838 s += f;
839
840 f = test_LQ_update_dim(hilb3, 2 * 512.0 * GSL_DBL_EPSILON);
841 gsl_test(f, " LQ_update hilbert(3)");
842 s += f;
843
844 f = test_LQ_update_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
845 gsl_test(f, " LQ_update hilbert(4)");
846 s += f;
847
848 f = test_LQ_update_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
849 gsl_test(f, " LQ_update hilbert(12)");
850 s += f;
851
852 f = test_LQ_update_dim(vander2, 8.0 * GSL_DBL_EPSILON);
853 gsl_test(f, " LQ_update vander(2)");
854 s += f;
855
856 f = test_LQ_update_dim(vander3, 64.0 * GSL_DBL_EPSILON);
857 gsl_test(f, " LQ_update vander(3)");
858 s += f;
859
860 f = test_LQ_update_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
861 gsl_test(f, " LQ_update vander(4)");
862 s += f;
863
864 f = test_LQ_update_dim(vander12, 0.0005); /* FIXME: bad accuracy */
865 gsl_test(f, " LQ_update vander(12)");
866 s += f;
867
868 return s;
869 }
870