t/Eigen.t

   1 package Math::GSL::Eigen::Test;
   2 use base q{Test::Class};
   3 use Test::More;
   4 use Math::GSL          qw/:all/;
   5 use Math::GSL::Test    qw/:all/;
   6 use Math::GSL::Eigen   qw/:all/;
   7 use Math::GSL::Matrix  qw/:all/;
   8 use Math::GSL::Vector  qw/:all/;
   9 use Math::GSL::Complex qw/:all/;
  10 use Math::GSL::Machine qw/:all/;
  11 use Data::Dumper;
  12 use strict;
  13
  14 sub make_fixture : Test(setup) {
  15     my $self = shift;
  16     $self->{eigen} = gsl_eigen_symm_alloc(2);
  17 }
  18
  19 sub teardown : Test(teardown) {
  20 }
  21
  22 sub GSL_EIGEN_SYMM_ALLOC : Tests {
  23     my $self = shift;
  24     isa_ok($self->{eigen}, 'Math::GSL::Eigen');
  25 }
  26
  27 sub GSL_EIGEN_SYMMV_ALLOC : Tests {
  28     my $eigen = gsl_eigen_symmv_alloc(5);
  29     isa_ok($eigen, 'Math::GSL::Eigen');
  30 }
  31
  32 sub GSL_EIGEN_HERM_ALLOC : Tests {
  33     my $eigen = gsl_eigen_herm_alloc(5);
  34     isa_ok($eigen, 'Math::GSL::Eigen');
  35 }
  36
  37 sub GSL_EIGEN_HERMV_ALLOC : Tests {
  38     my $eigen = gsl_eigen_hermv_alloc(5);
  39     isa_ok($eigen, 'Math::GSL::Eigen');
  40 }
  41
  42 sub GSL_EIGEN_NONSYMM_ALLOC : Tests {
  43     my $eigen = gsl_eigen_nonsymm_alloc(5);
  44     isa_ok($eigen, 'Math::GSL::Eigen');
  45 }
  46
  47 sub GSL_EIGEN_NONSYMMV_ALLOC : Tests {
  48     my $eigen = gsl_eigen_nonsymmv_alloc(5);
  49     isa_ok($eigen, 'Math::GSL::Eigen');
  50 }
  51
  52 sub GSL_EIGEN_GENSYMM_ALLOC : Tests {
  53     my $eigen = gsl_eigen_gensymm_alloc(5);
  54     isa_ok($eigen, 'Math::GSL::Eigen');
  55 }
  56
  57 sub GSL_EIGEN_GENSYMMV_ALLOC : Tests {
  58     my $eigen = gsl_eigen_gensymmv_alloc(5);
  59     isa_ok($eigen, 'Math::GSL::Eigen');
  60 }
  61
  62 sub GSL_EIGEN_GENHERM_ALLOC : Tests {
  63     my $eigen = gsl_eigen_genherm_alloc(5);
  64     isa_ok($eigen, 'Math::GSL::Eigen');
  65 }
  66
  67 sub GSL_EIGEN_GENHERMV_ALLOC : Tests {
  68     my $eigen = gsl_eigen_genhermv_alloc(5);
  69     isa_ok($eigen, 'Math::GSL::Eigen');
  70 }
  71
  72 sub GSL_EIGEN_GEN_ALLOC : Tests {
  73     my $eigen = gsl_eigen_gen_alloc(5);
  74     isa_ok($eigen, 'Math::GSL::Eigen');
  75 }
  76
  77 sub GSL_EIGEN_GENV_ALLOC : Tests {
  78     my $eigen = gsl_eigen_genv_alloc(5);
  79     isa_ok($eigen, 'Math::GSL::Eigen');
  80 }
  81
  82 sub GSL_EIGEN_SYMM : Tests {
  83     my $self = shift;
  84     my $m = gsl_matrix_alloc(2,2);
  85     gsl_matrix_set_identity($m);
  86     my $v = gsl_vector_alloc(2);
  87     ok_status(gsl_eigen_symm($m, $v, $self->{eigen}));
  88     map { is(gsl_vector_get($v, $_), 1) } (0..1);
  89 }
  90
  91 sub GSL_EIGEN_SYMMV : Tests {
  92     my $w = gsl_eigen_symmv_alloc(2);
  93     my $m = gsl_matrix_alloc(2,2);
  94     gsl_matrix_set($m, 0, 0, 2);
  95     gsl_matrix_set($m, 0, 1, 1);
  96     gsl_matrix_set($m, 1, 0, 1);
  97     gsl_matrix_set($m, 1, 1, 2);
  98     my $eval = gsl_vector_alloc(2);
  99     my $evec = gsl_matrix_alloc(2,2);
 100     ok_status(gsl_eigen_symmv($m, $eval, $evec, $w));
 101     is(gsl_vector_get($eval, 0), 3);
 102     is(gsl_vector_get($eval, 1), 1);
 103     my $x = gsl_matrix_get($evec, 0, 0);
 104     is (gsl_matrix_get($evec, 0, 1), -$x); #this is the eigenvector for the eigenvalue 1, which is the second eigenvalue in the $eval vector, but the GSL documentation says the first eigenvector should correspond to the first eigenvalue... where'e the error?
 105     is (sqrt($x**2+$x**2), 1);
 106
 107     $x = gsl_matrix_get($evec, 1, 0);
 108     is (gsl_matrix_get($evec, 1, 1), $x);
 109     is (sqrt($x**2+$x**2), 1);
 110
 111     my $v1 = gsl_vector_alloc(2);
 112     my $v2 = gsl_vector_alloc(2);
 113     gsl_matrix_get_col($v1, $evec, 0);
 114     gsl_matrix_get_col($v2, $evec, 1);
 115     gsl_vector_mul($v1, $v2);
 116     is(gsl_vector_get($v1, 0) + gsl_vector_get($v1, 1) , 0);
 117 }
 118
 119 sub GSL_EIGEN_SYMMV_SORT : Tests {
 120     my $w = gsl_eigen_symmv_alloc(2);
 121     my $m = gsl_matrix_alloc(2,2);
 122     gsl_matrix_set($m, 0, 0, 2);
 123     gsl_matrix_set($m, 0, 1, 1);
 124     gsl_matrix_set($m, 1, 0, 1);
 125     gsl_matrix_set($m, 1, 1, 2);
 126     my $eval = gsl_vector_alloc(2);
 127     my $evec = gsl_matrix_alloc(2,2);
 128     ok_status(gsl_eigen_symmv($m, $eval, $evec, $w));
 129     ok_status(gsl_eigen_symmv_sort ($eval, $evec, $GSL_EIGEN_SORT_VAL_ASC));
 130     is(gsl_vector_get($eval, 0), 1);
 131     is(gsl_vector_get($eval, 1), 3);
 132     my $x = gsl_matrix_get($evec, 0, 0);
 133     ok_similar(gsl_matrix_get($evec, 0, 1), -$x);
 134     ok_similar(sqrt($x**2+$x**2), 1);
 135
 136     $x = gsl_matrix_get($evec, 1, 0);
 137     ok_similar(gsl_matrix_get($evec, 1, 1), $x);
 138     ok_similar(sqrt($x**2+$x**2), 1);
 139
 140     my $v1 = gsl_vector_alloc(2);
 141     my $v2 = gsl_vector_alloc(2);
 142     gsl_matrix_get_col($v1, $evec, 0);
 143     gsl_matrix_get_col($v2, $evec, 1);
 144     gsl_vector_mul($v1, $v2);
 145     is(gsl_vector_get($v1, 0) + gsl_vector_get($v1, 1) , 0);
 146 }
 147
 148 sub GSL_EIGEN_HERM : Tests {
 149     my $matrix  = gsl_matrix_complex_alloc (2, 2);
 150     my $complex = gsl_complex_rect(3,0);
 151     gsl_matrix_complex_set($matrix, 0, 0, $complex);
 152
 153     $complex = gsl_complex_rect(2,1);
 154     gsl_matrix_complex_set($matrix, 0, 1, $complex);
 155
 156     $complex = gsl_complex_rect(2,-1);
 157     gsl_matrix_complex_set($matrix, 1, 0, $complex);
 158
 159     $complex = gsl_complex_rect(1,0);
 160     gsl_matrix_complex_set($matrix, 1, 1, $complex);
 161
 162     my $eigen  = gsl_eigen_herm_alloc(2);
 163     my $vector = gsl_vector_alloc(2);
 164     is(gsl_eigen_herm($matrix, $vector, $eigen), 0);
 165     is (gsl_vector_get($vector, 0), 2+sqrt(6));
 166     is (gsl_vector_get($vector, 1), 2-sqrt(6));
 167 }
 168 sub GSL_EIGEN_HERMV : Tests {
 169     my $matrix  = gsl_matrix_complex_alloc (2, 2);
 170     my $complex = gsl_complex_rect(3,0);
 171     gsl_matrix_complex_set($matrix, 0, 0, $complex);
 172
 173     $complex = gsl_complex_rect(2,1);
 174     gsl_matrix_complex_set($matrix, 0, 1, $complex);
 175
 176     $complex = gsl_complex_rect(2,-1);
 177     gsl_matrix_complex_set($matrix, 1, 0, $complex);
 178
 179     $complex = gsl_complex_rect(1,0);
 180     gsl_matrix_complex_set($matrix, 1, 1, $complex);
 181
 182     my $eigen  = gsl_eigen_hermv_alloc(2);
 183     my $vector = gsl_vector_alloc(2);
 184     my $evec = gsl_matrix_complex_alloc(2,2);
 185     is(gsl_eigen_hermv($matrix, $vector, $evec, $eigen), 0);
 186     is (gsl_vector_get($vector, 0), 2+sqrt(6));
 187     is (gsl_vector_get($vector, 1), 2-sqrt(6));
 188
 189 #    my $x = gsl_matrix_complex_get($evec, 1, 1);
 190 #    ok_similar( [ gsl_parts($x) ], [-0.750532688285823, 0 ], "last row");
 191 #    ok_similar( [ gsl_parts($x) ], [ 2/(-1+sqrt(6)), 1/(-1+sqrt(6)) ], "first value");
 192
 193 }
 194
 195 sub GSL_EIGEN_NONSYMM : Tests {
 196     my $matrix  = gsl_matrix_alloc (2, 2);
 197     gsl_matrix_set($matrix, 0, 0, -12);
 198     gsl_matrix_set($matrix, 1, 0, 7);
 199     gsl_matrix_set($matrix, 0, 1, 65);
 200     gsl_matrix_set($matrix, 1, 1, 59);
 201
 202     my $eigen  = gsl_eigen_nonsymm_alloc(2);
 203     my $vector = gsl_vector_complex_alloc(2);
 204     is(gsl_eigen_nonsymm($matrix, $vector, $eigen), 0);
 205
 206     my $x = gsl_vector_complex_real($vector);
 207     my $y = gsl_vector_complex_imag($vector);
 208
 209     # this interface seems hokey
 210     is_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
 211     is_similar( gsl_vector_get($y->{vector}, 1), 0 );
 212 }
 213
 214 sub GSL_EIGEN_NONSYMM_Z : Tests {
 215     my $matrix  = gsl_matrix_alloc (2, 2);
 216     gsl_matrix_set($matrix, 0, 0, -12);
 217     gsl_matrix_set($matrix, 1, 0, 7);
 218     gsl_matrix_set($matrix, 0, 1, 65);
 219     gsl_matrix_set($matrix, 1, 1, 59);
 220
 221     my $eigen  = gsl_eigen_nonsymm_alloc(2);
 222     my $vector = gsl_vector_complex_alloc(2);
 223     my $Z = gsl_matrix_alloc(2,2);
 224     is(gsl_eigen_nonsymm($matrix, $vector, $eigen), 0);
 225     is(gsl_eigen_nonsymm_Z($matrix,$vector, $Z, $eigen), 0);
 226     ok(is_similar(gsl_matrix_get($Z, 0, 0), 0.9958842418254068860784291, 0.005));
 227     ok_similar(gsl_matrix_get($Z, 0, 1), 0.09063430301952179629793610, "Z matrix", 0.1);
 228     ok_similar(gsl_matrix_get($Z, 1, 1), 0.9958842418254068860784291, "Z matrix", 0.005);
 229     ok_similar(gsl_matrix_get($Z, 1, 0), 0.09063430301952179629793610, "Z matrix", 0.1);
 230
 231     my $x = gsl_vector_complex_real($vector);
 232     my $y = gsl_vector_complex_imag($vector);
 233
 234     # this interface seems hokey
 235     is_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
 236     is_similar( gsl_vector_get($y->{vector}, 1), 0 );
 237 }
 238
 239 sub GSL_EIGEN_NONSYMMV_Z : Tests {
 240     my $matrix  = gsl_matrix_alloc (2, 2);
 241     gsl_matrix_set($matrix, 0, 0, -12);
 242     gsl_matrix_set($matrix, 1, 0, 7);
 243     gsl_matrix_set($matrix, 0, 1, 65);
 244     gsl_matrix_set($matrix, 1, 1, 59);
 245
 246     my $evec = gsl_matrix_complex_alloc(2,2);
 247     my $eigen  = gsl_eigen_nonsymmv_alloc(2);
 248     my $vector = gsl_vector_complex_alloc(2);
 249     my $Z = gsl_matrix_alloc(2,2);
 250     ok_status(gsl_eigen_nonsymmv_Z($matrix,$vector, $evec, $Z, $eigen));
 251
 252
 253     my $x = gsl_vector_complex_real($vector);
 254     my $y = gsl_vector_complex_imag($vector);
 255
 256     # this interface seems hokey
 257     is_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
 258     is_similar( gsl_vector_get($y->{vector}, 1), 0 );
 259
 260     $x = gsl_matrix_complex_get($evec, 1, 0);
 261     ok_similar(gsl_imag($x), 0, "evec matrix");
 262     ok_similar(gsl_real($x), 7/((71/2)+(.5*sqrt(6861))), "evec matrix", 0.01);
 263
 264     $x = gsl_matrix_complex_get($evec, 0, 0);
 265     ok_similar(gsl_imag($x), 0, "evec matrix");
 266     ok_similar(gsl_real($x), 7/((71/2)-(.5*sqrt(6861))), "evec matrix", 0.19);
 267
 268     $x = gsl_matrix_complex_get($evec, 0, 1);
 269     $y = gsl_matrix_complex_get($evec, 1, 1);
 270
 271     ok_similar(gsl_imag($x), 0, "evec matrix");
 272     ok_similar(gsl_imag($y), 0, "evec matrix");
 273
 274     local $TODO = "matlab differences";
 275
 276     ok_similar(gsl_real($x), 1); # this is the value I get with maple
 277     ok_similar(gsl_real($y), 1); # this is the value I get with maple
 278
 279     ok_similar([ gsl_matrix_get($Z, 0, 0)], [0.9958842418254068860784291] );
 280     ok_similar([ gsl_matrix_get($Z, 0, 1)], [0.09063430301952179629793610] );
 281     ok_similar([ gsl_matrix_get($Z, 1, 1)], [0.9958842418254068860784291] );
 282     ok_similar([ gsl_matrix_get($Z, 1, 0)], [0.09063430301952179629793610] );
 283
 284 }
 285 Test::Class->runtests;