examples/matrix/nonsymmetric_eigen

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