t/Linalg.t

   1 package Math::GSL::Linalg::Test;
   2 use base q{Test::Class};
   3 use Test::More tests => 53;
   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::Vector      qw/:all/;
  12 use Math::GSL::Machine     qw/:all/;
  13 use Math::GSL::Complex     qw/:all/;
  14 use Math::GSL::Permutation qw/:all/;
  15 use Data::Dumper;
  16 use strict;
  17
  18 sub make_fixture : Test(setup) {
  19     my $self = shift;
  20     $self->{matrix} = gsl_matrix_alloc(4, 4);
  21 }
  22
  23 sub teardown : Test(teardown) {
  24     unlink 'linalg' if -f 'linalg';
  25 }
  26
  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);
  37
  38     my ($result, $signum) = gsl_linalg_LU_decomp($base->raw, $permutation);
  39     is_deeply( [ $result, $signum ], [ 0, 1] );
  40
  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);
  45     my $line;
  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);
  52     map {
  53     my @got = $first->row($permutations[$_])->as_list;
  54     my @results = $R->row($_)->as_list;
  55     ok_similar([@got], [@results], "resulting row $_"); } (0..3);
  56 }
  57
  58 sub GSL_LINALG_LU_SOLVE : Tests {
  59     my $self = shift;
  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);
  64
  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);
  69
  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);
  74
  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);
  79
  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);
  85
  86     my $x = gsl_vector_alloc(4);
  87
  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);
  93     ok_similar(
  94         [ map { gsl_vector_get($x, $_) } (1..3) ],
  95         [ 3-10*$value, -2*$value+2, 13*$value-3 ]
  96     );
  97 }
  98
  99 sub GSL_LINALG_LU_SVX : Tests {
 100     my $self = shift;
 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);
 105
 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);
 110
 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);
 115
 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);
 120
 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);
 126
 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);
 132     ok_similar(
 133         [ map { gsl_vector_get($x, $_) } (1..3) ],
 134         [ 3-10*$value, -2*$value+2, 13*$value-3 ]
 135     );
 136 }
 137
 138 sub GSL_LINALG_LU_INVERT : Tests {
 139     my $self = shift;
 140     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 141
 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);
 146
 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);
 151
 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);
 156
 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);
 162
 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);
 167
 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);
 172
 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);
 177
 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);
 182 }
 183
 184 sub GSL_LINALG_LU_DET : Tests {
 185     my $self = shift;
 186     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 187
 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);
 192
 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);
 197
 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);
 202
 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);
 207 }
 208
 209 sub GSL_LINALG_LU_LNDET : Tests {
 210     my $self = shift;
 211     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 212
 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);
 217
 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);
 222
 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);
 227
 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));
 232 }
 233
 234 sub GSL_LINALG_QR_DECOMP : Tests {
 235     my $matrix = gsl_matrix_alloc(3,5);
 236     my ($i, $j);
 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));
 240      }
 241     }
 242
 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);
 249
 250     ok_status(gsl_linalg_QR_decomp($matrix, $tau));
 251     is(gsl_linalg_QR_unpack($matrix, $tau, $q, $r), 0);
 252   # compute a = q r
 253   gsl_blas_dgemm ($CblasNoTrans, $CblasNoTrans, 1.0, $q, $r, 0.0, $a);
 254
 255   my ($aij, $mij);
 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));
 261     }
 262   }
 263
 264 }
 265
 266 sub GSL_LINALG_CHOLESKY_DECOMP : Tests {
 267     my $self = shift;
 268     map{ gsl_matrix_set($self->{matrix}, 0, $_, $_+1)} (0..3);
 269
 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);
 274
 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);
 279
 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);
 284
 285     ok_status(gsl_linalg_cholesky_decomp($self->{matrix}));
 286     my $v = gsl_matrix_diagonal($self->{matrix});
 287     ok_similar(
 288         [ map { gsl_vector_get($v->{vector}, $_)} (0..3) ], [(1)x 4 ]
 289     );
 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);
 296 }
 297 sub GSL_LINALG_HESSENBERG_DECOMP_UNPACK_UNPACK_ACCUM_SET_ZERO : Tests {
 298     my $self = shift;
 299
 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);
 304
 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);
 309
 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);
 314
 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);
 335
 336     my $V = gsl_matrix_alloc(4,4);
 337     ok_status( gsl_linalg_hessenberg_unpack_accum($self->{matrix}, $tau, $V), 0); #I don't know how to test the result of this function...
 338
 339     ok_status( gsl_linalg_hessenberg_set_zero($self->{matrix}) );
 340     for(my $line = 2; $line<4; $line++) {
 341         map { ok_similar(gsl_matrix_get($self->{matrix}, $line, $_), 0, 1e-8, "Set zero") } (0..$line-2);
 342     }
 343 }
 344
 345 sub GSL_LINALG_BIDIAG_DECOMP_UNPACK_UNPACK2_UNPACK_B : Tests {
 346     my $self = shift;
 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);
 351
 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);
 356
 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);
 361
 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);
 368
 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 ] );
 378
 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...
 381 }
 382
 383 Test::Class->runtests;