1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
#!/usr/bin/perl -w
use Math::GSL qw/:all/;
use Math::GSL::Eigen qw/:all/;
use Math::GSL::Matrix qw/:all/;
use Math::GSL::Vector qw/:all/;
use Math::GSL::Complex qw/:all/;
use Math::GSL::Errno qw/:all/;
use Math::Complex;
use Data::Dumper;
use strict;
my $matrix = Math::GSL::Matrix->new(2,2)
->set_row(0, [0,-1] )
->set_row(1, [1, 0] );
print <<STUFF;
Finding eigenvalue/eigenvectors for
[ 0 -1 ]
[ 1 0 ]
STUFF
my $evec = gsl_matrix_complex_alloc(2,2);
my $eigen = gsl_eigen_nonsymmv_alloc(2);
my $vector = gsl_vector_complex_alloc(2);
# this actually computes the eigenvalues and vectors
gsl_eigen_nonsymmv($matrix->raw,$vector, $evec, $eigen);
my $x = gsl_vector_complex_real($vector);
my $y = gsl_vector_complex_imag($vector);
my $eig1 = as_complex( gsl_vector_get($x->{vector}, 0) , gsl_vector_get($y->{vector}, 0) );
my $eig2 = as_complex( gsl_vector_get($x->{vector}, 1) , gsl_vector_get($y->{vector}, 1) );
my $u1 = as_complex( gsl_matrix_complex_get($evec, 0, 0) );
my $u2 = as_complex( gsl_matrix_complex_get($evec, 1, 0) );
my $v1 = as_complex( gsl_matrix_complex_get($evec, 0, 1) );
my $v2 = as_complex( gsl_matrix_complex_get($evec, 1, 1) );
print <<ANSWER;
Eigenvectors:
u = ($u1,$u2)
v = ($v1,$v2)
Eigenvalues:
lambda_0 = $eig1
lambda_1 = $eig2
ANSWER
sub as_complex {
my ($w,$v) = @_;
my ($x,$y);
if( ref $w ) {
($x,$y) = (gsl_real($w),gsl_imag($w));
} else {
($x,$y) = ($w,$v);
}
my $z = Math::Complex->make( $x, $y);
return qq{$z};
}
|
