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