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