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     my $v = Math::GSL::Vector->new(2);
  86     gsl_matrix_set_identity($m);
  87     ok_status(gsl_eigen_symm($m, $v->raw, $self->{eigen}));
  88     ok_similar( [ sort $v->as_list ], [ 1, 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_similar(gsl_vector_get($eval, 0), 3);
 102     is_similar(gsl_vector_get($eval, 1), 1);
 103     my $x = gsl_matrix_get($evec, 0, 0);
 104
 105     #this is the eigenvector for the eigenvalue 1, which is the second
 106     #eigenvalue in the $eval vector, but the GSL documentation says the first
 107     #eigenvector should correspond to the first eigenvalue... where'e the error?
 108
 109     is_similar(gsl_matrix_get($evec, 0, 1), -$x);
 110
 111     is_similar (sqrt($x**2+$x**2), 1);
 112
 113     $x = gsl_matrix_get($evec, 1, 0);
 114     is_similar(gsl_matrix_get($evec, 1, 1), $x);
 115     is_similar(sqrt($x**2+$x**2), 1);
 116
 117     my $v1 = gsl_vector_alloc(2);
 118     my $v2 = gsl_vector_alloc(2);
 119     gsl_matrix_get_col($v1, $evec, 0);
 120     gsl_matrix_get_col($v2, $evec, 1);
 121     gsl_vector_mul($v1, $v2);
 122     is_similar(gsl_vector_get($v1, 0) + gsl_vector_get($v1, 1) , 0);
 123 }
 124
 125 sub GSL_EIGEN_SYMMV_SORT : Tests {
 126     my $w = gsl_eigen_symmv_alloc(2);
 127     my $m = gsl_matrix_alloc(2,2);
 128     gsl_matrix_set($m, 0, 0, 2);
 129     gsl_matrix_set($m, 0, 1, 1);
 130     gsl_matrix_set($m, 1, 0, 1);
 131     gsl_matrix_set($m, 1, 1, 2);
 132     my $eval = gsl_vector_alloc(2);
 133     my $evec = gsl_matrix_alloc(2,2);
 134     ok_status(gsl_eigen_symmv($m, $eval, $evec, $w));
 135     ok_status(gsl_eigen_symmv_sort ($eval, $evec, $GSL_EIGEN_SORT_VAL_ASC));
 136     is_similar(gsl_vector_get($eval, 0), 1);
 137     is_similar(gsl_vector_get($eval, 1), 3);
 138     my $x = gsl_matrix_get($evec, 0, 0);
 139     ok_similar(gsl_matrix_get($evec, 0, 1), -$x);
 140     ok_similar(sqrt($x**2+$x**2), 1);
 141
 142     $x = gsl_matrix_get($evec, 1, 0);
 143     ok_similar(gsl_matrix_get($evec, 1, 1), $x);
 144     ok_similar(sqrt($x**2+$x**2), 1);
 145
 146     my $v1 = gsl_vector_alloc(2);
 147     my $v2 = gsl_vector_alloc(2);
 148     gsl_matrix_get_col($v1, $evec, 0);
 149     gsl_matrix_get_col($v2, $evec, 1);
 150     gsl_vector_mul($v1, $v2);
 151     is_similar(gsl_vector_get($v1, 0) + gsl_vector_get($v1, 1) , 0);
 152 }
 153
 154 sub GSL_EIGEN_HERM : Tests {
 155     my $matrix  = gsl_matrix_complex_alloc (2, 2);
 156     my $complex = gsl_complex_rect(3,0);
 157     gsl_matrix_complex_set($matrix, 0, 0, $complex);
 158
 159     $complex = gsl_complex_rect(2,1);
 160     gsl_matrix_complex_set($matrix, 0, 1, $complex);
 161
 162     $complex = gsl_complex_rect(2,-1);
 163     gsl_matrix_complex_set($matrix, 1, 0, $complex);
 164
 165     $complex = gsl_complex_rect(1,0);
 166     gsl_matrix_complex_set($matrix, 1, 1, $complex);
 167
 168     my $eigen  = gsl_eigen_herm_alloc(2);
 169     my $vector = gsl_vector_alloc(2);
 170     ok_status(gsl_eigen_herm($matrix, $vector, $eigen));
 171     is_similar(gsl_vector_get($vector, 0), 2+sqrt(6));
 172     is_similar(gsl_vector_get($vector, 1), 2-sqrt(6));
 173 }
 174 sub GSL_EIGEN_HERMV : Tests {
 175     my $matrix  = gsl_matrix_complex_alloc (2, 2);
 176     my $complex = gsl_complex_rect(3,0);
 177     gsl_matrix_complex_set($matrix, 0, 0, $complex);
 178
 179     $complex = gsl_complex_rect(2,1);
 180     gsl_matrix_complex_set($matrix, 0, 1, $complex);
 181
 182     $complex = gsl_complex_rect(2,-1);
 183     gsl_matrix_complex_set($matrix, 1, 0, $complex);
 184
 185     $complex = gsl_complex_rect(1,0);
 186     gsl_matrix_complex_set($matrix, 1, 1, $complex);
 187
 188     my $eigen  = gsl_eigen_hermv_alloc(2);
 189     my $vector = gsl_vector_alloc(2);
 190     my $evec = gsl_matrix_complex_alloc(2,2);
 191     ok_status(gsl_eigen_hermv($matrix, $vector, $evec, $eigen));
 192     is_similar(gsl_vector_get($vector, 0), 2+sqrt(6));
 193     is_similar(gsl_vector_get($vector, 1), 2-sqrt(6));
 194
 195 #    my $x = gsl_matrix_complex_get($evec, 1, 1);
 196 #    ok_similar( [ gsl_parts($x) ], [-0.750532688285823, 0 ], "last row");
 197 #    ok_similar( [ gsl_parts($x) ], [ 2/(-1+sqrt(6)), 1/(-1+sqrt(6)) ], "first value");
 198
 199 }
 200
 201 sub GSL_EIGEN_NONSYMM : Tests {
 202     my $matrix  = gsl_matrix_alloc (2, 2);
 203     gsl_matrix_set($matrix, 0, 0, -12);
 204     gsl_matrix_set($matrix, 1, 0, 7);
 205     gsl_matrix_set($matrix, 0, 1, 65);
 206     gsl_matrix_set($matrix, 1, 1, 59);
 207
 208     my $eigen  = gsl_eigen_nonsymm_alloc(2);
 209     my $vector = gsl_vector_complex_alloc(2);
 210     ok_status(gsl_eigen_nonsymm($matrix, $vector, $eigen));
 211
 212     my $x = gsl_vector_complex_real($vector);
 213     my $y = gsl_vector_complex_imag($vector);
 214
 215     # this interface seems hokey
 216     is_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
 217     is_similar( gsl_vector_get($y->{vector}, 1), 0 );
 218 }
 219
 220 sub GSL_EIGEN_NONSYMM_Z : Tests {
 221     my $matrix  = gsl_matrix_alloc (2, 2);
 222     gsl_matrix_set($matrix, 0, 0, -12);
 223     gsl_matrix_set($matrix, 1, 0, 7);
 224     gsl_matrix_set($matrix, 0, 1, 65);
 225     gsl_matrix_set($matrix, 1, 1, 59);
 226
 227     my $eigen  = gsl_eigen_nonsymm_alloc(2);
 228     my $vector = gsl_vector_complex_alloc(2);
 229     my $Z = gsl_matrix_alloc(2,2);
 230     ok_status(gsl_eigen_nonsymm($matrix, $vector, $eigen));
 231     ok_status(gsl_eigen_nonsymm_Z($matrix,$vector, $Z, $eigen));
 232     ok_similar(gsl_matrix_get($Z, 0, 0), 0.9958842418254068860784291, "Z matrix", 0.005);
 233     ok_similar(gsl_matrix_get($Z, 0, 1), 0.09063430301952179629793610, "Z matrix", 0.1);
 234     ok_similar(gsl_matrix_get($Z, 1, 1), 0.9958842418254068860784291, "Z matrix", 0.005);
 235     ok_similar(gsl_matrix_get($Z, 1, 0), 0.09063430301952179629793610, "Z matrix", 0.1);
 236
 237     my $x = gsl_vector_complex_real($vector);
 238     my $y = gsl_vector_complex_imag($vector);
 239
 240     # this interface seems hokey
 241     is_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
 242     is_similar( gsl_vector_get($y->{vector}, 1), 0 );
 243 }
 244
 245 sub GSL_EIGEN_NONSYMMV : Tests {
 246     my $matrix  = gsl_matrix_alloc (2, 2);
 247     gsl_matrix_set($matrix, 0, 0, -12);
 248     gsl_matrix_set($matrix, 1, 0, 7);
 249     gsl_matrix_set($matrix, 0, 1, 65);
 250     gsl_matrix_set($matrix, 1, 1, 59);
 251
 252     my $evec = gsl_matrix_complex_alloc(2,2);
 253     my $eigen  = gsl_eigen_nonsymmv_alloc(2);
 254     my $vector = gsl_vector_complex_alloc(2);
 255     ok_status(gsl_eigen_nonsymmv($matrix,$vector, $evec, $eigen));
 256
 257
 258     my $x = gsl_vector_complex_real($vector);
 259     my $y = gsl_vector_complex_imag($vector);
 260
 261     # this interface seems hokey
 262     is_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
 263     is_similar( gsl_vector_get($y->{vector}, 1), 0 );
 264
 265     $x = gsl_matrix_complex_get($evec, 1, 0);
 266     ok_similar(gsl_imag($x), 0, "evec matrix");
 267     ok_similar(gsl_real($x), 7/((71/2)+(.5*sqrt(6861))), "evec matrix", 0.01);
 268
 269     $x = gsl_matrix_complex_get($evec, 0, 0);
 270     ok_similar(gsl_imag($x), 0, "evec matrix");
 271     ok_similar(gsl_real($x), 7/((71/2)-(.5*sqrt(6861))), "evec matrix", 0.19);
 272
 273     $x = gsl_matrix_complex_get($evec, 0, 1);
 274     $y = gsl_matrix_complex_get($evec, 1, 1);
 275
 276     ok_similar(gsl_imag($x), 0, "evec matrix");
 277     ok_similar(gsl_imag($y), 0, "evec matrix");
 278
 279     local $TODO = "matlab differences";
 280     ok_similar(gsl_real($x), 1); # this is the value I get with maple
 281     ok_similar(gsl_real($y), 1); # this is the value I get with maple
 282
 283
 284 }
 285 Test::Class->runtests;