1 package Math::GSL::Eigen::Test;
2 use base q{Test::Class};
3 use Test::More tests => 59;
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/;
16 sub make_fixture : Test(setup) {
18 $self->{eigen} = gsl_eigen_symm_alloc(2);
21 sub teardown : Test(teardown) {
24 sub GSL_EIGEN_SYMM_ALLOC : Tests {
26 isa_ok($self->{eigen}, 'Math::GSL::Eigen');
29 sub GSL_EIGEN_SYMMV_ALLOC : Tests {
30 my $eigen = gsl_eigen_symmv_alloc(5);
31 isa_ok($eigen, 'Math::GSL::Eigen');
34 sub GSL_EIGEN_HERM_ALLOC : Tests {
35 my $eigen = gsl_eigen_herm_alloc(5);
36 isa_ok($eigen, 'Math::GSL::Eigen');
39 sub GSL_EIGEN_HERMV_ALLOC : Tests {
40 my $eigen = gsl_eigen_hermv_alloc(5);
41 isa_ok($eigen, 'Math::GSL::Eigen');
44 sub GSL_EIGEN_NONSYMM_ALLOC : Tests {
45 my $eigen = gsl_eigen_nonsymm_alloc(5);
46 isa_ok($eigen, 'Math::GSL::Eigen');
49 sub GSL_EIGEN_NONSYMMV_ALLOC : Tests {
50 my $eigen = gsl_eigen_nonsymmv_alloc(5);
51 isa_ok($eigen, 'Math::GSL::Eigen');
54 sub GSL_EIGEN_GENSYMM_ALLOC : Tests {
55 my $eigen = gsl_eigen_gensymm_alloc(5);
56 isa_ok($eigen, 'Math::GSL::Eigen');
59 sub GSL_EIGEN_GENSYMMV_ALLOC : Tests {
60 my $eigen = gsl_eigen_gensymmv_alloc(5);
61 isa_ok($eigen, 'Math::GSL::Eigen');
64 sub GSL_EIGEN_GENHERM_ALLOC : Tests {
65 my $eigen = gsl_eigen_genherm_alloc(5);
66 isa_ok($eigen, 'Math::GSL::Eigen');
69 sub GSL_EIGEN_GENHERMV_ALLOC : Tests {
70 my $eigen = gsl_eigen_genhermv_alloc(5);
71 isa_ok($eigen, 'Math::GSL::Eigen');
74 sub GSL_EIGEN_GEN_ALLOC : Tests {
75 my $eigen = gsl_eigen_gen_alloc(5);
76 isa_ok($eigen, 'Math::GSL::Eigen');
79 sub GSL_EIGEN_GENV_ALLOC : Tests {
80 my $eigen = gsl_eigen_genv_alloc(5);
81 isa_ok($eigen, 'Math::GSL::Eigen');
84 sub GSL_EIGEN_SYMM : Tests {
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 ] );
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 ok_similar(gsl_vector_get($eval, 0), 3);
104 ok_similar(gsl_vector_get($eval, 1), 1);
105 my $x = gsl_matrix_get($evec, 0, 0);
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?
111 ok_similar(gsl_matrix_get($evec, 0, 1), -$x);
113 ok_similar (sqrt($x**2+$x**2), 1);
115 $x = gsl_matrix_get($evec, 1, 0);
116 ok_similar(gsl_matrix_get($evec, 1, 1), $x);
117 ok_similar(sqrt($x**2+$x**2), 1);
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 ok_similar(gsl_vector_get($v1, 0) + gsl_vector_get($v1, 1) , 0);
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 ok_similar(gsl_vector_get($eval, 0), 1);
139 ok_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);
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);
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 ok_similar(gsl_vector_get($v1, 0) + gsl_vector_get($v1, 1) , 0);
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);
161 $complex = gsl_complex_rect(2,1);
162 gsl_matrix_complex_set($matrix, 0, 1, $complex);
164 $complex = gsl_complex_rect(2,-1);
165 gsl_matrix_complex_set($matrix, 1, 0, $complex);
167 $complex = gsl_complex_rect(1,0);
168 gsl_matrix_complex_set($matrix, 1, 1, $complex);
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 ok_similar(gsl_vector_get($vector, 0), 2+sqrt(6));
174 ok_similar(gsl_vector_get($vector, 1), 2-sqrt(6));
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);
181 $complex = gsl_complex_rect(2,1);
182 gsl_matrix_complex_set($matrix, 0, 1, $complex);
184 $complex = gsl_complex_rect(2,-1);
185 gsl_matrix_complex_set($matrix, 1, 0, $complex);
187 $complex = gsl_complex_rect(1,0);
188 gsl_matrix_complex_set($matrix, 1, 1, $complex);
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 ok_similar(gsl_vector_get($vector, 0), 2+sqrt(6));
195 ok_similar(gsl_vector_get($vector, 1), 2-sqrt(6));
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");
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);
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));
214 my $x = gsl_vector_complex_real($vector);
215 my $y = gsl_vector_complex_imag($vector);
217 # this interface seems hokey
218 ok_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
219 ok_similar( gsl_vector_get($y->{vector}, 1), 0 );
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);
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);
239 my $x = gsl_vector_complex_real($vector);
240 my $y = gsl_vector_complex_imag($vector);
242 # this interface seems hokey
243 ok_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
244 ok_similar( gsl_vector_get($y->{vector}, 1), 0 );
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);
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));
260 my $x = gsl_vector_complex_real($vector);
261 my $y = gsl_vector_complex_imag($vector);
263 # this interface seems hokey
264 ok_similar( gsl_vector_get($x->{vector}, 1), (47/2)+(0.5*sqrt(6861)) );
265 ok_similar( gsl_vector_get($y->{vector}, 1), 0 );
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);
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);
275 $x = gsl_matrix_complex_get($evec, 0, 1);
276 $y = gsl_matrix_complex_get($evec, 1, 1);
278 ok_similar(gsl_imag($x), 0, "evec matrix");
279 ok_similar(gsl_imag($y), 0, "evec matrix");
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
287 Test::Class->runtests;