t/Eigen.t

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