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 Math::GSL::Test qw/:all/;
  11 use Data::Dumper;
  12 use Math::GSL::Errno qw/:all/;
  13 use strict;
  14
  15 sub make_fixture : Test(setup) {
  16 }
  17
  18 sub teardown : Test(teardown) {
  19 }
  20
  21 #sub GSL_BLAS_SDSDOT : Tests {
  22 # my $vec1 = gsl_vector_float_alloc(4);
  23 # my $vec2 = gsl_vector_float_alloc(4);
  24 # map { gsl_vector_float_set($vec1, $_, ($_+1)**2) } (0..3);
  25 # map { gsl_vector_float_set($vec2, $_, $_+1) } (0..3);
  26 # my ($x, $result)= gsl_blas_sdsdot(2, $vec1, $vec2);
  27 # this part fail because the vectors should be initiated with gsl_vector_float_alloc...
  28 # however, the gsl_vector_float_alloc function seems to be deprecated, how should I use BLAS level1 function then?
  29 # there's no test suite yet for the BLAS functions in GSL...
  30 #}
  31
  32 sub GSL_BLAS_DDOT : Tests {
  33   my $vec1 = Math::GSL::Vector->new([1,2,3,4,5]);
  34   my $vec2 = Math::GSL::Vector->new([5,4,3,2,1]);
  35   my ($x, $result) = gsl_blas_ddot($vec1->raw, $vec2->raw);
  36   ok_status($x);
  37   is($result,35);
  38 }
  39
  40 sub GSL_BLAS_ZDOTU : Tests {
  41   my $vec1 = gsl_vector_complex_alloc(2);
  42   my $vec2 = gsl_vector_complex_alloc(2);
  43   my $c = gsl_complex_rect(2,1);
  44   gsl_vector_complex_set($vec1,0,$c);
  45   gsl_vector_complex_set($vec2,0,$c);
  46   $c = gsl_complex_rect(1,1);
  47   gsl_vector_complex_set($vec1,1,$c);
  48   gsl_vector_complex_set($vec2,1,$c);
  49   ok_status(gsl_blas_zdotu($vec1, $vec2, $c));
  50   is(gsl_real($c), 3);
  51   is(gsl_imag($c), 6);
  52 }
  53
  54 sub GSL_BLAS_ZDOTC : Tests {
  55   my $vec1 = gsl_vector_complex_alloc(2);
  56   my $vec2 = gsl_vector_complex_alloc(2);
  57   my $c = gsl_complex_rect(2,1);
  58   gsl_vector_complex_set($vec1,0,$c);
  59   gsl_vector_complex_set($vec2,0,$c);
  60   $c = gsl_complex_rect(1,1);
  61   gsl_vector_complex_set($vec1,1,$c);
  62   gsl_vector_complex_set($vec2,1,$c);
  63   ok_status(gsl_blas_zdotc($vec1, $vec2, $c));
  64   is(gsl_real($c), 7);
  65   is(gsl_imag($c), 0);
  66 }
  67
  68 sub GSL_BLAS_DNRM2 : Tests {
  69   my $vec = Math::GSL::Vector->new([3,4]);
  70   is(gsl_blas_dnrm2($vec->raw), 5);
  71 }
  72
  73
  74 sub GSL_BLAS_DZNRM2 : Tests {
  75   my $vec = gsl_vector_complex_alloc(2);
  76   my $c = gsl_complex_rect(2,1);
  77   gsl_vector_complex_set($vec,0,$c);
  78   $c = gsl_complex_rect(1,1);
  79   gsl_vector_complex_set($vec,1,$c);
  80   is(gsl_blas_dznrm2($vec), sqrt(7));
  81 }
  82
  83 sub GSL_BLAS_DASUM : Tests {
  84   my $vec = Math::GSL::Vector->new([2,-3,4]);
  85   is(gsl_blas_dasum($vec->raw), 9);
  86 }
  87
  88 sub GSL_BLAS_DZASUM : Tests {
  89   my $vec = gsl_vector_complex_alloc(2);
  90   my $c = gsl_complex_rect(2,1);
  91   gsl_vector_complex_set($vec,0,$c);
  92   $c = gsl_complex_rect(1,1);
  93   gsl_vector_complex_set($vec,1,$c);
  94   is(gsl_blas_dzasum($vec), 5);
  95 }
  96
  97 sub GSL_BLAS_DSWAP : Tests {
  98   my $vec1 = Math::GSL::Vector->new([0,1,2]);
  99   my $vec2 = Math::GSL::Vector->new([2,1,0]);
 100   gsl_blas_dswap($vec1->raw, $vec2->raw);
 101   ok_similar( [0 .. 2], [ $vec2->as_list ] );
 102   ok_similar( [2, 1,0], [ $vec1->as_list ] );
 103 }
 104
 105 sub GSL_BLAS_ZSWAP : Tests {
 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   gsl_vector_complex_set($vec2,0, $c);
 114   $c = gsl_complex_rect(1,1);
 115   gsl_vector_complex_set($vec2,1, $c);
 116
 117   ok_status(gsl_blas_zswap($vec1, $vec2));
 118   $c = gsl_vector_complex_get($vec1,0);
 119   local $TODO = "Problem with the output of gsl_vector_complex_get";
 120   #is( gsl_real($c), 3);
 121   #is( gsl_imag($c), 3);
 122 }
 123
 124 sub GSL_BLAS_DCOPY : Tests {
 125  my $vec1 = Math::GSL::Vector->new([0,1,2]);
 126  my $vec2 = Math::GSL::Vector->new(3);
 127  ok_status(gsl_blas_dcopy($vec1->raw, $vec2->raw));
 128  my @got = $vec2->as_list;
 129  map { is($got[$_], $_) } (0..2);
 130 }
 131
 132 sub GSL_BLAS_DAXPY : Tests {
 133  my $vec1 = Math::GSL::Vector->new([0,1,2]);
 134  my $vec2 = Math::GSL::Vector->new([2,3,4]);
 135  ok_status(gsl_blas_daxpy(2,$vec1->raw, $vec2->raw));
 136  is_similar( [ $vec2->as_list ], [ 2, 5, 8 ] );
 137 }
 138
 139 sub GSL_BLAS_DSCAL : Tests {
 140  my $vec = Math::GSL::Vector->new([0,1,2]);
 141  gsl_blas_dscal(4, $vec->raw);
 142  my @got = $vec->as_list;
 143  map { is($got[$_], $_*4) } (0..2);
 144 }
 145
 146 sub GSL_BLAS_DROT : Tests {
 147  my $x = Math::GSL::Vector->new([1,2,3]);
 148  my $y = Math::GSL::Vector->new([0,1,2]);
 149  ok_status(gsl_blas_drot($x->raw,$y->raw,2,3));
 150  ok_similar( [$x->as_list], [ 2,7,12], 'first vector');
 151  ok_similar( [$y->as_list], [-3,-4,-5], 'second vector');
 152 }
 153
 154 sub GSL_BLAS_DGER : Tests {
 155  my $x = Math::GSL::Vector->new([1,2,3]);
 156  my $y = Math::GSL::Vector->new([0,1,2]);
 157  my $A = Math::GSL::Matrix->new(3,3);
 158  gsl_matrix_set_zero($A->raw);
 159  ok_status(gsl_blas_dger(2, $x->raw, $y->raw, $A->raw));
 160  ok_similar([$A->row(0)->as_list], [0,2,4]);
 161  ok_similar([$A->row(1)->as_list], [0,4,8]);
 162  ok_similar([$A->row(2)->as_list], [0,6,12]);
 163 }
 164
 165 sub GSL_BLAS_ZGERU : Tests {
 166  my $x = gsl_vector_complex_alloc(2);
 167  my $y = gsl_vector_complex_alloc(2);
 168  my $A = gsl_matrix_complex_alloc(2,2);
 169  my $alpha = gsl_complex_rect(2,2);
 170  gsl_vector_complex_set($x, 0, $alpha);
 171  $alpha = gsl_complex_rect(1,2);
 172  gsl_vector_complex_set($x, 1, $alpha);
 173  gsl_vector_complex_set($y, 0, $alpha);
 174  $alpha = gsl_complex_rect(3,2);
 175  gsl_vector_complex_set($y, 1, $alpha);
 176  $alpha = gsl_complex_rect(0,0);
 177  for (my $line=0; $line<2; $line++) {
 178  map { gsl_matrix_complex_set($A, $line, $_, $alpha) } (0..1); }
 179  $alpha = gsl_complex_rect(1,0);
 180  ok_status(gsl_blas_zgeru($alpha, $x, $y, $A));
 181
 182  $alpha= gsl_matrix_complex_get($A, 0,0);
 183  ok_similar([gsl_parts($alpha)], [-2, 6]);
 184  $alpha= gsl_matrix_complex_get($A, 1,0);
 185  ok_similar([gsl_parts($alpha)], [-3, 4]);
 186  $alpha= gsl_matrix_complex_get($A, 1,0);
 187  ok_similar([gsl_parts($alpha)], [-3, 4]);
 188  $alpha= gsl_matrix_complex_get($A, 0,1);
 189  ok_similar([gsl_parts($alpha)], [2, 10]);
 190  $alpha= gsl_matrix_complex_get($A, 1,1);
 191  ok_similar([gsl_parts($alpha)], [-1, 8]);
 192 }
 193
 194 sub GSL_BLAS_DGEMV : Tests {
 195  my $x = Math::GSL::Vector->new([1,2,3]);
 196  my $y = Math::GSL::Vector->new([0,1,2]);
 197  my $A = Math::GSL::Matrix->new(3,3);
 198  gsl_matrix_set_identity($A->raw);
 199  ok_status(gsl_blas_dgemv($CblasNoTrans, 2, $A->raw, $x->raw,2, $y->raw));
 200  ok_similar( [ $y->as_list ], [  2, 6, 10 ] );
 201 }
 202
 203 sub GSL_BLAS_DTRMV : Tests {
 204  my $x = Math::GSL::Vector->new([1,2,3]);
 205  my $A = Math::GSL::Matrix->new(3,3);
 206  gsl_matrix_set($A->raw, 0,0,3);
 207  gsl_matrix_set($A->raw, 1,1,3);
 208  gsl_matrix_set($A->raw, 2,2,3);
 209  gsl_matrix_set($A->raw, 0,1,2);
 210  gsl_matrix_set($A->raw, 0,2,3);
 211  gsl_matrix_set($A->raw, 1,2,4);
 212  ok_status(gsl_blas_dtrmv($CblasLower, $CblasNoTrans, $CblasNonUnit, $A->raw, $x->raw));
 213  is_similar( [ $x->as_list ], [ 3, 6, 9 ] );
 214 }
 215
 216 sub GSL_BLAS_DTRSV : Tests {
 217  my $x = Math::GSL::Vector->new([40,40,40,40]);
 218  my $A = Math::GSL::Matrix->new(4,4);
 219  map { gsl_matrix_set($A->raw, $_,0,$_+1); } (0..3);
 220  map { gsl_matrix_set($A->raw, $_,1,$_+2); } (0..2);
 221  gsl_matrix_set($A->raw, 3,1,1);
 222  map { gsl_matrix_set($A->raw, $_,2,$_+3); } (0..1);
 223  map { gsl_matrix_set($A->raw, $_,2,$_-1); } (2..3);
 224  map { gsl_matrix_set($A->raw, $_,3,$_); } (1..3);
 225  gsl_matrix_set($A->raw, 0,3,4);
 226  ok_status(gsl_blas_dtrsv($CblasLower, $CblasNoTrans, $CblasNonUnit, $A->raw, $x->raw));
 227  ok_similar([$x->as_list], [40,-40/3,-80/3,-160/9]);
 228 }
 229
 230 sub GSL_BLAS_DROTG : Tests {
 231  my $a = [1];
 232  my $b = [2];
 233
 234  my ($status, $c, $s) = gsl_blas_drotg($a, $b);
 235  ok_similar( [$c, $s ], [ 1/sqrt(5) , 2/sqrt(5)  ] );
 236 }
 237
 238 sub GSL_BLAS_DSYMV : Tests {
 239  my $x = Math::GSL::Vector->new([1,2,3]);
 240  my $y = Math::GSL::Vector->new([3,2,1]);
 241  my $A = Math::GSL::Matrix->new(3,3);
 242  map { gsl_matrix_set($A->raw, $_,0,$_+1); } (0..2);
 243  gsl_matrix_set($A->raw, 0, 1, 2);
 244  gsl_matrix_set($A->raw, 1, 1, 1);
 245  gsl_matrix_set($A->raw, 2, 1, 2);
 246  gsl_matrix_set($A->raw, 0, 2, 3);
 247  gsl_matrix_set($A->raw, 1, 2, 2);
 248  gsl_matrix_set($A->raw, 2, 2, 1);
 249  ok_status(gsl_blas_dsymv($CblasLower, 2, $A->raw, $x->raw, 3, $y->raw));
 250  ok_similar( [$y->as_list], [37,26,23]);
 251 }
 252
 253 sub GSL_BLAS_DSYR : Tests {
 254     my $x = Math::GSL::Vector->new([1,2,3]);
 255     my $A = Math::GSL::Matrix->new(3,3);
 256     gsl_matrix_set_zero($A->raw);
 257
 258     ok_status(gsl_blas_dsyr($CblasLower, 2, $x->raw, $A->raw));
 259     ok_similar([ $A->row(0)->as_list ], [2,0,0]);
 260     ok_similar([ $A->row(1)->as_list ], [4,8,0]);
 261     ok_similar([ $A->row(2)->as_list ], [6,12,18]);
 262 }
 263
 264 sub GSL_BLAS_ZHER : Tests {
 265  my $x = gsl_vector_complex_alloc(2);
 266  my $A = gsl_matrix_complex_alloc(2,2);
 267  my $alpha = gsl_complex_rect(0,0);
 268  for(my $line=0; $line<2; $line++){
 269  map { gsl_matrix_complex_set($A, $_, $line, $alpha) } (0..1); }
 270  $alpha = gsl_complex_rect(1,2);
 271  gsl_vector_complex_set($x, 0, $alpha);
 272  $alpha = gsl_complex_rect(2,2);
 273  gsl_vector_complex_set($x, 1, $alpha);
 274  $alpha = gsl_complex_rect(1,1);
 275  ok_status(gsl_blas_zher($CblasLower, 2, $x, $A));
 276  my @got = gsl_parts(gsl_matrix_complex_get($A,0,0));
 277  ok_similar([@got], [10, 0]);
 278  @got = gsl_parts(gsl_matrix_complex_get($A,1,0));
 279  ok_similar([@got], [12, -4]);
 280  @got = gsl_parts(gsl_matrix_complex_get($A,1,1));
 281  ok_similar([@got], [16, 0]);
 282  @got = gsl_parts(gsl_matrix_complex_get($A,0,1));
 283  ok_similar([@got], [0, 0]);
 284 }
 285
 286 sub GSL_BLAS_DSYR2 : Tests {
 287  my $x = Math::GSL::Vector->new([1,2,3]);
 288  my $y = Math::GSL::Vector->new([3,2,1]);
 289  my $A = Math::GSL::Matrix->new(3,3);
 290  gsl_matrix_set_zero($A->raw);
 291  map { gsl_matrix_set($A->raw, $_, 0, ($_+1)**2) } (0..2);
 292  map { gsl_matrix_set($A->raw, $_, 1, ($_+1)**2) } (1..2);
 293  gsl_matrix_set($A->raw, 0, 1, (1)**2);
 294  gsl_matrix_set($A->raw, 0, 1, 4);
 295  gsl_matrix_set($A->raw, 0, 2, 9);
 296  gsl_matrix_set($A->raw, 1, 2, 9);
 297  gsl_matrix_set($A->raw, 2, 2, 3);
 298  ok_status(gsl_blas_dsyr2($CblasLower, 2, $x->raw, $y->raw, $A->raw));
 299  my @got = $A->row(0)->as_list;
 300  ok_similar([@got], [13, 4, 9]);
 301  @got = $A->row(1)->as_list;
 302  ok_similar([@got], [20, 20, 9]);
 303  @got = $A->row(2)->as_list;
 304  ok_similar([@got], [29, 25, 15]);
 305 }
 306
 307 sub GSL_BLAS_DGEMM : Tests {
 308  my $A = Math::GSL::Matrix->new(2,2);
 309  gsl_matrix_set($A->raw, 0,0,1);
 310  gsl_matrix_set($A->raw, 1,0,3);
 311  gsl_matrix_set($A->raw, 0,1,4);
 312  gsl_matrix_set($A->raw, 1,1,2);
 313  my $B = Math::GSL::Matrix->new(2,2);
 314  gsl_matrix_set($B->raw, 0,0,2);
 315  gsl_matrix_set($B->raw, 1,0,5);
 316  gsl_matrix_set($B->raw, 0,1,1);
 317  gsl_matrix_set($B->raw, 1,1,3);
 318  my $C = Math::GSL::Matrix->new(2,2);
 319  gsl_matrix_set_zero($C->raw);
 320  ok_status(gsl_blas_dgemm($CblasNoTrans, $CblasNoTrans, 1, $A->raw, $B->raw, 1, $C->raw));
 321  my @got = $C->row(0)->as_list;
 322  ok_similar([@got], [22, 13]);
 323  @got = $C->row(1)->as_list;
 324  ok_similar([@got], [16, 9]);
 325 }
 326
 327 sub GSL_BLAS_DSYMM : Tests {
 328  my $A = Math::GSL::Matrix->new(2,2);
 329  gsl_matrix_set($A->raw, 0,0,1);
 330  gsl_matrix_set($A->raw, 1,0,3);
 331  gsl_matrix_set($A->raw, 0,1,4);
 332  gsl_matrix_set($A->raw, 1,1,2);
 333  my $B = Math::GSL::Matrix->new(2,2);
 334  gsl_matrix_set($B->raw, 0,0,2);
 335  gsl_matrix_set($B->raw, 1,0,5);
 336  gsl_matrix_set($B->raw, 0,1,1);
 337  gsl_matrix_set($B->raw, 1,1,3);
 338  my $C = Math::GSL::Matrix->new(2,2);
 339  gsl_matrix_set_zero($C->raw);
 340  ok_status(gsl_blas_dsymm($CblasLeft, $CblasUpper, 1, $A->raw, $B->raw, 1, $C->raw));
 341  my @got = $C->row(0)->as_list;
 342  ok_similar([@got], [22, 13]);
 343  @got = $C->row(1)->as_list;
 344  ok_similar([@got], [18, 10]);
 345 }
 346
 347 sub GSL_BLAS_ZGEMM : Tests {
 348  my $A = gsl_matrix_complex_alloc(2,2);
 349  my $alpha = gsl_complex_rect(1,2);
 350  gsl_matrix_complex_set($A, 0,0,$alpha);
 351  $alpha = gsl_complex_rect(3,1);
 352  gsl_matrix_complex_set($A, 0,1,$alpha);
 353  $alpha = gsl_complex_rect(2,2);
 354  gsl_matrix_complex_set($A, 1,0,$alpha);
 355  $alpha = gsl_complex_rect(0,2);
 356  gsl_matrix_complex_set($A, 1,1,$alpha);
 357
 358  my $B = gsl_matrix_complex_alloc(2,2);
 359  $alpha = gsl_complex_rect(1,1);
 360  gsl_matrix_complex_set($B, 0,0,$alpha);
 361  $alpha = gsl_complex_rect(2,2);
 362  gsl_matrix_complex_set($B, 0,1,$alpha);
 363  $alpha = gsl_complex_rect(1,2);
 364  gsl_matrix_complex_set($B, 1,0,$alpha);
 365  $alpha = gsl_complex_rect(1,3);
 366  gsl_matrix_complex_set($B, 1,1,$alpha);
 367
 368  my $C = gsl_matrix_complex_alloc(2,2);
 369  $alpha = gsl_complex_rect(0,0);
 370  map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 371  map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 372
 373  $alpha = gsl_complex_rect(2,0);
 374  my $beta = gsl_complex_rect(1,0);
 375  ok_status(gsl_blas_zgemm($CblasNoTrans, $CblasNoTrans, $alpha, $A, $B, $beta, $C));
 376  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [0,20]);
 377  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [-4,32]);
 378  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [-8,12]);
 379  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [-12,20]);
 380 }
 381
 382 sub GSL_BLAS_ZSYMM : Tests {
 383  my $A = gsl_matrix_complex_alloc(2,2);
 384  my $alpha = gsl_complex_rect(1,2);
 385  gsl_matrix_complex_set($A, 0,0,$alpha);
 386  $alpha = gsl_complex_rect(3,1);
 387  gsl_matrix_complex_set($A, 0,1,$alpha);
 388  $alpha = gsl_complex_rect(3,1);
 389  gsl_matrix_complex_set($A, 1,0,$alpha);
 390  $alpha = gsl_complex_rect(0,2);
 391  gsl_matrix_complex_set($A, 1,1,$alpha);
 392
 393  my $B = gsl_matrix_complex_alloc(2,2);
 394  $alpha = gsl_complex_rect(1,1);
 395  gsl_matrix_complex_set($B, 0,0,$alpha);
 396  $alpha = gsl_complex_rect(2,2);
 397  gsl_matrix_complex_set($B, 0,1,$alpha);
 398  $alpha = gsl_complex_rect(2,2);
 399  gsl_matrix_complex_set($B, 1,0,$alpha);
 400  $alpha = gsl_complex_rect(1,3);
 401  gsl_matrix_complex_set($B, 1,1,$alpha);
 402
 403  my $C = gsl_matrix_complex_alloc(2,2);
 404  $alpha = gsl_complex_rect(0,0);
 405  map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 406  map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 407
 408  $alpha = gsl_complex_rect(2,0);
 409  my $beta = gsl_complex_rect(1,0);
 410  ok_status(gsl_blas_zsymm($CblasLeft, $CblasUpper, $alpha, $A, $B, $beta, $C));
 411  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [6,22]);
 412  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [-4,32]);
 413  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [-4,16]);
 414  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [-4,20]);
 415 }
 416 sub GSL_BLAS_ZHEMM : Tests {
 417  my $A = gsl_matrix_complex_alloc(2,2);
 418  my $alpha = gsl_complex_rect(3,0);
 419  gsl_matrix_complex_set($A, 0,0,$alpha);
 420  $alpha = gsl_complex_rect(2,1);
 421  gsl_matrix_complex_set($A, 0,1,$alpha);
 422  $alpha = gsl_complex_rect(2,-1);
 423  gsl_matrix_complex_set($A, 1,0,$alpha);
 424  $alpha = gsl_complex_rect(1,0);
 425  gsl_matrix_complex_set($A, 1,1,$alpha);
 426
 427  my $B = gsl_matrix_complex_alloc(2,2);
 428  $alpha = gsl_complex_rect(1,0);
 429  gsl_matrix_complex_set($B, 0,0,$alpha);
 430  $alpha = gsl_complex_rect(2,2);
 431  gsl_matrix_complex_set($B, 0,1,$alpha);
 432  $alpha = gsl_complex_rect(2,-2);
 433  gsl_matrix_complex_set($B, 1,0,$alpha);
 434  $alpha = gsl_complex_rect(2,0);
 435  gsl_matrix_complex_set($B, 1,1,$alpha);
 436
 437  my $C = gsl_matrix_complex_alloc(2,2);
 438  $alpha = gsl_complex_rect(0,0);
 439  map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 440  map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 441
 442  $alpha = gsl_complex_rect(2,0);
 443  my $beta = gsl_complex_rect(1,0);
 444  ok_status(gsl_blas_zhemm($CblasLeft, $CblasUpper, $alpha, $A, $B, $beta, $C));
 445  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [18,-4]);
 446  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [20,16]);
 447  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [8,-6]);
 448  ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [16,4]);
 449 }
 450
 451 sub GSL_BLAS_DTRMM : Tests {
 452  my $A = Math::GSL::Matrix->new(2,2);
 453  gsl_matrix_set($A->raw, 0,0,1);
 454  gsl_matrix_set($A->raw, 1,0,3);
 455  gsl_matrix_set($A->raw, 0,1,4);
 456  gsl_matrix_set($A->raw, 1,1,2);
 457  my $B = Math::GSL::Matrix->new(2,2);
 458  gsl_matrix_set($B->raw, 0,0,2);
 459  gsl_matrix_set($B->raw, 1,0,5);
 460  gsl_matrix_set($B->raw, 0,1,1);
 461  gsl_matrix_set($B->raw, 1,1,3);
 462  my $C = Math::GSL::Matrix->new(2,2);
 463  gsl_matrix_set_zero($C->raw);
 464  ok_status(gsl_blas_dtrmm($CblasLeft, $CblasUpper, $CblasNoTrans, $CblasUnit, 1, $A->raw, $B->raw));
 465  my @got = $B->row(0)->as_list;
 466  ok_similar([@got], [22, 13]);
 467  @got = $B->row(1)->as_list;
 468  ok_similar([@got], [5, 3]);
 469 }
 470
 471 sub GSL_BLAS_ZTRMM : Tests {
 472  my $A = gsl_matrix_complex_alloc(2,2);
 473  my $alpha = gsl_complex_rect(3,1);
 474  gsl_matrix_complex_set($A, 0,1,$alpha);
 475
 476  my $B = gsl_matrix_complex_alloc(2,2);
 477  $alpha = gsl_complex_rect(1,0);
 478  gsl_matrix_complex_set($B, 0,0,$alpha);
 479  $alpha = gsl_complex_rect(2,2);
 480  gsl_matrix_complex_set($B, 0,1,$alpha);
 481  $alpha = gsl_complex_rect(1,-2);
 482  gsl_matrix_complex_set($B, 1,0,$alpha);
 483  $alpha = gsl_complex_rect(4,2);
 484  gsl_matrix_complex_set($B, 1,1,$alpha);
 485
 486  $alpha = gsl_complex_rect(1,0);
 487  ok_status(gsl_blas_ztrmm($CblasLeft, $CblasUpper, $CblasNoTrans, $CblasUnit, $alpha, $A, $B));
 488  ok_similar([gsl_parts(gsl_matrix_complex_get($B, 0, 0))], [6, -5]);
 489  ok_similar([gsl_parts(gsl_matrix_complex_get($B, 0, 1))], [12, 12]);
 490  ok_similar([gsl_parts(gsl_matrix_complex_get($B, 1, 0))], [1, -2]);
 491  ok_similar([gsl_parts(gsl_matrix_complex_get($B, 1, 1))], [4, 2]);
 492 }
 493
 494 sub GSL_BLAS_DSYRK : Tests {
 495  my $A = Math::GSL::Matrix->new(2,2);
 496  gsl_matrix_set($A->raw, 0, 0, 1);
 497  gsl_matrix_set($A->raw, 0, 1, 4);
 498  gsl_matrix_set($A->raw, 1, 0, 4);
 499  gsl_matrix_set($A->raw, 1, 1, 3);
 500  my $C = Math::GSL::Matrix->new(2,2);
 501  gsl_matrix_set_zero($C->raw);
 502  ok_status(gsl_blas_dsyrk ($CblasUpper, $CblasNoTrans, 1, $A->raw, 1, $C->raw));
 503  ok_similar([$C->row(0)->as_list], [17,16]);
 504  ok_similar([$C->row(1)->as_list], [0,25]);
 505 }
 506
 507 sub GSL_BLAS_ZSYRK : Tests {
 508   my $A = gsl_matrix_complex_alloc(2,2);
 509  my $alpha = gsl_complex_rect(3,1);
 510  gsl_matrix_complex_set($A, 0,0,$alpha);
 511  $alpha = gsl_complex_rect(2,1);
 512  gsl_matrix_complex_set($A, 0,1,$alpha);
 513  gsl_matrix_complex_set($A, 1,0,$alpha);
 514  $alpha = gsl_complex_rect(1,1);
 515  gsl_matrix_complex_set($A, 1,1,$alpha);
 516
 517  my $C = gsl_matrix_complex_alloc(2,2);
 518  $alpha = gsl_complex_rect(0,0);
 519  map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 520  map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 521
 522  $alpha = gsl_complex_rect(1,0);
 523  my $beta = gsl_complex_rect(1,0);
 524  ok_status(gsl_blas_zsyrk($CblasUpper, $CblasNoTrans, $alpha, $A, $beta, $C));
 525  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [11, 10]);
 526  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [6, 8]);
 527  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 528  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [3, 6]);
 529 }
 530
 531 sub GSL_BLAS_ZHERK : Tests {
 532   my $A = gsl_matrix_complex_alloc(2,2);
 533  my $alpha = gsl_complex_rect(3,0);
 534  gsl_matrix_complex_set($A, 0,0,$alpha);
 535  $alpha = gsl_complex_rect(2,1);
 536  gsl_matrix_complex_set($A, 0,1,$alpha);
 537  $alpha = gsl_complex_rect(2,-1);
 538  gsl_matrix_complex_set($A, 1,0,$alpha);
 539  $alpha = gsl_complex_rect(1,0);
 540  gsl_matrix_complex_set($A, 1,1,$alpha);
 541
 542  my $C = gsl_matrix_complex_alloc(2,2);
 543  $alpha = gsl_complex_rect(0,0);
 544  map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 545  map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 546
 547  ok_status(gsl_blas_zherk ($CblasUpper, $CblasNoTrans, 1, $A, 1, $C));
 548  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [14, 0]);
 549  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [8, 4]);
 550  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 551  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [6, 0]);
 552 }
 553
 554 sub GSL_BLAS_ZHER2K : Tests {
 555  my $A = gsl_matrix_complex_alloc(2,2);
 556  my $alpha = gsl_complex_rect(3,0);
 557  gsl_matrix_complex_set($A, 0,0,$alpha);
 558  $alpha = gsl_complex_rect(2,1);
 559  gsl_matrix_complex_set($A, 0,1,$alpha);
 560  $alpha = gsl_complex_rect(2,-1);
 561  gsl_matrix_complex_set($A, 1,0,$alpha);
 562  $alpha = gsl_complex_rect(1,0);
 563  gsl_matrix_complex_set($A, 1,1,$alpha);
 564
 565  my $B = gsl_matrix_complex_alloc(2,2);
 566  $alpha = gsl_complex_rect(6,0);
 567  gsl_matrix_complex_set($A, 0,0,$alpha);
 568  $alpha = gsl_complex_rect(3,1);
 569  gsl_matrix_complex_set($A, 0,1,$alpha);
 570  $alpha = gsl_complex_rect(3,-1);
 571  gsl_matrix_complex_set($A, 1,0,$alpha);
 572  $alpha = gsl_complex_rect(5,0);
 573  gsl_matrix_complex_set($A, 1,1,$alpha);
 574
 575  my $C = gsl_matrix_complex_alloc(2,2);
 576  $alpha = gsl_complex_rect(0,0);
 577  map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 578  map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 579
 580  $alpha = gsl_complex_rect(1,0);
 581
 582  ok_status(gsl_blas_zher2k($CblasUpper, $CblasNoTrans, $alpha, $A, $B, 1, $C));
 583  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 584  local $TODO = "These results follow the formula given by the documentation, don't know why it fails";
 585  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [50, 0]);
 586  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [34, 15]);
 587  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [24, 0]);
 588 }
 589
 590 sub GSL_BLAS_DSYR2K : Tests {
 591  my $A = Math::GSL::Matrix->new(2,2);
 592  gsl_matrix_set($A->raw, 0, 0, 1);
 593  gsl_matrix_set($A->raw, 0, 1, 4);
 594  gsl_matrix_set($A->raw, 1, 0, 4);
 595  gsl_matrix_set($A->raw, 1, 1, 3);
 596  my $B = Math::GSL::Matrix->new(2,2);
 597  gsl_matrix_set($B->raw, 0, 0, 2);
 598  gsl_matrix_set($B->raw, 0, 1, 5);
 599  gsl_matrix_set($B->raw, 1, 0, 5);
 600  gsl_matrix_set($B->raw, 1, 1, 1);
 601  my $C = Math::GSL::Matrix->new(2,2);
 602  gsl_matrix_set_zero($C->raw);
 603  ok_status(gsl_blas_dsyr2k ($CblasUpper, $CblasNoTrans, 1, $A->raw, $B->raw, 1, $C->raw, ));
 604  ok_similar([$C->row(0)->as_list], [44,32]);
 605  ok_similar([$C->row(1)->as_list], [0,46]);
 606 }
 607
 608 sub GSL_BLAS_ZSYR2K : Tests {
 609  my $A = gsl_matrix_complex_alloc(2,2);
 610  my $alpha = gsl_complex_rect(3,0);
 611  gsl_matrix_complex_set($A, 0,0,$alpha);
 612  $alpha = gsl_complex_rect(2,1);
 613  gsl_matrix_complex_set($A, 0,1,$alpha);
 614  $alpha = gsl_complex_rect(2,1);
 615  gsl_matrix_complex_set($A, 1,0,$alpha);
 616  $alpha = gsl_complex_rect(1,0);
 617  gsl_matrix_complex_set($A, 1,1,$alpha);
 618
 619  my $B = gsl_matrix_complex_alloc(2,2);
 620  $alpha = gsl_complex_rect(6,0);
 621  gsl_matrix_complex_set($A, 0,0,$alpha);
 622  $alpha = gsl_complex_rect(3,1);
 623  gsl_matrix_complex_set($A, 0,1,$alpha);
 624  $alpha = gsl_complex_rect(3,1);
 625  gsl_matrix_complex_set($A, 1,0,$alpha);
 626  $alpha = gsl_complex_rect(5,0);
 627  gsl_matrix_complex_set($A, 1,1,$alpha);
 628
 629  my $C = gsl_matrix_complex_alloc(2,2);
 630  $alpha = gsl_complex_rect(0,0);
 631  map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 632  map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 633
 634  $alpha = gsl_complex_rect(1,0);
 635  my $beta = gsl_complex_rect(1,0);
 636
 637  ok_status(gsl_blas_zsyr2k($CblasUpper, $CblasNoTrans, $alpha, $A, $B, $beta, $C));
 638  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 639  local $TODO = "These results follow the formula given by the documentation, don't know why it fails";
 640  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [46, 10]);
 641  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [34, 15]);
 642  ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [24, 4]);
 643
 644
 645 }
 646 1;