lib/Math/GSL/BLAS/Test.pm

   1 package Math::GSL::BLAS::Test;
   2 use base q{Test::Class};
   3 use Test::More;
   4 use Math::GSL::BLAS qw/:all/;
   5 use Math::GSL::Vector qw/:all/;
   6 use Math::GSL::Complex qw/:all/;
   7 use Math::GSL::Matrix qw/:all/;
   8 use Math::GSL::CBLAS qw/:all/;
   9 use Math::GSL qw/:all/;
  10 use Data::Dumper;
  11 use Math::GSL::Errno qw/:all/;
  12 use strict;
  13
  14 sub make_fixture : Test(setup) {
  15 }
  16
  17 sub teardown : Test(teardown) {
  18 }
  19
  20 #sub GSL_BLAS_SDSDOT : Tests {
  21 # my $vec1 = gsl_vector_float_alloc(4);
  22 # my $vec2 = gsl_vector_float_alloc(4);
  23 # map { gsl_vector_float_set($vec1, $_, ($_+1)**2) } (0..3);
  24 # map { gsl_vector_float_set($vec2, $_, $_+1) } (0..3);
  25 # my ($x, $result)= gsl_blas_sdsdot(2, $vec1, $vec2);
  26 # this part fail because the vectors should be initiated with gsl_vector_float_alloc...
  27 # however, the gsl_vector_float_alloc function seems to be deprecated, how should I use BLAS level1 function then?
  28 # there's no test suite yet for the BLAS functions in GSL...
  29 #}
  30
  31 sub GSL_BLAS_DDOT : Tests {
  32   my $vec1 = Math::GSL::Vector->new([1,2,3,4,5]);
  33   my $vec2 = Math::GSL::Vector->new([5,4,3,2,1]);
  34   my ($x, $result) = gsl_blas_ddot($vec1->raw, $vec2->raw);
  35   is($x, 0);
  36   is($result,35);
  37 }
  38
  39 sub GSL_BLAS_ZDOTU : Tests {
  40   my $vec1 = gsl_vector_complex_alloc(2);
  41   my $vec2 = gsl_vector_complex_alloc(2);
  42   my $c = gsl_complex_rect(2,1);
  43   gsl_vector_complex_set($vec1,0,$c);
  44   gsl_vector_complex_set($vec2,0,$c);
  45   $c = gsl_complex_rect(1,1);
  46   gsl_vector_complex_set($vec1,1,$c);
  47   gsl_vector_complex_set($vec2,1,$c);
  48   is(gsl_blas_zdotu($vec1, $vec2, $c),0);
  49   is(gsl_real($c), 3);
  50   is(gsl_imag($c), 6);
  51 }
  52
  53 sub GSL_BLAS_ZDOTC : Tests {
  54   my $vec1 = gsl_vector_complex_alloc(2);
  55   my $vec2 = gsl_vector_complex_alloc(2);
  56   my $c = gsl_complex_rect(2,1);
  57   gsl_vector_complex_set($vec1,0,$c);
  58   gsl_vector_complex_set($vec2,0,$c);
  59   $c = gsl_complex_rect(1,1);
  60   gsl_vector_complex_set($vec1,1,$c);
  61   gsl_vector_complex_set($vec2,1,$c);
  62   is(gsl_blas_zdotc($vec1, $vec2, $c),0);
  63   is(gsl_real($c), 7);
  64   is(gsl_imag($c), 0);
  65 }
  66
  67 sub GSL_BLAS_DNRM2 : Tests {
  68   my $vec = Math::GSL::Vector->new([3,4]);
  69   is(gsl_blas_dnrm2($vec->raw), 5);
  70 }
  71
  72
  73 sub GSL_BLAS_DZNRM2 : Tests {
  74   my $vec = gsl_vector_complex_alloc(2);
  75   my $c = gsl_complex_rect(2,1);
  76   gsl_vector_complex_set($vec,0,$c);
  77   $c = gsl_complex_rect(1,1);
  78   gsl_vector_complex_set($vec,1,$c);
  79   is(gsl_blas_dznrm2($vec), sqrt(7));
  80 }
  81
  82 sub GSL_BLAS_DASUM : Tests {
  83   my $vec = Math::GSL::Vector->new([2,-3,4]);
  84   is(gsl_blas_dasum($vec->raw), 9);
  85 }
  86
  87 sub GSL_BLAS_DZASUM : Tests {
  88   my $vec = gsl_vector_complex_alloc(2);
  89   my $c = gsl_complex_rect(2,1);
  90   gsl_vector_complex_set($vec,0,$c);
  91   $c = gsl_complex_rect(1,1);
  92   gsl_vector_complex_set($vec,1,$c);
  93   is(gsl_blas_dzasum($vec), 5);
  94 }
  95
  96 sub GSL_BLAS_DSWAP : Tests {
  97   my $vec1 = Math::GSL::Vector->new([0,1,2]);
  98   my $vec2 = Math::GSL::Vector->new([2,1,0]);
  99   gsl_blas_dswap($vec1->raw, $vec2->raw);
 100   my @got = $vec2->as_list;
 101   map { is($got[$_], $_) } (0..2);
 102   @got = $vec1->as_list;
 103   is($got[0], 2);
 104   is($got[1], 1);
 105   is($got[2], 0);
 106 }
 107
 108 sub GSL_BLAS_ZSWAP : Tests {
 109   local $TODO = "Problem with the output of gsl_vector_complex_get";
 110   my $vec1 = gsl_vector_complex_alloc(2);
 111   my $vec2 = gsl_vector_complex_alloc(2);
 112   my $c = gsl_complex_rect(5,4);
 113   gsl_vector_complex_set($vec1,0,$c);
 114   $c = gsl_complex_rect(2,2);
 115   gsl_vector_complex_set($vec1,1, $c);
 116   $c = gsl_complex_rect(3,3);
 117   print Dumper [ $c ];
 118   gsl_vector_complex_set($vec2,0, $c);
 119   $c = gsl_complex_rect(1,1);
 120   gsl_vector_complex_set($vec2,1, $c);
 121
 122   is(gsl_blas_zswap($vec1, $vec2), 0);
 123   $c = gsl_vector_complex_get($vec1,0);
 124   print Dumper [ $c ];
 125 #  is( gsl_real($c), 3);
 126 #  is( gsl_imag($x), 3);
 127 }
 128
 129 sub GSL_BLAS_DCOPY : Tests {
 130  my $vec1 = Math::GSL::Vector->new([0,1,2]);
 131  my $vec2 = Math::GSL::Vector->new(3);
 132  is(gsl_blas_dcopy($vec1->raw, $vec2->raw),0);
 133  my @got = $vec2->as_list;
 134  map { is($got[$_], $_) } (0..2);
 135 }
 136
 137 sub GSL_BLAS_DAXPY : Tests {
 138  my $vec1 = Math::GSL::Vector->new([0,1,2]);
 139  my $vec2 = Math::GSL::Vector->new([2,3,4]);
 140  is(gsl_blas_daxpy(2,$vec1->raw, $vec2->raw),0);
 141  my @got = $vec2->as_list;
 142  is($got[0], 2);
 143  is($got[1], 5);
 144  is($got[2], 8);
 145 }
 146
 147 sub GSL_BLAS_DSCAL : Tests {
 148  my $vec = Math::GSL::Vector->new([0,1,2]);
 149  gsl_blas_dscal(4, $vec->raw);
 150  my @got = $vec->as_list;
 151  map { is($got[$_], $_*4) } (0..2);
 152 }
 153
 154 sub GSL_BLAS_DROT : Tests {
 155  my $x = Math::GSL::Vector->new([1,2,3]);
 156  my $y = Math::GSL::Vector->new([0,1,2]);
 157  is(gsl_blas_drot($x->raw,$y->raw,2,3),0);
 158  ok_similar( [$x->as_list], [ 2,7,12], 'first vector');
 159  ok_similar( [$y->as_list], [-3,-4,-5], 'second vector');
 160 }
 161
 162 sub GSL_BLAS_DGER : Tests {
 163  my $x = Math::GSL::Vector->new([1,2,3]);
 164  my $y = Math::GSL::Vector->new([0,1,2]);
 165  my $A = Math::GSL::Matrix->new(3,3);
 166  gsl_matrix_set_zero($A->raw);
 167  is(gsl_blas_dger(2, $x->raw, $y->raw, $A->raw),0);
 168  my @got = $A->as_list_row(0);
 169  map { is($got[$_], 0) } (0..2);
 170  @got = $A->as_list_row(1);
 171  map { is($got[$_], ($_+1)*2) } (0..2);
 172  @got = $A->as_list_row(2);
 173  map { is($got[$_], ($_+1)*4) } (0..2);
 174 }
 175
 176 sub GSL_BLAS_ZGERU : Tests {
 177  my $x = gsl_vector_complex_alloc(2);
 178  my $y = gsl_vector_complex_alloc(2);
 179  my $A = gsl_matrix_complex_alloc(2,2);
 180  my $alpha = gsl_complex_rect(2,2);
 181  gsl_vector_complex_set($x, 0, $alpha);
 182  $alpha = gsl_complex_rect(1,2);
 183  gsl_vector_complex_set($x, 1, $alpha);
 184  gsl_vector_complex_set($y, 0, $alpha);
 185  $alpha = gsl_complex_rect(3,2);
 186  gsl_vector_complex_set($y, 1, $alpha);
 187  $alpha = gsl_complex_rect(0,0);
 188  for (my $line=0; $line<2; $line++) {
 189  map { gsl_matrix_complex_set($A, $line, $_, $alpha) } (0..1); }
 190  $alpha = gsl_complex_rect(1,0);
 191  is(gsl_blas_zgeru($alpha, $x, $y, $A),0);
 192
 193  $alpha= gsl_matrix_complex_get($A, 0,0);
 194  ok_similar([gsl_parts($alpha)], [-2, 6]);
 195  $alpha= gsl_matrix_complex_get($A, 1,0);
 196  ok_similar([gsl_parts($alpha)], [-3, 4]);
 197  $alpha= gsl_matrix_complex_get($A, 1,0);
 198  ok_similar([gsl_parts($alpha)], [-3, 4]);
 199  $alpha= gsl_matrix_complex_get($A, 0,1);
 200  ok_similar([gsl_parts($alpha)], [2, 10]);
 201  $alpha= gsl_matrix_complex_get($A, 1,1);
 202  ok_similar([gsl_parts($alpha)], [-1, 8]);
 203 }
 204 1;