lib/Math/GSL/Linalg/Test.pm

   1 package Math::GSL::Linalg::Test;
   2 use base q{Test::Class};
   3 use Test::More;
   4 use Math::GSL::Linalg qw/:all/;
   5 use Math::GSL::Matrix qw/:all/;
   6 use Math::GSL::Permutation qw/:all/;
   7 use Math::GSL::Vector qw/:all/;
   8 use Math::GSL qw/:all/;
   9 use Data::Dumper;
  10 use Math::GSL::Errno qw/:all/;
  11 use strict;
  12
  13 sub make_fixture : Test(setup) {
  14     my $self = shift;
  15     $self->{matrix} = gsl_matrix_alloc(4, 4);
  16 }
  17
  18 sub teardown : Test(teardown) {
  19     unlink 'linalg' if -f 'linalg';
  20 }
  21
  22 sub GSL_LINALG_LU_DECOMP : Tests {
  23     my $self = shift;
  24     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
  25     map { gsl_matrix_set($self->{matrix}, 1, $_, $_+5) } (0..3);
  26     map { gsl_matrix_set($self->{matrix}, 2, $_, $_+9) } (0..3);
  27     map { gsl_matrix_set($self->{matrix}, 3, $_, $_+13) } (0..3);
  28
  29     my $permutation = gsl_permutation_alloc(4);
  30     gsl_permutation_init($permutation);
  31     my ($result, $signum) = gsl_linalg_LU_decomp($self->{matrix}, $permutation);
  32     is_deeply( [ $result, $signum ], [ 0, 1] );
  33     map { is( gsl_matrix_get($self->{matrix}, 0, $_), $_+1) } (0..3); # I have no idea why these tests fail, I got my values for the LU decompositon from maple and they are valid...
  34     ok_similar( [ map { gsl_matrix_get($self->{matrix}, 0, $_) } (2..3) ],
  35                 [ 0, 0 ]
  36               );
  37     is (gsl_matrix_get($self->{matrix}, 3, 3),0);
  38     is (gsl_matrix_get($self->{matrix}, 1, 0),5);
  39     is (gsl_matrix_get($self->{matrix}, 2, 0),9);
  40     is (gsl_matrix_get($self->{matrix}, 2, 1),2);
  41     is (gsl_matrix_get($self->{matrix}, 3, 0),13);
  42     is (gsl_matrix_get($self->{matrix}, 3, 1),3);
  43     is (gsl_matrix_get($self->{matrix}, 3, 2),0);
  44 }
  45
  46 sub GSL_LINALG_LU_SOLVE : Tests {
  47     my $self = shift;
  48     gsl_matrix_set($self->{matrix}, 0, 0, 1);
  49     gsl_matrix_set($self->{matrix}, 0, 1, 1);
  50     gsl_matrix_set($self->{matrix}, 0, 2, 2);
  51     gsl_matrix_set($self->{matrix}, 0, 3, 1);
  52
  53     gsl_matrix_set($self->{matrix}, 1, 0, 2);
  54     gsl_matrix_set($self->{matrix}, 1, 1, 3);
  55     gsl_matrix_set($self->{matrix}, 1, 2, -1);
  56     gsl_matrix_set($self->{matrix}, 1, 3, 2);
  57
  58     gsl_matrix_set($self->{matrix}, 2, 0, 5);
  59     gsl_matrix_set($self->{matrix}, 2, 1, -1);
  60     gsl_matrix_set($self->{matrix}, 2, 2, 1);
  61     gsl_matrix_set($self->{matrix}, 2, 3, -1);
  62
  63     gsl_matrix_set($self->{matrix}, 3, 0, 1);
  64     gsl_matrix_set($self->{matrix}, 3, 1, 0);
  65     gsl_matrix_set($self->{matrix}, 3, 2, 7);
  66     gsl_matrix_set($self->{matrix}, 3, 3, 1);
  67
  68     my $b = gsl_vector_alloc(4);
  69     gsl_vector_set($b, 0, 4);
  70     gsl_vector_set($b, 1, 1);
  71     gsl_vector_set($b, 2, 2);
  72     gsl_vector_set($b, 3, 11);
  73
  74     my $x = gsl_vector_alloc(4);
  75
  76     my $permutation = gsl_permutation_alloc(4);
  77     gsl_permutation_init($permutation);
  78     gsl_linalg_LU_decomp($self->{matrix}, $permutation);
  79     gsl_linalg_LU_solve($self->{matrix}, $permutation, $b, $x);
  80     my $value = gsl_vector_get($x, 0);
  81     ok_similar(
  82         [ map { gsl_vector_get($x, $_) } (1..3) ],
  83         [ 3-10*$value, -2*$value+2, 13*$value-3 ]
  84     );
  85 }
  86
  87 sub GSL_LINALG_LU_SVX : Tests {
  88     my $self = shift;
  89     gsl_matrix_set($self->{matrix}, 0, 0, 1);
  90     gsl_matrix_set($self->{matrix}, 0, 1, 1);
  91     gsl_matrix_set($self->{matrix}, 0, 2, 2);
  92     gsl_matrix_set($self->{matrix}, 0, 3, 1);
  93
  94     gsl_matrix_set($self->{matrix}, 1, 0, 2);
  95     gsl_matrix_set($self->{matrix}, 1, 1, 3);
  96     gsl_matrix_set($self->{matrix}, 1, 2, -1);
  97     gsl_matrix_set($self->{matrix}, 1, 3, 2);
  98
  99     gsl_matrix_set($self->{matrix}, 2, 0, 5);
 100     gsl_matrix_set($self->{matrix}, 2, 1, -1);
 101     gsl_matrix_set($self->{matrix}, 2, 2, 1);
 102     gsl_matrix_set($self->{matrix}, 2, 3, -1);
 103
 104     gsl_matrix_set($self->{matrix}, 3, 0, 1);
 105     gsl_matrix_set($self->{matrix}, 3, 1, 0);
 106     gsl_matrix_set($self->{matrix}, 3, 2, 7);
 107     gsl_matrix_set($self->{matrix}, 3, 3, 1);
 108
 109     my $x = gsl_vector_alloc(4);
 110     gsl_vector_set($x, 0, 4);
 111     gsl_vector_set($x, 1, 1);
 112     gsl_vector_set($x, 2, 2);
 113     gsl_vector_set($x, 3, 11);
 114
 115     my $permutation = gsl_permutation_alloc(4);
 116     gsl_permutation_init($permutation);
 117     gsl_linalg_LU_decomp($self->{matrix}, $permutation);
 118     gsl_linalg_LU_svx($self->{matrix}, $permutation, $x);
 119     my $value = gsl_vector_get($x, 0);
 120     ok_similar(
 121         [ map { gsl_vector_get($x, $_) } (1..3) ],
 122         [ 3-10*$value, -2*$value+2, 13*$value-3 ]
 123     );
 124 }
 125
 126 sub GSL_LINALG_LU_INVERT : Tests {
 127     my $self = shift;
 128     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 129
 130     gsl_matrix_set($self->{matrix}, 1, 0, 2);
 131     gsl_matrix_set($self->{matrix}, 1, 1, 3);
 132     gsl_matrix_set($self->{matrix}, 1, 2, 4);
 133     gsl_matrix_set($self->{matrix}, 1, 3, 1);
 134
 135     gsl_matrix_set($self->{matrix}, 2, 0, 3);
 136     gsl_matrix_set($self->{matrix}, 2, 1, 4);
 137     gsl_matrix_set($self->{matrix}, 2, 2, 1);
 138     gsl_matrix_set($self->{matrix}, 2, 3, 2);
 139
 140     gsl_matrix_set($self->{matrix}, 3, 0, 4);
 141     gsl_matrix_set($self->{matrix}, 3, 1, 1);
 142     gsl_matrix_set($self->{matrix}, 3, 2, 2);
 143     gsl_matrix_set($self->{matrix}, 3, 3, 3);
 144
 145     my $inverse = gsl_matrix_alloc(4,4);
 146     my $permutation = gsl_permutation_alloc(4);
 147     gsl_permutation_init($permutation);
 148     gsl_linalg_LU_decomp($self->{matrix}, $permutation);
 149     gsl_linalg_LU_invert($self->{matrix}, $permutation, $inverse);
 150
 151     is (gsl_matrix_get($inverse, 0, 0), -9/40);
 152     is (gsl_matrix_get($inverse, 0, 1), 1/40);
 153     is (gsl_matrix_get($inverse, 0, 2), 1/40);
 154     is (gsl_matrix_get($inverse, 0, 3), 11/40);
 155
 156     is (gsl_matrix_get($inverse, 1, 0), 1/40);
 157     is (gsl_matrix_get($inverse, 1, 1), 1/40);
 158     is (gsl_matrix_get($inverse, 1, 2), 11/40);
 159     is (gsl_matrix_get($inverse, 1, 3), -9/40);
 160
 161     is (gsl_matrix_get($inverse, 2, 0), 1/40);
 162     is (gsl_matrix_get($inverse, 2, 1), 11/40);
 163     is (gsl_matrix_get($inverse, 2, 2), -9/40);
 164     is (gsl_matrix_get($inverse, 2, 3), 1/40);
 165
 166     is (gsl_matrix_get($inverse, 3, 0), 11/40);
 167     is (gsl_matrix_get($inverse, 3, 1), -9/40);
 168     is (gsl_matrix_get($inverse, 3, 2), 1/40);
 169     is (gsl_matrix_get($inverse, 3, 3), 1/40);
 170 }
 171
 172 sub GSL_LINALG_LU_DET : Tests {
 173     my $self = shift;
 174     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 175
 176     gsl_matrix_set($self->{matrix}, 1, 0, 2);
 177     gsl_matrix_set($self->{matrix}, 1, 1, 3);
 178     gsl_matrix_set($self->{matrix}, 1, 2, 4);
 179     gsl_matrix_set($self->{matrix}, 1, 3, 1);
 180
 181     gsl_matrix_set($self->{matrix}, 2, 0, 3);
 182     gsl_matrix_set($self->{matrix}, 2, 1, 4);
 183     gsl_matrix_set($self->{matrix}, 2, 2, 1);
 184     gsl_matrix_set($self->{matrix}, 2, 3, 2);
 185
 186     gsl_matrix_set($self->{matrix}, 3, 0, 4);
 187     gsl_matrix_set($self->{matrix}, 3, 1, 1);
 188     gsl_matrix_set($self->{matrix}, 3, 2, 2);
 189     gsl_matrix_set($self->{matrix}, 3, 3, 3);
 190
 191     my $permutation = gsl_permutation_alloc(4);
 192     gsl_permutation_init($permutation);
 193     my ($result, $signum) = gsl_linalg_LU_decomp($self->{matrix}, $permutation);
 194     ok_similar(gsl_linalg_LU_det($self->{matrix}, $signum), 160);
 195 }
 196
 197 sub GSL_LINALG_LU_LNDET : Tests {
 198     my $self = shift;
 199     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 200
 201     gsl_matrix_set($self->{matrix}, 1, 0, 2);
 202     gsl_matrix_set($self->{matrix}, 1, 1, 3);
 203     gsl_matrix_set($self->{matrix}, 1, 2, 4);
 204     gsl_matrix_set($self->{matrix}, 1, 3, 1);
 205
 206     gsl_matrix_set($self->{matrix}, 2, 0, 3);
 207     gsl_matrix_set($self->{matrix}, 2, 1, 4);
 208     gsl_matrix_set($self->{matrix}, 2, 2, 1);
 209     gsl_matrix_set($self->{matrix}, 2, 3, 2);
 210
 211     gsl_matrix_set($self->{matrix}, 3, 0, 4);
 212     gsl_matrix_set($self->{matrix}, 3, 1, 1);
 213     gsl_matrix_set($self->{matrix}, 3, 2, 2);
 214     gsl_matrix_set($self->{matrix}, 3, 3, 3);
 215
 216     my $permutation = gsl_permutation_alloc(4);
 217     gsl_permutation_init($permutation);
 218     gsl_linalg_LU_decomp($self->{matrix}, $permutation);
 219     ok_similar(gsl_linalg_LU_lndet($self->{matrix}), log(160));
 220 }
 221
 222 sub GSL_LINALG_QR_DECOMP : Tests {
 223     local $TODO ="the values doesn't seem to fit the value I got from maple. Probably the same problem than gsl_linalg_LU_decomp...";
 224     my $matrix = gsl_matrix_alloc(4,3);
 225     gsl_matrix_set($matrix, 0, 0, -3);
 226     gsl_matrix_set($matrix, 1, 0, 2);
 227     gsl_matrix_set($matrix, 2, 0, -5);
 228     gsl_matrix_set($matrix, 3, 0, 1);
 229
 230     gsl_matrix_set($matrix, 0, 1, 2);
 231     gsl_matrix_set($matrix, 1, 1, 1);
 232     gsl_matrix_set($matrix, 2, 1, 2);
 233     gsl_matrix_set($matrix, 3, 1, -3);
 234
 235     gsl_matrix_set($matrix, 0, 2, 4);
 236     gsl_matrix_set($matrix, 1, 2, -1);
 237     gsl_matrix_set($matrix, 2, 2, 4);
 238     gsl_matrix_set($matrix, 3, 2, 2);
 239
 240     my $tau = gsl_vector_alloc(3);
 241     gsl_linalg_QR_decomp($matrix, $tau);
 242
 243     is(gsl_matrix_get($matrix, 0, 0), sqrt(29));
 244     is(gsl_matrix_get($matrix, 1, 0), (-8/29)*sqrt(29));
 245     is(gsl_matrix_get($matrix, 2, 0), (35/29)*sqrt(29));
 246     is(gsl_matrix_get($matrix, 3, 0), (-1/29)*sqrt(29));
 247 }
 248 1;