1 package Math
::GSL
::Linalg
::Test
;
2 use Math
::GSL
::Test qw
/:all/;
3 use base
q{Test::Class};
5 use Math
::GSL
::Linalg qw
/:all/;
6 use Math
::GSL
::Matrix qw
/:all/;
7 use Math
::GSL
::Permutation qw
/:all/;
8 use Math
::GSL
::Vector qw
/:all/;
9 use Math
::GSL
::CBLAS qw
/:all/;
10 use Math
::GSL
::BLAS qw
/:all/;
11 use Math
::GSL
::Machine qw
/:all/;
12 use Math
::GSL
::Complex qw
/:all/;
13 use Math
::GSL qw
/:all/;
15 use Math
::GSL
::Errno qw
/:all/;
18 sub make_fixture
: Test
(setup
) {
20 $self->{matrix
} = gsl_matrix_alloc
(4, 4);
23 sub teardown
: Test
(teardown
) {
24 unlink 'linalg' if -f
'linalg';
27 sub GSL_LINALG_LU_DECOMP
: Tests
{
28 my $base = Math
::GSL
::Matrix
->new(4,4);
29 $base->set_row(0, [0,1,2,3])
30 ->set_row(1, [5,6,7,8])
31 ->set_row(2, [9,10,11,12])
32 ->set_row(3, [13,14,15,16]);
33 my $permutation = gsl_permutation_alloc
(4);
34 gsl_permutation_init
($permutation);
35 my $first = Math
::GSL
::Matrix
->new(4,4);
36 gsl_matrix_memcpy
($first->raw, $base->raw);
38 my ($result, $signum) = gsl_linalg_LU_decomp
($base->raw, $permutation);
39 is_deeply
( [ $result, $signum ], [ 0, 1] );
41 my $U = Math
::GSL
::Matrix
->new(4,4);
42 my $R = Math
::GSL
::Matrix
->new(4,4);
43 my $L = Math
::GSL
::Matrix
->new(4,4);
44 gsl_matrix_set_identity
($L->raw);
46 for ($line=3; $line>-1; $line--) {
47 map { gsl_matrix_set
($U->raw, $line, $_, gsl_matrix_get
($base->raw, $line, $_)) } ($line..3) };
48 for ($line=3; $line>1; $line--) {
49 map { gsl_matrix_set
($L->raw, $line, $_, gsl_matrix_get
($base->raw, $line, $_)) } (0..$line-1) };
50 gsl_blas_dgemm
($CblasNoTrans, $CblasNoTrans, 1, $L->raw, $U->raw, 1, $R->raw);
51 my @permutations = map { gsl_permutation_get
($permutation, $_) } (0..3);
53 my @got = $first->row($permutations[$_])->as_list;
54 my @results = $R->row($_)->as_list;
55 ok_similar
([@got], [@results], "resulting row $_"); } (0..3);
58 sub GSL_LINALG_LU_SOLVE
: Tests
{
60 gsl_matrix_set
($self->{matrix
}, 0, 0, 1);
61 gsl_matrix_set
($self->{matrix
}, 0, 1, 1);
62 gsl_matrix_set
($self->{matrix
}, 0, 2, 2);
63 gsl_matrix_set
($self->{matrix
}, 0, 3, 1);
65 gsl_matrix_set
($self->{matrix
}, 1, 0, 2);
66 gsl_matrix_set
($self->{matrix
}, 1, 1, 3);
67 gsl_matrix_set
($self->{matrix
}, 1, 2, -1);
68 gsl_matrix_set
($self->{matrix
}, 1, 3, 2);
70 gsl_matrix_set
($self->{matrix
}, 2, 0, 5);
71 gsl_matrix_set
($self->{matrix
}, 2, 1, -1);
72 gsl_matrix_set
($self->{matrix
}, 2, 2, 1);
73 gsl_matrix_set
($self->{matrix
}, 2, 3, -1);
75 gsl_matrix_set
($self->{matrix
}, 3, 0, 1);
76 gsl_matrix_set
($self->{matrix
}, 3, 1, 0);
77 gsl_matrix_set
($self->{matrix
}, 3, 2, 7);
78 gsl_matrix_set
($self->{matrix
}, 3, 3, 1);
80 my $b = gsl_vector_alloc
(4);
81 gsl_vector_set
($b, 0, 4);
82 gsl_vector_set
($b, 1, 1);
83 gsl_vector_set
($b, 2, 2);
84 gsl_vector_set
($b, 3, 11);
86 my $x = gsl_vector_alloc
(4);
88 my $permutation = gsl_permutation_alloc
(4);
89 gsl_permutation_init
($permutation);
90 gsl_linalg_LU_decomp
($self->{matrix
}, $permutation);
91 gsl_linalg_LU_solve
($self->{matrix
}, $permutation, $b, $x);
92 my $value = gsl_vector_get
($x, 0);
94 [ map { gsl_vector_get
($x, $_) } (1..3) ],
95 [ 3-10*$value, -2*$value+2, 13*$value-3 ]
99 sub GSL_LINALG_LU_SVX
: Tests
{
101 gsl_matrix_set
($self->{matrix
}, 0, 0, 1);
102 gsl_matrix_set
($self->{matrix
}, 0, 1, 1);
103 gsl_matrix_set
($self->{matrix
}, 0, 2, 2);
104 gsl_matrix_set
($self->{matrix
}, 0, 3, 1);
106 gsl_matrix_set
($self->{matrix
}, 1, 0, 2);
107 gsl_matrix_set
($self->{matrix
}, 1, 1, 3);
108 gsl_matrix_set
($self->{matrix
}, 1, 2, -1);
109 gsl_matrix_set
($self->{matrix
}, 1, 3, 2);
111 gsl_matrix_set
($self->{matrix
}, 2, 0, 5);
112 gsl_matrix_set
($self->{matrix
}, 2, 1, -1);
113 gsl_matrix_set
($self->{matrix
}, 2, 2, 1);
114 gsl_matrix_set
($self->{matrix
}, 2, 3, -1);
116 gsl_matrix_set
($self->{matrix
}, 3, 0, 1);
117 gsl_matrix_set
($self->{matrix
}, 3, 1, 0);
118 gsl_matrix_set
($self->{matrix
}, 3, 2, 7);
119 gsl_matrix_set
($self->{matrix
}, 3, 3, 1);
121 my $x = gsl_vector_alloc
(4);
122 gsl_vector_set
($x, 0, 4);
123 gsl_vector_set
($x, 1, 1);
124 gsl_vector_set
($x, 2, 2);
125 gsl_vector_set
($x, 3, 11);
127 my $permutation = gsl_permutation_alloc
(4);
128 gsl_permutation_init
($permutation);
129 gsl_linalg_LU_decomp
($self->{matrix
}, $permutation);
130 gsl_linalg_LU_svx
($self->{matrix
}, $permutation, $x);
131 my $value = gsl_vector_get
($x, 0);
133 [ map { gsl_vector_get
($x, $_) } (1..3) ],
134 [ 3-10*$value, -2*$value+2, 13*$value-3 ]
138 sub GSL_LINALG_LU_INVERT
: Tests
{
140 map { gsl_matrix_set
($self->{matrix
}, 0, $_, $_+1) } (0..3);
142 gsl_matrix_set
($self->{matrix
}, 1, 0, 2);
143 gsl_matrix_set
($self->{matrix
}, 1, 1, 3);
144 gsl_matrix_set
($self->{matrix
}, 1, 2, 4);
145 gsl_matrix_set
($self->{matrix
}, 1, 3, 1);
147 gsl_matrix_set
($self->{matrix
}, 2, 0, 3);
148 gsl_matrix_set
($self->{matrix
}, 2, 1, 4);
149 gsl_matrix_set
($self->{matrix
}, 2, 2, 1);
150 gsl_matrix_set
($self->{matrix
}, 2, 3, 2);
152 gsl_matrix_set
($self->{matrix
}, 3, 0, 4);
153 gsl_matrix_set
($self->{matrix
}, 3, 1, 1);
154 gsl_matrix_set
($self->{matrix
}, 3, 2, 2);
155 gsl_matrix_set
($self->{matrix
}, 3, 3, 3);
157 my $inverse = gsl_matrix_alloc
(4,4);
158 my $permutation = gsl_permutation_alloc
(4);
159 gsl_permutation_init
($permutation);
160 gsl_linalg_LU_decomp
($self->{matrix
}, $permutation);
161 gsl_linalg_LU_invert
($self->{matrix
}, $permutation, $inverse);
163 is_similar
(gsl_matrix_get
($inverse, 0, 0), -9/40);
164 is_similar
(gsl_matrix_get
($inverse, 0, 1), 1/40);
165 is_similar
(gsl_matrix_get
($inverse, 0, 2), 1/40);
166 is_similar
(gsl_matrix_get
($inverse, 0, 3), 11/40);
168 is_similar
(gsl_matrix_get
($inverse, 1, 0), 1/40);
169 is_similar
(gsl_matrix_get
($inverse, 1, 1), 1/40);
170 is_similar
(gsl_matrix_get
($inverse, 1, 2), 11/40);
171 is_similar
(gsl_matrix_get
($inverse, 1, 3), -9/40);
173 is_similar
(gsl_matrix_get
($inverse, 2, 0), 1/40);
174 is_similar
(gsl_matrix_get
($inverse, 2, 1), 11/40);
175 is_similar
(gsl_matrix_get
($inverse, 2, 2), -9/40);
176 is_similar
(gsl_matrix_get
($inverse, 2, 3), 1/40);
178 is_similar
(gsl_matrix_get
($inverse, 3, 0), 11/40);
179 is_similar
(gsl_matrix_get
($inverse, 3, 1), -9/40);
180 is_similar
(gsl_matrix_get
($inverse, 3, 2), 1/40);
181 is_similar
(gsl_matrix_get
($inverse, 3, 3), 1/40);
184 sub GSL_LINALG_LU_DET
: Tests
{
186 map { gsl_matrix_set
($self->{matrix
}, 0, $_, $_+1) } (0..3);
188 gsl_matrix_set
($self->{matrix
}, 1, 0, 2);
189 gsl_matrix_set
($self->{matrix
}, 1, 1, 3);
190 gsl_matrix_set
($self->{matrix
}, 1, 2, 4);
191 gsl_matrix_set
($self->{matrix
}, 1, 3, 1);
193 gsl_matrix_set
($self->{matrix
}, 2, 0, 3);
194 gsl_matrix_set
($self->{matrix
}, 2, 1, 4);
195 gsl_matrix_set
($self->{matrix
}, 2, 2, 1);
196 gsl_matrix_set
($self->{matrix
}, 2, 3, 2);
198 gsl_matrix_set
($self->{matrix
}, 3, 0, 4);
199 gsl_matrix_set
($self->{matrix
}, 3, 1, 1);
200 gsl_matrix_set
($self->{matrix
}, 3, 2, 2);
201 gsl_matrix_set
($self->{matrix
}, 3, 3, 3);
203 my $permutation = gsl_permutation_alloc
(4);
204 gsl_permutation_init
($permutation);
205 my ($result, $signum) = gsl_linalg_LU_decomp
($self->{matrix
}, $permutation);
206 ok_similar
(gsl_linalg_LU_det
($self->{matrix
}, $signum), 160);
209 sub GSL_LINALG_LU_LNDET
: Tests
{
211 map { gsl_matrix_set
($self->{matrix
}, 0, $_, $_+1) } (0..3);
213 gsl_matrix_set
($self->{matrix
}, 1, 0, 2);
214 gsl_matrix_set
($self->{matrix
}, 1, 1, 3);
215 gsl_matrix_set
($self->{matrix
}, 1, 2, 4);
216 gsl_matrix_set
($self->{matrix
}, 1, 3, 1);
218 gsl_matrix_set
($self->{matrix
}, 2, 0, 3);
219 gsl_matrix_set
($self->{matrix
}, 2, 1, 4);
220 gsl_matrix_set
($self->{matrix
}, 2, 2, 1);
221 gsl_matrix_set
($self->{matrix
}, 2, 3, 2);
223 gsl_matrix_set
($self->{matrix
}, 3, 0, 4);
224 gsl_matrix_set
($self->{matrix
}, 3, 1, 1);
225 gsl_matrix_set
($self->{matrix
}, 3, 2, 2);
226 gsl_matrix_set
($self->{matrix
}, 3, 3, 3);
228 my $permutation = gsl_permutation_alloc
(4);
229 gsl_permutation_init
($permutation);
230 gsl_linalg_LU_decomp
($self->{matrix
}, $permutation);
231 ok_similar
(gsl_linalg_LU_lndet
($self->{matrix
}), log(160));
234 sub GSL_LINALG_QR_DECOMP
: Tests
{
235 my $matrix = gsl_matrix_alloc
(3,5);
237 for($i=0; $i<3; $i++) {
238 for($j=0; $j<5; $j++) {
239 gsl_matrix_set
($matrix, $i, $j, 1.0/($i+$j+1.0));
243 my $tau = gsl_vector_alloc
(3);
244 my $q = gsl_matrix_alloc
(3,3);
245 my $r = gsl_matrix_alloc
(3,5);
246 my $a = gsl_matrix_alloc
(3,5);
247 my $save = gsl_matrix_alloc
(3, 5);
248 gsl_matrix_memcpy
($save, $matrix);
250 ok_status
(gsl_linalg_QR_decomp
($matrix, $tau));
251 is
(gsl_linalg_QR_unpack
($matrix, $tau, $q, $r), 0);
253 gsl_blas_dgemm
($CblasNoTrans, $CblasNoTrans, 1.0, $q, $r, 0.0, $a);
256 for($i=0; $i<3; $i++) {
257 for($j=0; $j<5; $j++) {
258 $aij = gsl_matrix_get
($a, $i, $j);
259 $mij = gsl_matrix_get
($save, $i, $j);
260 ok
(is_similar_relative
($aij, $mij, 2 * 8.0 * $GSL_DBL_EPSILON));
266 sub GSL_LINALG_CHOLESKY_DECOMP
: Tests
{
268 map{ gsl_matrix_set
($self->{matrix
}, 0, $_, $_+1)} (0..3);
270 gsl_matrix_set
($self->{matrix
}, 1, 0, 2);
271 gsl_matrix_set
($self->{matrix
}, 1, 1, 5);
272 gsl_matrix_set
($self->{matrix
}, 1, 2, 8);
273 gsl_matrix_set
($self->{matrix
}, 1, 3, 11);
275 gsl_matrix_set
($self->{matrix
}, 2, 0, 3);
276 gsl_matrix_set
($self->{matrix
}, 2, 1, 8);
277 gsl_matrix_set
($self->{matrix
}, 2, 2, 14);
278 gsl_matrix_set
($self->{matrix
}, 2, 3, 20);
280 gsl_matrix_set
($self->{matrix
}, 3, 0, 4);
281 gsl_matrix_set
($self->{matrix
}, 3, 1, 11);
282 gsl_matrix_set
($self->{matrix
}, 3, 2, 20);
283 gsl_matrix_set
($self->{matrix
}, 3, 3, 30);
285 ok_status
(gsl_linalg_cholesky_decomp
($self->{matrix
}));
286 my $v = gsl_matrix_diagonal
($self->{matrix
});
288 [ map { gsl_vector_get
($v->{vector
}, $_)} (0..3) ], [(1)x
4 ]
290 is
(gsl_matrix_get
($self->{matrix
}, 1, 0), 2);
291 is
(gsl_matrix_get
($self->{matrix
}, 2, 0), 3);
292 is
(gsl_matrix_get
($self->{matrix
}, 2, 1), 2);
293 is
(gsl_matrix_get
($self->{matrix
}, 3, 0), 4);
294 is
(gsl_matrix_get
($self->{matrix
}, 3, 1), 3);
295 is
(gsl_matrix_get
($self->{matrix
}, 3, 2), 2);
297 sub GSL_LINALG_HESSENBERG_DECOMP_UNPACK_UNPACK_ACCUM_SET_ZERO
: Tests
{
300 gsl_matrix_set
($self->{matrix
}, 1, 0, 3);
301 gsl_matrix_set
($self->{matrix
}, 1, 1, 2);
302 gsl_matrix_set
($self->{matrix
}, 1, 2, 5);
303 gsl_matrix_set
($self->{matrix
}, 1, 3, -4);
305 gsl_matrix_set
($self->{matrix
}, 1, 0, 5);
306 gsl_matrix_set
($self->{matrix
}, 1, 1, -3);
307 gsl_matrix_set
($self->{matrix
}, 1, 2, 6);
308 gsl_matrix_set
($self->{matrix
}, 1, 3, 9);
310 gsl_matrix_set
($self->{matrix
}, 2, 0, -2);
311 gsl_matrix_set
($self->{matrix
}, 2, 1, 1);
312 gsl_matrix_set
($self->{matrix
}, 2, 2, 5);
313 gsl_matrix_set
($self->{matrix
}, 2, 3, 8);
315 gsl_matrix_set
($self->{matrix
}, 3, 0, -6);
316 gsl_matrix_set
($self->{matrix
}, 3, 1, 7);
317 gsl_matrix_set
($self->{matrix
}, 3, 2, 2);
318 gsl_matrix_set
($self->{matrix
}, 3, 3, -8);
319 my $tau = gsl_vector_alloc
(4);
320 ok_status
(gsl_linalg_hessenberg_decomp
($self->{matrix
}, $tau));
321 my $U = gsl_matrix_alloc
(4,4);
322 ok_status
(gsl_linalg_hessenberg_unpack
($self->{matrix
}, $tau, $U));
323 is
(gsl_matrix_get
($U, 0, 0), 1);
324 map { is
(gsl_matrix_get
($U, $_, 0), 0) } (1..3);
325 map { is
(gsl_matrix_get
($U, 0, $_), 0) } (1..3);
326 is_similar
(gsl_matrix_get
($U, 1, 1), -0.620173672946042309);
327 is_similar
(gsl_matrix_get
($U, 1, 2), -0.268847804615518438);
328 is_similar
(gsl_matrix_get
($U, 1, 3), 0.736956900597335762);
329 is_similar
(gsl_matrix_get
($U, 2, 1), 0.248069469178416908);
330 is_similar
(gsl_matrix_get
($U, 2, 2), -0.958442423454322956);
331 is_similar
(gsl_matrix_get
($U, 2, 3), -0.140888819231843737);
332 is_similar
(gsl_matrix_get
($U, 3, 1), 0.744208407535250749);
333 is_similar
(gsl_matrix_get
($U, 3, 2), 0.0954409706385089263);
334 is_similar
(gsl_matrix_get
($U, 3, 3), 0.661093690241727594);
336 my $V = gsl_matrix_alloc
(4,4);
337 is
(gsl_linalg_hessenberg_unpack_accum
($self->{matrix
}, $tau, $V), 0); #I don't know how to test the result of this function...
339 is
(gsl_linalg_hessenberg_set_zero
($self->{matrix
}), 0);
340 for(my $line = 2; $line<4; $line++) {
341 map { is
(gsl_matrix_get
($self->{matrix
}, $line, $_), 0, "Set zero") } (0..$line-2);
345 sub GSL_LINALG_BIDIAG_DECOMP_UNPACK_UNPACK2_UNPACK_B
: Tests
{
347 gsl_matrix_set
($self->{matrix
}, 1, 0, 3);
348 gsl_matrix_set
($self->{matrix
}, 1, 1, 2);
349 gsl_matrix_set
($self->{matrix
}, 1, 2, 5);
350 gsl_matrix_set
($self->{matrix
}, 1, 3, -4);
352 gsl_matrix_set
($self->{matrix
}, 1, 0, 5);
353 gsl_matrix_set
($self->{matrix
}, 1, 1, -3);
354 gsl_matrix_set
($self->{matrix
}, 1, 2, 6);
355 gsl_matrix_set
($self->{matrix
}, 1, 3, 9);
357 gsl_matrix_set
($self->{matrix
}, 2, 0, -2);
358 gsl_matrix_set
($self->{matrix
}, 2, 1, 1);
359 gsl_matrix_set
($self->{matrix
}, 2, 2, 5);
360 gsl_matrix_set
($self->{matrix
}, 2, 3, 8);
362 gsl_matrix_set
($self->{matrix
}, 3, 0, -6);
363 gsl_matrix_set
($self->{matrix
}, 3, 1, 7);
364 gsl_matrix_set
($self->{matrix
}, 3, 2, 2);
365 gsl_matrix_set
($self->{matrix
}, 3, 3, -8);
366 my $tau_U = gsl_vector_alloc
(4);
367 my $tau_V = gsl_vector_alloc
(3);
369 ok_status
(gsl_linalg_bidiag_decomp
($self->{matrix
}, $tau_U, $tau_V));
370 my $U = gsl_matrix_alloc
(4,4);
371 my $V = gsl_matrix_alloc
(4,4);
372 my $diag = gsl_vector_alloc
(4);
373 my $superdiag = gsl_vector_alloc
(3);
374 ok_status
(gsl_linalg_bidiag_unpack
($self->{matrix
}, $tau_U, $U, $tau_V, $V, $diag, $superdiag));
375 is
(gsl_matrix_get
($V, 0, 0), 1);
376 ok_similar
( [ map { gsl_matrix_get
($V, $_, 0) } (1..3) ], [ (0) x
3 ] );
377 ok_similar
( [ map { gsl_matrix_get
($V, 0, $_) } (1..3) ], [ (0) x
3 ] );
379 is_similar
(gsl_matrix_get
($U, 1, 1), -0.609437002705849772);
380 is_similar
(gsl_matrix_get
($U, 1, 2), -0.758604748961341558); #doesn't fit the data I've got...