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   ok_status($x,$GSL_SUCCESS);
  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   ok_similar( [0 .. 2], [ $vec2->as_list ] );
 101   ok_similar( [2, 1,0], [ $vec1->as_list ] );
 102 }
 103
 104 sub GSL_BLAS_ZSWAP : Tests {
 105   local $TODO = "Problem with the output of gsl_vector_complex_get";
 106   my $vec1 = gsl_vector_complex_alloc(2);
 107   my $vec2 = gsl_vector_complex_alloc(2);
 108   my $c = gsl_complex_rect(5,4);
 109   gsl_vector_complex_set($vec1,0,$c);
 110   $c = gsl_complex_rect(2,2);
 111   gsl_vector_complex_set($vec1,1, $c);
 112   $c = gsl_complex_rect(3,3);
 113   print Dumper [ $c ];
 114   gsl_vector_complex_set($vec2,0, $c);
 115   $c = gsl_complex_rect(1,1);
 116   gsl_vector_complex_set($vec2,1, $c);
 117
 118   is(gsl_blas_zswap($vec1, $vec2), 0);
 119   $c = gsl_vector_complex_get($vec1,0);
 120   print Dumper [ $c ];
 121 #  is( gsl_real($c), 3);
 122 #  is( gsl_imag($x), 3);
 123 }
 124
 125 sub GSL_BLAS_DCOPY : Tests {
 126  my $vec1 = Math::GSL::Vector->new([0,1,2]);
 127  my $vec2 = Math::GSL::Vector->new(3);
 128  is(gsl_blas_dcopy($vec1->raw, $vec2->raw),0);
 129  my @got = $vec2->as_list;
 130  map { is($got[$_], $_) } (0..2);
 131 }
 132
 133 sub GSL_BLAS_DAXPY : Tests {
 134  my $vec1 = Math::GSL::Vector->new([0,1,2]);
 135  my $vec2 = Math::GSL::Vector->new([2,3,4]);
 136  is(gsl_blas_daxpy(2,$vec1->raw, $vec2->raw),0);
 137  my @got = $vec2->as_list;
 138  is($got[0], 2);
 139  is($got[1], 5);
 140  is($got[2], 8);
 141 }
 142
 143 sub GSL_BLAS_DSCAL : Tests {
 144  my $vec = Math::GSL::Vector->new([0,1,2]);
 145  gsl_blas_dscal(4, $vec->raw);
 146  my @got = $vec->as_list;
 147  map { is($got[$_], $_*4) } (0..2);
 148 }
 149
 150 sub GSL_BLAS_DROT : Tests {
 151  my $x = Math::GSL::Vector->new([1,2,3]);
 152  my $y = Math::GSL::Vector->new([0,1,2]);
 153  is(gsl_blas_drot($x->raw,$y->raw,2,3),0);
 154  ok_similar( [$x->as_list], [ 2,7,12], 'first vector');
 155  ok_similar( [$y->as_list], [-3,-4,-5], 'second vector');
 156 }
 157
 158 sub GSL_BLAS_DGER : Tests {
 159 local $TODO = "how do you compute a rank-1 update?";
 160  my $x = Math::GSL::Vector->new([1,2,3]);
 161  my $y = Math::GSL::Vector->new([0,1,2]);
 162  my $A = Math::GSL::Matrix->new(3,3);
 163  gsl_matrix_set_zero($A->raw);
 164  is(gsl_blas_dger(2, $x->raw, $y->raw, $A->raw),0);
 165 # ok_similar([$A->row(0)->as_list], [0,2,4]);
 166 # ok_similar([$A->row(1)->as_list], [0,4,8]);
 167 # ok_similar([$A->row(2)->as_list], [0,6,12]);
 168 }
 169
 170 sub GSL_BLAS_ZGERU : Tests {
 171  my $x = gsl_vector_complex_alloc(2);
 172  my $y = gsl_vector_complex_alloc(2);
 173  my $A = gsl_matrix_complex_alloc(2,2);
 174  my $alpha = gsl_complex_rect(2,2);
 175  gsl_vector_complex_set($x, 0, $alpha);
 176  $alpha = gsl_complex_rect(1,2);
 177  gsl_vector_complex_set($x, 1, $alpha);
 178  gsl_vector_complex_set($y, 0, $alpha);
 179  $alpha = gsl_complex_rect(3,2);
 180  gsl_vector_complex_set($y, 1, $alpha);
 181  $alpha = gsl_complex_rect(0,0);
 182  for (my $line=0; $line<2; $line++) {
 183  map { gsl_matrix_complex_set($A, $line, $_, $alpha) } (0..1); }
 184  $alpha = gsl_complex_rect(1,0);
 185  is(gsl_blas_zgeru($alpha, $x, $y, $A),0);
 186
 187  $alpha= gsl_matrix_complex_get($A, 0,0);
 188  ok_similar([gsl_parts($alpha)], [-2, 6]);
 189  $alpha= gsl_matrix_complex_get($A, 1,0);
 190  ok_similar([gsl_parts($alpha)], [-3, 4]);
 191  $alpha= gsl_matrix_complex_get($A, 1,0);
 192  ok_similar([gsl_parts($alpha)], [-3, 4]);
 193  $alpha= gsl_matrix_complex_get($A, 0,1);
 194  ok_similar([gsl_parts($alpha)], [2, 10]);
 195  $alpha= gsl_matrix_complex_get($A, 1,1);
 196  ok_similar([gsl_parts($alpha)], [-1, 8]);
 197 }
 198
 199 sub GSL_BLAS_DGEMV : Tests {
 200  my $x = Math::GSL::Vector->new([1,2,3]);
 201  my $y = Math::GSL::Vector->new([0,1,2]);
 202  my $A = Math::GSL::Matrix->new(3,3);
 203  gsl_matrix_set_identity($A->raw);
 204  is(gsl_blas_dgemv($CblasNoTrans, 2, $A->raw, $x->raw,2, $y->raw),0);
 205  my @got = $y->as_list;
 206  is($got[0], 2);
 207  is($got[1], 6);
 208  is($got[2], 10);
 209 }
 210
 211 sub GSL_BLAS_DTRMV : Tests {
 212  local $TODO = "Problem with the output of gsl_vector_complex_get";
 213  my $x = Math::GSL::Vector->new([1,2,3]);
 214  my $A = Math::GSL::Matrix->new(3,3);
 215  gsl_matrix_set($A->raw, 0,0,3);
 216  gsl_matrix_set($A->raw, 1,1,3);
 217  gsl_matrix_set($A->raw, 2,2,3);
 218  gsl_matrix_set($A->raw, 0,1,2);
 219  gsl_matrix_set($A->raw, 0,2,3);
 220  gsl_matrix_set($A->raw, 1,2,4);
 221  is(gsl_blas_dtrmv($CblasLower, $CblasNoTrans, $CblasNonUnit, $A->raw, $x->raw),0);
 222  my @got = $x->as_list;
 223  is($got[0], 3);
 224  is($got[1], 6);
 225  is($got[2], 9);
 226 }
 227
 228 sub GSL_BLAS_DTRSV : Tests {
 229  local $TODO = "I don't understand what the result is supposed to be...";
 230  my $x = Math::GSL::Vector->new([40,40,40,40]);
 231  my $A = Math::GSL::Matrix->new(4,4);
 232  map { gsl_matrix_set($A->raw, $_,0,$_+1); } (0..3);
 233  map { gsl_matrix_set($A->raw, $_,1,$_+2); } (0..2);
 234  gsl_matrix_set($A->raw, 3,1,1);
 235  map { gsl_matrix_set($A->raw, $_,2,$_+3); } (0..1);
 236  map { gsl_matrix_set($A->raw, $_,2,$_-1); } (2..3);
 237  map { gsl_matrix_set($A->raw, $_,3,$_); } (1..3);
 238  gsl_matrix_set($A->raw, 0,3,4);
 239  is(gsl_blas_dtrsv($CblasLower, $CblasNoTrans, $CblasNonUnit, $A->raw, $x->raw),0);
 240  my @got = $A->row(0)->as_list;
 241 }
 242
 243 sub GSL_BLAS_DROTG : Tests {
 244  my $a = [1];
 245  my $b = [2];
 246  my $c = [0];
 247  my $s = [0];
 248  ok_status(gsl_blas_drotg($a, $b, $c, $s), $GSL_SUCCESS);
 249  local $TODO = "need a typemap for outputting arrays";
 250  ok_similar( [$c->[0], $s->[0]], [ 1/sqrt(5) , 2/sqrt(5)  ] );
 251 }
 252
 253 sub GSL_BLAS_DSYMV : Tests {
 254  my $x = Math::GSL::Vector->new([1,2,3]);
 255  my $y = Math::GSL::Vector->new([3,2,1]);
 256  my $A = Math::GSL::Matrix->new(3,3);
 257  map { gsl_matrix_set($A->raw, $_,0,$_+1); } (0..2);
 258  gsl_matrix_set($A->raw, 0, 1, 2);
 259  gsl_matrix_set($A->raw, 1, 1, 1);
 260  gsl_matrix_set($A->raw, 2, 1, 2);
 261  gsl_matrix_set($A->raw, 0, 2, 3);
 262  gsl_matrix_set($A->raw, 1, 2, 2);
 263  gsl_matrix_set($A->raw, 2, 2, 1);
 264  is(gsl_blas_dsymv($CblasLower, 2, $A->raw, $x->raw, 3, $y->raw),0);
 265  my @got = $y->as_list;
 266  ok_similar( [@got], [37,26,23]);
 267 }
 268
 269 sub GSL_BLAS_DSYR : Tests {
 270  local $TODO = "How do you compute a rank 1 update?";
 271  my $x = Math::GSL::Vector->new([1,2,3]);
 272  my $A = Math::GSL::Matrix->new(3,3);
 273  gsl_matrix_set_zero($A->raw);
 274
 275  is(gsl_blas_dsyr($CblasLower, 2, $x->raw, $A->raw),0);
 276  my @got = $A->row(0)->as_list;
 277 # ok_similar([ @got ], [2,4,6]);
 278  @got = $A->row(1)->as_list;
 279 # ok_similar([ @got ], [0,8,12]);
 280  @got = $A->row(2)->as_list;
 281 # ok_similar([ @got ], [0,0,18]);
 282 }
 283
 284 sub GSL_BLAS_ZHER : Tests {
 285  my $x = gsl_vector_complex_alloc(2);
 286  my $A = gsl_matrix_complex_alloc(2,2);
 287  my $alpha = gsl_complex_rect(0,0);
 288  for(my $line=0; $line<2; $line++){
 289  map { gsl_matrix_complex_set($A, $_, $line, $alpha) } (0..1); }
 290  $alpha = gsl_complex_rect(1,2);
 291  gsl_vector_complex_set($x, 0, $alpha);
 292  $alpha = gsl_complex_rect(2,2);
 293  gsl_vector_complex_set($x, 1, $alpha);
 294  $alpha = gsl_complex_rect(1,1);
 295  is(gsl_blas_zher($CblasLower, 2, $x, $A),0);
 296  my @got = gsl_parts(gsl_matrix_complex_get($A,0,0));
 297  ok_similar([@got], [10, 0]);
 298  @got = gsl_parts(gsl_matrix_complex_get($A,1,0));
 299  ok_similar([@got], [12, -4]);
 300  @got = gsl_parts(gsl_matrix_complex_get($A,1,1));
 301  ok_similar([@got], [16, 0]);
 302  @got = gsl_parts(gsl_matrix_complex_get($A,0,1));
 303  ok_similar([@got], [0, 0]);
 304 }
 305
 306 sub GSL_BLAS_DSYR2 : Tests {
 307  local $TODO = "How do you compute a rank 2 update?";
 308  my $x = Math::GSL::Vector->new([1,2,3]);
 309  my $y = Math::GSL::Vector->new([3,2,1]);
 310  my $A = Math::GSL::Matrix->new(3,3);
 311  gsl_matrix_set_zero($A->raw);
 312  map { gsl_matrix_set($A->raw, $_, 0, ($_+1)**2) } (0..2);
 313  map { gsl_matrix_set($A->raw, $_, 1, ($_+1)**2) } (1..2);
 314  gsl_matrix_set($A->raw, $_, 1, ($_+1)**2);
 315  gsl_matrix_set($A->raw, 0, 1, 4);
 316  gsl_matrix_set($A->raw, 0, 2, 9);
 317  gsl_matrix_set($A->raw, 1, 2, 9);
 318  gsl_matrix_set($A->raw, 2, 2, 3);
 319  is(gsl_blas_dsyr2($CblasLower, 2, $x->raw, $y->raw, $A->raw),0);
 320  my @got = $A->row(0)->as_list;
 321 # ok_similar([@got], [13, 20, 29]);
 322  @got = $A->row(1)->as_list;
 323 # ok_similar([@got], [4, 20, 25]);
 324  @got = $A->row(2)->as_list;
 325 # ok_similar([@got], [9, 9, 15]);
 326 }
 327
 328 sub GSL_BLAS_DGEMM : Tests {
 329  my $A = Math::GSL::Matrix->new(2,2);
 330  gsl_matrix_set($A->raw, 0,0,1);
 331  gsl_matrix_set($A->raw, 1,0,3);
 332  gsl_matrix_set($A->raw, 0,1,4);
 333  gsl_matrix_set($A->raw, 1,1,2);
 334  my $B = Math::GSL::Matrix->new(2,2);
 335  gsl_matrix_set($B->raw, 0,0,2);
 336  gsl_matrix_set($B->raw, 1,0,5);
 337  gsl_matrix_set($B->raw, 0,1,1);
 338  gsl_matrix_set($B->raw, 1,1,3);
 339  my $C = Math::GSL::Matrix->new(2,2);
 340  gsl_matrix_set_zero($C->raw);
 341  is(gsl_blas_dgemm($CblasNoTrans, $CblasNoTrans, 1, $A->raw, $B->raw, 1, $C->raw),0);
 342  my @got = $C->row(0)->as_list;
 343  ok_similar([@got], [22, 13]);
 344  @got = $C->row(1)->as_list;
 345  ok_similar([@got], [16, 9]);
 346 }
 347
 348 sub GSL_BLAS_DSYMM : Tests {
 349  my $A = Math::GSL::Matrix->new(2,2);
 350  gsl_matrix_set($A->raw, 0,0,1);
 351  gsl_matrix_set($A->raw, 1,0,3);
 352  gsl_matrix_set($A->raw, 0,1,4);
 353  gsl_matrix_set($A->raw, 1,1,2);
 354  my $B = Math::GSL::Matrix->new(2,2);
 355  gsl_matrix_set($B->raw, 0,0,2);
 356  gsl_matrix_set($B->raw, 1,0,5);
 357  gsl_matrix_set($B->raw, 0,1,1);
 358  gsl_matrix_set($B->raw, 1,1,3);
 359  my $C = Math::GSL::Matrix->new(2,2);
 360  gsl_matrix_set_zero($C->raw);
 361  is(gsl_blas_dsymm($CblasLeft, $CblasUpper, 1, $A->raw, $B->raw, 1, $C->raw),0);
 362  my @got = $C->row(0)->as_list;
 363  ok_similar([@got], [22, 13]);
 364  @got = $C->row(1)->as_list;
 365  ok_similar([@got], [18, 10]);
 366 }
 367
 368 sub GSL_BLAS_ZGEMM : Tests {
 369  my $A = gsl_matrix_complex_alloc(2,2);
 370  my $alpha = gsl_complex_rect(1,2);
 371  gsl_matrix_complex_set($A, 0,0,$alpha);
 372  $alpha = gsl_complex_rect(3,1);
 373  gsl_matrix_complex_set($A, 0,1,$alpha);
 374  $alpha = gsl_complex_rect(2,2);
 375  gsl_matrix_complex_set($A, 1,0,$alpha);
 376  $alpha = gsl_complex_rect(0,2);
 377  gsl_matrix_complex_set($A, 1,1,$alpha);
 378
 379  my $B = gsl_matrix_complex_alloc(2,2);
 380  $alpha = gsl_complex_rect(1,1);
 381  gsl_matrix_complex_set($B, 0,0,$alpha);
 382  $alpha = gsl_complex_rect(2,2);
 383  gsl_matrix_complex_set($B, 0,1,$alpha);
 384  $alpha = gsl_complex_rect(1,2);
 385  gsl_matrix_complex_set($B, 1,0,$alpha);
 386  $alpha = gsl_complex_rect(1,3);
 387  gsl_matrix_complex_set($B, 1,1,$alpha);
 388
 389  my $C = gsl_matrix_complex_alloc(2,2);
 390  $alpha = gsl_complex_rect(0,0);
 391  map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 392  map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 393
 394  $alpha = gsl_complex_rect(2,0);
 395  my $beta = gsl_complex_rect(1,0);
 396  is(gsl_blas_zgemm($CblasNoTrans, $CblasNoTrans, $alpha, $A, $B, $beta, $C),0);
 397  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [0,20]);
 398  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [-4,32]);
 399  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [-8,12]);
 400  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [-12,20]);
 401 }
 402
 403 sub GSL_BLAS_ZSYMM : Tests {
 404  my $A = gsl_matrix_complex_alloc(2,2);
 405  my $alpha = gsl_complex_rect(1,2);
 406  gsl_matrix_complex_set($A, 0,0,$alpha);
 407  $alpha = gsl_complex_rect(3,1);
 408  gsl_matrix_complex_set($A, 0,1,$alpha);
 409  $alpha = gsl_complex_rect(3,1);
 410  gsl_matrix_complex_set($A, 1,0,$alpha);
 411  $alpha = gsl_complex_rect(0,2);
 412  gsl_matrix_complex_set($A, 1,1,$alpha);
 413
 414  my $B = gsl_matrix_complex_alloc(2,2);
 415  $alpha = gsl_complex_rect(1,1);
 416  gsl_matrix_complex_set($B, 0,0,$alpha);
 417  $alpha = gsl_complex_rect(2,2);
 418  gsl_matrix_complex_set($B, 0,1,$alpha);
 419  $alpha = gsl_complex_rect(2,2);
 420  gsl_matrix_complex_set($B, 1,0,$alpha);
 421  $alpha = gsl_complex_rect(1,3);
 422  gsl_matrix_complex_set($B, 1,1,$alpha);
 423
 424  my $C = gsl_matrix_complex_alloc(2,2);
 425  $alpha = gsl_complex_rect(0,0);
 426  map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 427  map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 428
 429  $alpha = gsl_complex_rect(2,0);
 430  my $beta = gsl_complex_rect(1,0);
 431  is(gsl_blas_zsymm($CblasLeft, $CblasUpper, $alpha, $A, $B, $beta, $C),0);
 432  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [6,22]);
 433  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [-4,32]);
 434  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [-4,16]);
 435  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [-4,20]);
 436 }
 437 1;