t/Linalg.t

   1 package Math::GSL::Linalg::Test;
   2 use base q{Test::Class};
   3 use Test::More tests => 55;
   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/;
  16 use Math::Complex;
  17 use Data::Dumper;
  18 use strict;
  19
  20 sub make_fixture : Test(setup) {
  21     my $self = shift;
  22     $self->{matrix} = gsl_matrix_alloc(4, 4);
  23 }
  24
  25 sub teardown : Test(teardown) {
  26     unlink 'linalg' if -f 'linalg';
  27 }
  28
  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);
  39
  40     my ($result, $signum) = gsl_linalg_LU_decomp($base->raw, $permutation);
  41     is_deeply( [ $result, $signum ], [ 0, 1] );
  42
  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);
  47     my $line;
  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);
  54     map {
  55     my @got = $first->row($permutations[$_])->as_list;
  56     my @results = $R->row($_)->as_list;
  57     ok_similar([@got], [@results], "resulting row $_"); } (0..3);
  58 }
  59
  60 sub GSL_LINALG_LU_SOLVE : Tests {
  61     my $self = shift;
  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);
  66
  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);
  71
  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);
  76
  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);
  81
  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);
  87
  88     my $x = gsl_vector_alloc(4);
  89
  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);
  95     ok_similar(
  96         [ map { gsl_vector_get($x, $_) } (1..3) ],
  97         [ 3-10*$value, -2*$value+2, 13*$value-3 ]
  98     );
  99 }
 100
 101 sub GSL_LINALG_LU_SVX : Tests {
 102     my $self = shift;
 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);
 107
 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);
 112
 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);
 117
 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);
 122
 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);
 128
 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);
 134     ok_similar(
 135         [ map { gsl_vector_get($x, $_) } (1..3) ],
 136         [ 3-10*$value, -2*$value+2, 13*$value-3 ]
 137     );
 138 }
 139
 140 sub GSL_LINALG_LU_INVERT : Tests {
 141     my $self = shift;
 142     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 143
 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);
 148
 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);
 153
 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);
 158
 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);
 164
 165     is_similar(gsl_matrix_get($inverse, 0, 0), -9/40);
 166     is_similar(gsl_matrix_get($inverse, 0, 1), 1/40);
 167     is_similar(gsl_matrix_get($inverse, 0, 2), 1/40);
 168     is_similar(gsl_matrix_get($inverse, 0, 3), 11/40);
 169
 170     is_similar(gsl_matrix_get($inverse, 1, 0), 1/40);
 171     is_similar(gsl_matrix_get($inverse, 1, 1), 1/40);
 172     is_similar(gsl_matrix_get($inverse, 1, 2), 11/40);
 173     is_similar(gsl_matrix_get($inverse, 1, 3), -9/40);
 174
 175     is_similar(gsl_matrix_get($inverse, 2, 0), 1/40);
 176     is_similar(gsl_matrix_get($inverse, 2, 1), 11/40);
 177     is_similar(gsl_matrix_get($inverse, 2, 2), -9/40);
 178     is_similar(gsl_matrix_get($inverse, 2, 3), 1/40);
 179
 180     is_similar(gsl_matrix_get($inverse, 3, 0), 11/40);
 181     is_similar(gsl_matrix_get($inverse, 3, 1), -9/40);
 182     is_similar(gsl_matrix_get($inverse, 3, 2), 1/40);
 183     is_similar(gsl_matrix_get($inverse, 3, 3), 1/40);
 184 }
 185
 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);
 193     my $copy = $m->copy;
 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');
 200 }
 201
 202 sub GSL_LINALG_LU_DET : Tests {
 203     my $self = shift;
 204     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 205
 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);
 210
 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);
 215
 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);
 220
 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);
 225 }
 226
 227 sub GSL_LINALG_LU_LNDET : Tests {
 228     my $self = shift;
 229     map { gsl_matrix_set($self->{matrix}, 0, $_, $_+1) } (0..3);
 230
 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);
 235
 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);
 240
 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);
 245
 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));
 250 }
 251
 252 sub GSL_LINALG_QR_DECOMP : Tests {
 253     my $matrix = gsl_matrix_alloc(3,5);
 254     my ($i, $j);
 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));
 258      }
 259     }
 260
 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);
 267
 268     ok_status(gsl_linalg_QR_decomp($matrix, $tau));
 269     is(gsl_linalg_QR_unpack($matrix, $tau, $q, $r), 0);
 270   # compute a = q r
 271   gsl_blas_dgemm ($CblasNoTrans, $CblasNoTrans, 1.0, $q, $r, 0.0, $a);
 272
 273   my ($aij, $mij);
 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));
 279     }
 280   }
 281
 282 }
 283
 284 sub GSL_LINALG_CHOLESKY_DECOMP : Tests {
 285     my $self = shift;
 286     map{ gsl_matrix_set($self->{matrix}, 0, $_, $_+1)} (0..3);
 287
 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);
 292
 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);
 297
 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);
 302
 303     ok_status(gsl_linalg_cholesky_decomp($self->{matrix}));
 304     my $v = gsl_matrix_diagonal($self->{matrix});
 305     ok_similar(
 306         [ map { gsl_vector_get($v->{vector}, $_)} (0..3) ], [(1)x 4 ]
 307     );
 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);
 314 }
 315 sub GSL_LINALG_HESSENBERG_DECOMP_UNPACK_UNPACK_ACCUM_SET_ZERO : Tests {
 316     my $self = shift;
 317
 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);
 322
 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);
 327
 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);
 332
 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     is(gsl_matrix_get($U, 0, 0), 1);
 342     map { is(gsl_matrix_get($U, $_, 0), 0) } (1..3);
 343     map { is(gsl_matrix_get($U, 0, $_), 0) } (1..3);
 344     is_similar(gsl_matrix_get($U, 1, 1), -0.620173672946042309);
 345     is_similar(gsl_matrix_get($U, 1, 2), -0.268847804615518438);
 346     is_similar(gsl_matrix_get($U, 1, 3), 0.736956900597335762);
 347     is_similar(gsl_matrix_get($U, 2, 1), 0.248069469178416908);
 348     is_similar(gsl_matrix_get($U, 2, 2), -0.958442423454322956);
 349     is_similar(gsl_matrix_get($U, 2, 3), -0.140888819231843737);
 350     is_similar(gsl_matrix_get($U, 3, 1), 0.744208407535250749);
 351     is_similar(gsl_matrix_get($U, 3, 2), 0.0954409706385089263);
 352     is_similar(gsl_matrix_get($U, 3, 3), 0.661093690241727594);
 353
 354     my $V = gsl_matrix_alloc(4,4);
 355     ok_status( gsl_linalg_hessenberg_unpack_accum($self->{matrix}, $tau, $V), 0); #I don't know how to test the result of this function...
 356
 357     ok_status( gsl_linalg_hessenberg_set_zero($self->{matrix}) );
 358     for(my $line = 2; $line<4; $line++) {
 359         map { ok_similar(gsl_matrix_get($self->{matrix}, $line, $_), 0, "Set zero", 1e-8) } (0..$line-2);
 360     }
 361 }
 362
 363 sub GSL_LINALG_BIDIAG_DECOMP_UNPACK_UNPACK2_UNPACK_B : Tests {
 364     my $self = shift;
 365     gsl_matrix_set($self->{matrix}, 1, 0, 3);
 366     gsl_matrix_set($self->{matrix}, 1, 1, 2);
 367     gsl_matrix_set($self->{matrix}, 1, 2, 5);
 368     gsl_matrix_set($self->{matrix}, 1, 3, -4);
 369
 370     gsl_matrix_set($self->{matrix}, 1, 0, 5);
 371     gsl_matrix_set($self->{matrix}, 1, 1, -3);
 372     gsl_matrix_set($self->{matrix}, 1, 2, 6);
 373     gsl_matrix_set($self->{matrix}, 1, 3, 9);
 374
 375     gsl_matrix_set($self->{matrix}, 2, 0, -2);
 376     gsl_matrix_set($self->{matrix}, 2, 1, 1);
 377     gsl_matrix_set($self->{matrix}, 2, 2, 5);
 378     gsl_matrix_set($self->{matrix}, 2, 3, 8);
 379
 380     gsl_matrix_set($self->{matrix}, 3, 0, -6);
 381     gsl_matrix_set($self->{matrix}, 3, 1, 7);
 382     gsl_matrix_set($self->{matrix}, 3, 2, 2);
 383     gsl_matrix_set($self->{matrix}, 3, 3, -8);
 384     my $tau_U = gsl_vector_alloc(4);
 385     my $tau_V = gsl_vector_alloc(3);
 386
 387     ok_status(gsl_linalg_bidiag_decomp($self->{matrix}, $tau_U, $tau_V));
 388     my $U = gsl_matrix_alloc(4,4);
 389     my $V = gsl_matrix_alloc(4,4);
 390     my $diag = gsl_vector_alloc(4);
 391     my $superdiag = gsl_vector_alloc(3);
 392     ok_status(gsl_linalg_bidiag_unpack($self->{matrix}, $tau_U, $U, $tau_V, $V, $diag, $superdiag));
 393     is(gsl_matrix_get($V, 0, 0), 1);
 394     ok_similar( [ map { gsl_matrix_get($V, $_, 0) } (1..3) ], [ (0) x 3 ] );
 395     ok_similar( [ map { gsl_matrix_get($V, 0, $_) } (1..3) ], [ (0) x 3 ] );
 396
 397     is_similar(gsl_matrix_get($U, 1, 1), -0.609437002705849772);
 398     is_similar(gsl_matrix_get($U, 1, 2), -0.758604748961341558); #doesn't fit the data I've got...
 399 }
 400
 401 Test::Class->runtests;