1 package Math::GSL::Linalg::Test;
2 use base q{Test::Class};
3 use Test::More tests => 78;
4 use Math::GSL qw/:all/;
5 use Math::GSL::BLAS qw/:all/;
6 use Math::GSL::Test qw/:all/;
7 use Math::GSL::CBLAS qw/:all/;
8 use Math::GSL::Errno qw/:all/;
9 use Math::GSL::Linalg qw/:all/;
10 use Math::GSL::Matrix qw/:all/;
11 use Math::GSL::MatrixComplex qw/:all/;
12 use Math::GSL::Vector qw/:all/;
13 use Math::GSL::Machine qw/:all/;
14 use Math::GSL::Complex qw/:all/;
15 use Math::GSL::Permutation qw/:all/;
20 sub make_fixture : Test(setup) {
22 $self->{matrix} = gsl_matrix_alloc(4, 4);
25 sub teardown : Test(teardown) {
26 unlink 'linalg' if -f 'linalg';
29 sub GSL_LINALG_LU_DECOMP : Tests {
30 my $base = Math::GSL::Matrix->new(4,4);
31 $base->set_row(0, [0,1,2,3])
32 ->set_row(1, [5,6,7,8])
33 ->set_row(2, [9,10,11,12])
34 ->set_row(3, [13,14,15,16]);
35 my $permutation = gsl_permutation_alloc(4);
36 gsl_permutation_init($permutation);
37 my $first = Math::GSL::Matrix->new(4,4);
38 gsl_matrix_memcpy($first->raw, $base->raw);
40 my ($result, $signum) = gsl_linalg_LU_decomp($base->raw, $permutation);
41 is_deeply( [ $result, $signum ], [ 0, 1] );
43 my $U = Math::GSL::Matrix->new(4,4);
44 my $R = Math::GSL::Matrix->new(4,4);
45 my $L = Math::GSL::Matrix->new(4,4);
46 gsl_matrix_set_identity($L->raw);
48 for ($line=3; $line>-1; $line--) {
49 map { gsl_matrix_set($U->raw, $line, $_, gsl_matrix_get($base->raw, $line, $_)) } ($line..3) };
50 for ($line=3; $line>1; $line--) {
51 map { gsl_matrix_set($L->raw, $line, $_, gsl_matrix_get($base->raw, $line, $_)) } (0..$line-1) };
52 gsl_blas_dgemm($CblasNoTrans, $CblasNoTrans, 1, $L->raw, $U->raw, 1, $R->raw);
53 my @permutations = map { gsl_permutation_get($permutation, $_) } (0..3);
55 my @got = $first->row($permutations[$_])->as_list;
56 my @results = $R->row($_)->as_list;
57 ok_similar([@got], [@results], "resulting row $_"); } (0..3);
60 sub GSL_LINALG_LU_SOLVE : Tests {
62 gsl_matrix_set($self->{matrix}, 0, 0, 1);
63 gsl_matrix_set($self->{matrix}, 0, 1, 1);
64 gsl_matrix_set($self->{matrix}, 0, 2, 2);
65 gsl_matrix_set($self->{matrix}, 0, 3, 1);
67 gsl_matrix_set($self->{matrix}, 1, 0, 2);
68 gsl_matrix_set($self->{matrix}, 1, 1, 3);
69 gsl_matrix_set($self->{matrix}, 1, 2, -1);
70 gsl_matrix_set($self->{matrix}, 1, 3, 2);
72 gsl_matrix_set($self->{matrix}, 2, 0, 5);
73 gsl_matrix_set($self->{matrix}, 2, 1, -1);
74 gsl_matrix_set($self->{matrix}, 2, 2, 1);
75 gsl_matrix_set($self->{matrix}, 2, 3, -1);
77 gsl_matrix_set($self->{matrix}, 3, 0, 1);
78 gsl_matrix_set($self->{matrix}, 3, 1, 0);
79 gsl_matrix_set($self->{matrix}, 3, 2, 7);
80 gsl_matrix_set($self->{matrix}, 3, 3, 1);
82 my $b = gsl_vector_alloc(4);
83 gsl_vector_set($b, 0, 4);
84 gsl_vector_set($b, 1, 1);
85 gsl_vector_set($b, 2, 2);
86 gsl_vector_set($b, 3, 11);
88 my $x = gsl_vector_alloc(4);
90 my $permutation = gsl_permutation_alloc(4);
91 gsl_permutation_init($permutation);
92 gsl_linalg_LU_decomp($self->{matrix}, $permutation);
93 gsl_linalg_LU_solve($self->{matrix}, $permutation, $b, $x);
94 my $value = gsl_vector_get($x, 0);
96 [ map { gsl_vector_get($x, $_) } (1..3) ],
97 [ 3-10*$value, -2*$value+2, 13*$value-3 ]
101 sub GSL_LINALG_LU_SVX : Tests {
103 gsl_matrix_set($self->{matrix}, 0, 0, 1);
104 gsl_matrix_set($self->{matrix}, 0, 1, 1);
105 gsl_matrix_set($self->{matrix}, 0, 2, 2);
106 gsl_matrix_set($self->{matrix}, 0, 3, 1);
108 gsl_matrix_set($self->{matrix}, 1, 0, 2);
109 gsl_matrix_set($self->{matrix}, 1, 1, 3);
110 gsl_matrix_set($self->{matrix}, 1, 2, -1);
111 gsl_matrix_set($self->{matrix}, 1, 3, 2);
113 gsl_matrix_set($self->{matrix}, 2, 0, 5);
114 gsl_matrix_set($self->{matrix}, 2, 1, -1);
115 gsl_matrix_set($self->{matrix}, 2, 2, 1);
116 gsl_matrix_set($self->{matrix}, 2, 3, -1);
118 gsl_matrix_set($self->{matrix}, 3, 0, 1);
119 gsl_matrix_set($self->{matrix}, 3, 1, 0);
120 gsl_matrix_set($self->{matrix}, 3, 2, 7);
121 gsl_matrix_set($self->{matrix}, 3, 3, 1);
123 my $x = gsl_vector_alloc(4);
124 gsl_vector_set($x, 0, 4);
125 gsl_vector_set($x, 1, 1);
126 gsl_vector_set($x, 2, 2);
127 gsl_vector_set($x, 3, 11);
129 my $permutation = gsl_permutation_alloc(4);
130 gsl_permutation_init($permutation);
131 gsl_linalg_LU_decomp($self->{matrix}, $permutation);
132 gsl_linalg_LU_svx($self->{matrix}, $permutation, $x);
133 my $value = gsl_vector_get($x, 0);
135 [ map { gsl_vector_get($x, $_) } (1..3) ],
136 [ 3-10*$value, -2*$value+2, 13*$value-3 ]
140 sub GSL_LINALG_LU_INVERT : Tests {
142 map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
144 gsl_matrix_set($self->{matrix}, 1, 0, 2);
145 gsl_matrix_set($self->{matrix}, 1, 1, 3);
146 gsl_matrix_set($self->{matrix}, 1, 2, 4);
147 gsl_matrix_set($self->{matrix}, 1, 3, 1);
149 gsl_matrix_set($self->{matrix}, 2, 0, 3);
150 gsl_matrix_set($self->{matrix}, 2, 1, 4);
151 gsl_matrix_set($self->{matrix}, 2, 2, 1);
152 gsl_matrix_set($self->{matrix}, 2, 3, 2);
154 gsl_matrix_set($self->{matrix}, 3, 0, 4);
155 gsl_matrix_set($self->{matrix}, 3, 1, 1);
156 gsl_matrix_set($self->{matrix}, 3, 2, 2);
157 gsl_matrix_set($self->{matrix}, 3, 3, 3);
159 my $inverse = gsl_matrix_alloc(4,4);
160 my $permutation = gsl_permutation_alloc(4);
161 gsl_permutation_init($permutation);
162 gsl_linalg_LU_decomp($self->{matrix}, $permutation);
163 gsl_linalg_LU_invert($self->{matrix}, $permutation, $inverse);
165 ok_similar(gsl_matrix_get($inverse, 0, 0), -9/40);
166 ok_similar(gsl_matrix_get($inverse, 0, 1), 1/40);
167 ok_similar(gsl_matrix_get($inverse, 0, 2), 1/40);
168 ok_similar(gsl_matrix_get($inverse, 0, 3), 11/40);
170 ok_similar(gsl_matrix_get($inverse, 1, 0), 1/40);
171 ok_similar(gsl_matrix_get($inverse, 1, 1), 1/40);
172 ok_similar(gsl_matrix_get($inverse, 1, 2), 11/40);
173 ok_similar(gsl_matrix_get($inverse, 1, 3), -9/40);
175 ok_similar(gsl_matrix_get($inverse, 2, 0), 1/40);
176 ok_similar(gsl_matrix_get($inverse, 2, 1), 11/40);
177 ok_similar(gsl_matrix_get($inverse, 2, 2), -9/40);
178 ok_similar(gsl_matrix_get($inverse, 2, 3), 1/40);
180 ok_similar(gsl_matrix_get($inverse, 3, 0), 11/40);
181 ok_similar(gsl_matrix_get($inverse, 3, 1), -9/40);
182 ok_similar(gsl_matrix_get($inverse, 3, 2), 1/40);
183 ok_similar(gsl_matrix_get($inverse, 3, 3), 1/40);
186 sub GSL_LINALG_COMPLEX_LU_DET : Tests(2) {
187 local $TODO = q{ $gsl_det is wrong or Math::Complex conversion is };
188 my $m = Math::GSL::MatrixComplex->new(2,2)
189 ->set_row(0, [ 1 , 2 ])
190 ->set_row(1, [ 3 , 4*i ]);
191 my $permutation = gsl_permutation_alloc(2);
192 gsl_permutation_init($permutation);
194 my ($result, $signum) = gsl_linalg_complex_LU_decomp($copy->raw, $permutation);
195 my $gsl_det = gsl_linalg_complex_LU_det($m->raw, $signum);
196 my $det = Math::Complex->make( gsl_real($gsl_det), gsl_imag($gsl_det) );
197 #warn Dumper [ $det ];
198 ok_similar( [ Re $det ], [ -6 ], 'real');
199 ok_similar( [ Im $det ], [ 4 ], 'imag');
202 sub GSL_LINALG_LU_DET : Tests {
204 map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
206 gsl_matrix_set($self->{matrix}, 1, 0, 2);
207 gsl_matrix_set($self->{matrix}, 1, 1, 3);
208 gsl_matrix_set($self->{matrix}, 1, 2, 4);
209 gsl_matrix_set($self->{matrix}, 1, 3, 1);
211 gsl_matrix_set($self->{matrix}, 2, 0, 3);
212 gsl_matrix_set($self->{matrix}, 2, 1, 4);
213 gsl_matrix_set($self->{matrix}, 2, 2, 1);
214 gsl_matrix_set($self->{matrix}, 2, 3, 2);
216 gsl_matrix_set($self->{matrix}, 3, 0, 4);
217 gsl_matrix_set($self->{matrix}, 3, 1, 1);
218 gsl_matrix_set($self->{matrix}, 3, 2, 2);
219 gsl_matrix_set($self->{matrix}, 3, 3, 3);
221 my $permutation = gsl_permutation_alloc(4);
222 gsl_permutation_init($permutation);
223 my ($result, $signum) = gsl_linalg_LU_decomp($self->{matrix}, $permutation);
224 ok_similar(gsl_linalg_LU_det($self->{matrix}, $signum), 160);
227 sub GSL_LINALG_LU_LNDET : Tests {
229 map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
231 gsl_matrix_set($self->{matrix}, 1, 0, 2);
232 gsl_matrix_set($self->{matrix}, 1, 1, 3);
233 gsl_matrix_set($self->{matrix}, 1, 2, 4);
234 gsl_matrix_set($self->{matrix}, 1, 3, 1);
236 gsl_matrix_set($self->{matrix}, 2, 0, 3);
237 gsl_matrix_set($self->{matrix}, 2, 1, 4);
238 gsl_matrix_set($self->{matrix}, 2, 2, 1);
239 gsl_matrix_set($self->{matrix}, 2, 3, 2);
241 gsl_matrix_set($self->{matrix}, 3, 0, 4);
242 gsl_matrix_set($self->{matrix}, 3, 1, 1);
243 gsl_matrix_set($self->{matrix}, 3, 2, 2);
244 gsl_matrix_set($self->{matrix}, 3, 3, 3);
246 my $permutation = gsl_permutation_alloc(4);
247 gsl_permutation_init($permutation);
248 gsl_linalg_LU_decomp($self->{matrix}, $permutation);
249 ok_similar(gsl_linalg_LU_lndet($self->{matrix}), log(160));
252 sub GSL_LINALG_QR_DECOMP : Tests {
253 my $matrix = gsl_matrix_alloc(3,5);
255 for($i=0; $i<3; $i++) {
256 for($j=0; $j<5; $j++) {
257 gsl_matrix_set($matrix, $i, $j, 1.0/($i+$j+1.0));
261 my $tau = gsl_vector_alloc(3);
262 my $q = gsl_matrix_alloc(3,3);
263 my $r = gsl_matrix_alloc(3,5);
264 my $a = gsl_matrix_alloc(3,5);
265 my $save = gsl_matrix_alloc(3, 5);
266 gsl_matrix_memcpy($save, $matrix);
268 ok_status(gsl_linalg_QR_decomp($matrix, $tau));
269 ok_similar(gsl_linalg_QR_unpack($matrix, $tau, $q, $r), 0);
271 gsl_blas_dgemm ($CblasNoTrans, $CblasNoTrans, 1.0, $q, $r, 0.0, $a);
274 for($i=0; $i<3; $i++) {
275 for($j=0; $j<5; $j++) {
276 $aij = gsl_matrix_get($a, $i, $j);
277 $mij = gsl_matrix_get($save, $i, $j);
278 ok(is_similar_relative($aij, $mij, 2 * 8.0 * $GSL_DBL_EPSILON));
284 sub GSL_LINALG_CHOLESKY_DECOMP : Tests {
286 map{ gsl_matrix_set($self->{matrix}, 0, $_, $_+1)} (0..3);
288 gsl_matrix_set($self->{matrix}, 1, 0, 2);
289 gsl_matrix_set($self->{matrix}, 1, 1, 5);
290 gsl_matrix_set($self->{matrix}, 1, 2, 8);
291 gsl_matrix_set($self->{matrix}, 1, 3, 11);
293 gsl_matrix_set($self->{matrix}, 2, 0, 3);
294 gsl_matrix_set($self->{matrix}, 2, 1, 8);
295 gsl_matrix_set($self->{matrix}, 2, 2, 14);
296 gsl_matrix_set($self->{matrix}, 2, 3, 20);
298 gsl_matrix_set($self->{matrix}, 3, 0, 4);
299 gsl_matrix_set($self->{matrix}, 3, 1, 11);
300 gsl_matrix_set($self->{matrix}, 3, 2, 20);
301 gsl_matrix_set($self->{matrix}, 3, 3, 30);
303 ok_status(gsl_linalg_cholesky_decomp($self->{matrix}));
304 my $v = gsl_matrix_diagonal($self->{matrix});
306 [ map { gsl_vector_get($v->{vector}, $_)} (0..3) ], [(1)x 4 ]
308 is(gsl_matrix_get($self->{matrix}, 1, 0), 2);
309 is(gsl_matrix_get($self->{matrix}, 2, 0), 3);
310 is(gsl_matrix_get($self->{matrix}, 2, 1), 2);
311 is(gsl_matrix_get($self->{matrix}, 3, 0), 4);
312 is(gsl_matrix_get($self->{matrix}, 3, 1), 3);
313 is(gsl_matrix_get($self->{matrix}, 3, 2), 2);
315 sub GSL_LINALG_HESSENBERG_DECOMP_UNPACK_UNPACK_ACCUM_SET_ZERO : Tests {
318 gsl_matrix_set($self->{matrix}, 1, 0, 3);
319 gsl_matrix_set($self->{matrix}, 1, 1, 2);
320 gsl_matrix_set($self->{matrix}, 1, 2, 5);
321 gsl_matrix_set($self->{matrix}, 1, 3, -4);
323 gsl_matrix_set($self->{matrix}, 1, 0, 5);
324 gsl_matrix_set($self->{matrix}, 1, 1, -3);
325 gsl_matrix_set($self->{matrix}, 1, 2, 6);
326 gsl_matrix_set($self->{matrix}, 1, 3, 9);
328 gsl_matrix_set($self->{matrix}, 2, 0, -2);
329 gsl_matrix_set($self->{matrix}, 2, 1, 1);
330 gsl_matrix_set($self->{matrix}, 2, 2, 5);
331 gsl_matrix_set($self->{matrix}, 2, 3, 8);
333 gsl_matrix_set($self->{matrix}, 3, 0, -6);
334 gsl_matrix_set($self->{matrix}, 3, 1, 7);
335 gsl_matrix_set($self->{matrix}, 3, 2, 2);
336 gsl_matrix_set($self->{matrix}, 3, 3, -8);
337 my $tau = gsl_vector_alloc(4);
338 ok_status(gsl_linalg_hessenberg_decomp($self->{matrix}, $tau));
339 my $U = gsl_matrix_alloc(4,4);
340 ok_status(gsl_linalg_hessenberg_unpack($self->{matrix}, $tau, $U));
341 ok_similar(gsl_matrix_get($U, 0, 0), 1);
343 ok_similar( [ map { gsl_matrix_get($U, $_, 0), 0 } (1..3) ], [ 0, 0, 0, 0, 0, 0 ]);
345 ok_similar( [ map { gsl_matrix_get($U, 0, $_), 0 } (1..3) ], [ 0, 0, 0, 0, 0, 0 ]);
347 ok_similar(gsl_matrix_get($U, 1, 1), -0.620173672946042309);
348 ok_similar(gsl_matrix_get($U, 1, 2), -0.268847804615518438);
349 ok_similar(gsl_matrix_get($U, 1, 3), 0.736956900597335762);
350 ok_similar(gsl_matrix_get($U, 2, 1), 0.248069469178416908);
351 ok_similar(gsl_matrix_get($U, 2, 2), -0.958442423454322956);
352 ok_similar(gsl_matrix_get($U, 2, 3), -0.140888819231843737);
353 ok_similar(gsl_matrix_get($U, 3, 1), 0.744208407535250749);
354 ok_similar(gsl_matrix_get($U, 3, 2), 0.0954409706385089263);
355 ok_similar(gsl_matrix_get($U, 3, 3), 0.661093690241727594);
357 my $V = gsl_matrix_alloc(4,4);
358 ok_status( gsl_linalg_hessenberg_unpack_accum($self->{matrix}, $tau, $V), 0); #I don't know how to test the result of this function...
360 ok_status( gsl_linalg_hessenberg_set_zero($self->{matrix}) );
361 for(my $line = 2; $line<4; $line++) {
362 map { ok_similar(gsl_matrix_get($self->{matrix}, $line, $_), 0, "Set zero", 1e-8) } (0..$line-2);
366 sub GSL_LINALG_BIDIAG_DECOMP_UNPACK_UNPACK2_UNPACK_B : Tests {
368 gsl_matrix_set($self->{matrix}, 1, 0, 3);
369 gsl_matrix_set($self->{matrix}, 1, 1, 2);
370 gsl_matrix_set($self->{matrix}, 1, 2, 5);
371 gsl_matrix_set($self->{matrix}, 1, 3, -4);
373 gsl_matrix_set($self->{matrix}, 1, 0, 5);
374 gsl_matrix_set($self->{matrix}, 1, 1, -3);
375 gsl_matrix_set($self->{matrix}, 1, 2, 6);
376 gsl_matrix_set($self->{matrix}, 1, 3, 9);
378 gsl_matrix_set($self->{matrix}, 2, 0, -2);
379 gsl_matrix_set($self->{matrix}, 2, 1, 1);
380 gsl_matrix_set($self->{matrix}, 2, 2, 5);
381 gsl_matrix_set($self->{matrix}, 2, 3, 8);
383 gsl_matrix_set($self->{matrix}, 3, 0, -6);
384 gsl_matrix_set($self->{matrix}, 3, 1, 7);
385 gsl_matrix_set($self->{matrix}, 3, 2, 2);
386 gsl_matrix_set($self->{matrix}, 3, 3, -8);
387 my $tau_U = gsl_vector_alloc(4);
388 my $tau_V = gsl_vector_alloc(3);
390 ok_status(gsl_linalg_bidiag_decomp($self->{matrix}, $tau_U, $tau_V));
391 my $U = gsl_matrix_alloc(4,4);
392 my $V = gsl_matrix_alloc(4,4);
393 my $diag = gsl_vector_alloc(4);
394 my $superdiag = gsl_vector_alloc(3);
395 ok_status(gsl_linalg_bidiag_unpack($self->{matrix}, $tau_U, $U, $tau_V, $V, $diag, $superdiag));
396 ok_similar(gsl_matrix_get($V, 0, 0), 1);
397 ok_similar( [ map { gsl_matrix_get($V, $_, 0) } (1..3) ], [ (0) x 3 ] );
398 ok_similar( [ map { gsl_matrix_get($V, 0, $_) } (1..3) ], [ (0) x 3 ] );
400 local $TODO = 'look into this';
401 ok_similar(gsl_matrix_get($U, 1, 1), -0.609437002705849772);
402 ok_similar(gsl_matrix_get($U, 1, 2), -0.758604748961341558); #doesn't fit the data I've got...
405 Test::Class->runtests;