t/BLAS.t

   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 BEGIN { gsl_set_error_handler_off(); }
  16
  17 sub make_fixture : Test(setup) {
  18 }
  19
  20 sub teardown : Test(teardown) {
  21 }
  22
  23 #sub GSL_BLAS_SDSDOT : Tests {
  24 # my $vec1 = gsl_vector_float_alloc(4);
  25 # my $vec2 = gsl_vector_float_alloc(4);
  26 # map { gsl_vector_float_set($vec1, $_, ($_+1)**2) } (0..3);
  27 # map { gsl_vector_float_set($vec2, $_, $_+1) } (0..3);
  28 # my ($x, $result)= gsl_blas_sdsdot(2, $vec1, $vec2);
  29 # this part fail because the vectors should be initiated with gsl_vector_float_alloc...
  30 # however, the gsl_vector_float_alloc function seems to be deprecated, how should I use BLAS level1 function then?
  31 # there's no test suite yet for the BLAS functions in GSL...
  32 #}
  33
  34 sub GSL_BLAS_DDOT : Tests {
  35     my $vec1 = Math::GSL::Vector->new([1,2,3,4,5]);
  36     my $vec2 = Math::GSL::Vector->new([5,4,3,2,1]);
  37     my ($x, $result) = gsl_blas_ddot($vec1->raw, $vec2->raw);
  38     ok_status($x);
  39     is_similar($result,35);
  40 }
  41
  42 sub GSL_BLAS_ZDOTU : Tests {
  43     my $vec1 = gsl_vector_complex_alloc(2);
  44     my $vec2 = gsl_vector_complex_alloc(2);
  45     my $c = gsl_complex_rect(2,1);
  46     gsl_vector_complex_set($vec1,0,$c);
  47     gsl_vector_complex_set($vec2,0,$c);
  48     $c = gsl_complex_rect(1,1);
  49     gsl_vector_complex_set($vec1,1,$c);
  50     gsl_vector_complex_set($vec2,1,$c);
  51     ok_status(gsl_blas_zdotu($vec1, $vec2, $c));
  52     is_similar([ gsl_parts($c) ], [3,6]);
  53 }
  54
  55 sub GSL_BLAS_ZDOTC : Tests {
  56     my $vec1 = gsl_vector_complex_alloc(2);
  57     my $vec2 = gsl_vector_complex_alloc(2);
  58     my $c = gsl_complex_rect(2,1);
  59     gsl_vector_complex_set($vec1,0,$c);
  60     gsl_vector_complex_set($vec2,0,$c);
  61     $c = gsl_complex_rect(1,1);
  62     gsl_vector_complex_set($vec1,1,$c);
  63     gsl_vector_complex_set($vec2,1,$c);
  64     ok_status(gsl_blas_zdotc($vec1, $vec2, $c));
  65     is_similar([gsl_parts($c)], [7,0]);
  66 }
  67
  68 sub GSL_BLAS_DNRM2 : Tests {
  69     my $vec = Math::GSL::Vector->new([3,4]);
  70     is_similar(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_similar([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_similar(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_similar(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     ok( defined $c,"gsl_vector_complex_get");
 121
 122     # goes boom
 123     #is_similar( [gsl_parts($c) ], [ 3,3 ] );
 124 }
 125
 126 sub GSL_BLAS_DCOPY : Tests {
 127     my $vec1 = Math::GSL::Vector->new([0,1,2]);
 128     my $vec2 = Math::GSL::Vector->new(3);
 129     ok_status(gsl_blas_dcopy($vec1->raw, $vec2->raw));
 130     ok_similar( [ $vec2->as_list ], [ 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     ok_status(gsl_blas_daxpy(2,$vec1->raw, $vec2->raw));
 137     is_similar( [ $vec2->as_list ], [ 2, 5, 8 ] );
 138 }
 139
 140 sub GSL_BLAS_DSCAL : Tests {
 141     my $vec = Math::GSL::Vector->new([0,1,2]);
 142     gsl_blas_dscal(4, $vec->raw);
 143     is_similar( [ $vec->as_list ], [0,4,8] );
 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
 269     for my $line (0 .. 1) {
 270         map { gsl_matrix_complex_set($A, $_, $line, $alpha) } (0..1);
 271     }
 272     $alpha = gsl_complex_rect(1,2);
 273     gsl_vector_complex_set($x, 0, $alpha);
 274     $alpha = gsl_complex_rect(2,2);
 275     gsl_vector_complex_set($x, 1, $alpha);
 276     $alpha = gsl_complex_rect(1,1);
 277     ok_status(gsl_blas_zher($CblasLower, 2, $x, $A));
 278
 279     ok_similar([ gsl_parts(gsl_matrix_complex_get($A,0,0)) ], [10, 0 ]);
 280     ok_similar([ gsl_parts(gsl_matrix_complex_get($A,1,0)) ], [12, -4]);
 281     ok_similar([ gsl_parts(gsl_matrix_complex_get($A,1,1)) ], [16, 0 ]);
 282     ok_similar([ gsl_parts(gsl_matrix_complex_get($A,0,1)) ], [0 , 0 ]);
 283 }
 284
 285 sub GSL_BLAS_DSYR2 : Tests {
 286     my $x = Math::GSL::Vector->new([1,2,3]);
 287     my $y = Math::GSL::Vector->new([3,2,1]);
 288     my $A = Math::GSL::Matrix->new(3,3);
 289     gsl_matrix_set_zero($A->raw);
 290     map { gsl_matrix_set($A->raw, $_, 0, ($_+1)**2) } (0..2);
 291     map { gsl_matrix_set($A->raw, $_, 1, ($_+1)**2) } (1..2);
 292     gsl_matrix_set($A->raw, 0, 1, (1)**2);
 293     gsl_matrix_set($A->raw, 0, 1, 4);
 294     gsl_matrix_set($A->raw, 0, 2, 9);
 295     gsl_matrix_set($A->raw, 1, 2, 9);
 296     gsl_matrix_set($A->raw, 2, 2, 3);
 297     ok_status(gsl_blas_dsyr2($CblasLower, 2, $x->raw, $y->raw, $A->raw));
 298     my @got = $A->row(0)->as_list;
 299     ok_similar([@got], [13, 4, 9]);
 300     @got = $A->row(1)->as_list;
 301     ok_similar([@got], [20, 20, 9]);
 302     @got = $A->row(2)->as_list;
 303     ok_similar([@got], [29, 25, 15]);
 304 }
 305
 306 sub GSL_BLAS_DGEMM : Tests {
 307     my $A = Math::GSL::Matrix->new(2,2);
 308     gsl_matrix_set($A->raw, 0,0,1);
 309     gsl_matrix_set($A->raw, 1,0,3);
 310     gsl_matrix_set($A->raw, 0,1,4);
 311     gsl_matrix_set($A->raw, 1,1,2);
 312     my $B = Math::GSL::Matrix->new(2,2);
 313     gsl_matrix_set($B->raw, 0,0,2);
 314     gsl_matrix_set($B->raw, 1,0,5);
 315     gsl_matrix_set($B->raw, 0,1,1);
 316     gsl_matrix_set($B->raw, 1,1,3);
 317     my $C = Math::GSL::Matrix->new(2,2);
 318     gsl_matrix_set_zero($C->raw);
 319     ok_status(gsl_blas_dgemm($CblasNoTrans, $CblasNoTrans, 1, $A->raw, $B->raw, 1, $C->raw));
 320     my @got = $C->row(0)->as_list;
 321     ok_similar([@got], [22, 13]);
 322     @got = $C->row(1)->as_list;
 323     ok_similar([@got], [16, 9]);
 324 }
 325
 326 sub GSL_BLAS_DSYMM : Tests {
 327     my $A = Math::GSL::Matrix->new(2,2);
 328     gsl_matrix_set($A->raw, 0,0,1);
 329     gsl_matrix_set($A->raw, 1,0,3);
 330     gsl_matrix_set($A->raw, 0,1,4);
 331     gsl_matrix_set($A->raw, 1,1,2);
 332     my $B = Math::GSL::Matrix->new(2,2);
 333     gsl_matrix_set($B->raw, 0,0,2);
 334     gsl_matrix_set($B->raw, 1,0,5);
 335     gsl_matrix_set($B->raw, 0,1,1);
 336     gsl_matrix_set($B->raw, 1,1,3);
 337     my $C = Math::GSL::Matrix->new(2,2);
 338     gsl_matrix_set_zero($C->raw);
 339     ok_status(gsl_blas_dsymm($CblasLeft, $CblasUpper, 1, $A->raw, $B->raw, 1, $C->raw));
 340     my @got = $C->row(0)->as_list;
 341     ok_similar([@got], [22, 13]);
 342     @got = $C->row(1)->as_list;
 343     ok_similar([@got], [18, 10]);
 344 }
 345
 346 sub GSL_BLAS_ZGEMM : Tests {
 347     my $A = gsl_matrix_complex_alloc(2,2);
 348     my $alpha = gsl_complex_rect(1,2);
 349     gsl_matrix_complex_set($A, 0,0,$alpha);
 350     $alpha = gsl_complex_rect(3,1);
 351     gsl_matrix_complex_set($A, 0,1,$alpha);
 352     $alpha = gsl_complex_rect(2,2);
 353     gsl_matrix_complex_set($A, 1,0,$alpha);
 354     $alpha = gsl_complex_rect(0,2);
 355     gsl_matrix_complex_set($A, 1,1,$alpha);
 356
 357     my $B = gsl_matrix_complex_alloc(2,2);
 358     $alpha = gsl_complex_rect(1,1);
 359     gsl_matrix_complex_set($B, 0,0,$alpha);
 360     $alpha = gsl_complex_rect(2,2);
 361     gsl_matrix_complex_set($B, 0,1,$alpha);
 362     $alpha = gsl_complex_rect(1,2);
 363     gsl_matrix_complex_set($B, 1,0,$alpha);
 364     $alpha = gsl_complex_rect(1,3);
 365     gsl_matrix_complex_set($B, 1,1,$alpha);
 366
 367     my $C = gsl_matrix_complex_alloc(2,2);
 368     $alpha = gsl_complex_rect(0,0);
 369     map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 370     map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 371
 372     $alpha = gsl_complex_rect(2,0);
 373     my $beta = gsl_complex_rect(1,0);
 374     ok_status(gsl_blas_zgemm($CblasNoTrans, $CblasNoTrans, $alpha, $A, $B, $beta, $C));
 375     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [0,20]);
 376     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [-4,32]);
 377     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [-8,12]);
 378     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [-12,20]);
 379 }
 380
 381 sub GSL_BLAS_ZSYMM : Tests {
 382     my $A = gsl_matrix_complex_alloc(2,2);
 383     my $alpha = gsl_complex_rect(1,2);
 384     gsl_matrix_complex_set($A, 0,0,$alpha);
 385     $alpha = gsl_complex_rect(3,1);
 386     gsl_matrix_complex_set($A, 0,1,$alpha);
 387     $alpha = gsl_complex_rect(3,1);
 388     gsl_matrix_complex_set($A, 1,0,$alpha);
 389     $alpha = gsl_complex_rect(0,2);
 390     gsl_matrix_complex_set($A, 1,1,$alpha);
 391
 392     my $B = gsl_matrix_complex_alloc(2,2);
 393     $alpha = gsl_complex_rect(1,1);
 394     gsl_matrix_complex_set($B, 0,0,$alpha);
 395     $alpha = gsl_complex_rect(2,2);
 396     gsl_matrix_complex_set($B, 0,1,$alpha);
 397     $alpha = gsl_complex_rect(2,2);
 398     gsl_matrix_complex_set($B, 1,0,$alpha);
 399     $alpha = gsl_complex_rect(1,3);
 400     gsl_matrix_complex_set($B, 1,1,$alpha);
 401
 402     my $C = gsl_matrix_complex_alloc(2,2);
 403     $alpha = gsl_complex_rect(0,0);
 404     map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 405     map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 406
 407     $alpha = gsl_complex_rect(2,0);
 408     my $beta = gsl_complex_rect(1,0);
 409     ok_status(gsl_blas_zsymm($CblasLeft, $CblasUpper, $alpha, $A, $B, $beta, $C));
 410     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [6,22]);
 411     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [-4,32]);
 412     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [-4,16]);
 413     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [-4,20]);
 414 }
 415 sub GSL_BLAS_ZHEMM : Tests {
 416     my $A = gsl_matrix_complex_alloc(2,2);
 417     my $alpha = gsl_complex_rect(3,0);
 418     gsl_matrix_complex_set($A, 0,0,$alpha);
 419     $alpha = gsl_complex_rect(2,1);
 420     gsl_matrix_complex_set($A, 0,1,$alpha);
 421     $alpha = gsl_complex_rect(2,-1);
 422     gsl_matrix_complex_set($A, 1,0,$alpha);
 423     $alpha = gsl_complex_rect(1,0);
 424     gsl_matrix_complex_set($A, 1,1,$alpha);
 425
 426     my $B = gsl_matrix_complex_alloc(2,2);
 427     $alpha = gsl_complex_rect(1,0);
 428     gsl_matrix_complex_set($B, 0,0,$alpha);
 429     $alpha = gsl_complex_rect(2,2);
 430     gsl_matrix_complex_set($B, 0,1,$alpha);
 431     $alpha = gsl_complex_rect(2,-2);
 432     gsl_matrix_complex_set($B, 1,0,$alpha);
 433     $alpha = gsl_complex_rect(2,0);
 434     gsl_matrix_complex_set($B, 1,1,$alpha);
 435
 436     my $C = gsl_matrix_complex_alloc(2,2);
 437     $alpha = gsl_complex_rect(0,0);
 438     map { gsl_matrix_complex_set($C, 0, $_, $alpha) } (0..1);
 439     map { gsl_matrix_complex_set($C, 1, $_, $alpha) } (0..1);
 440
 441     $alpha = gsl_complex_rect(2,0);
 442     my $beta = gsl_complex_rect(1,0);
 443     ok_status(gsl_blas_zhemm($CblasLeft, $CblasUpper, $alpha, $A, $B, $beta, $C));
 444     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [18,-4]);
 445     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [20,16]);
 446     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [8,-6]);
 447     ok_similar([ gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [16,4]);
 448 }
 449
 450 sub GSL_BLAS_DTRMM : Tests {
 451     my $A = Math::GSL::Matrix->new(2,2);
 452     gsl_matrix_set($A->raw, 0,0,1);
 453     gsl_matrix_set($A->raw, 1,0,3);
 454     gsl_matrix_set($A->raw, 0,1,4);
 455     gsl_matrix_set($A->raw, 1,1,2);
 456     my $B = Math::GSL::Matrix->new(2,2);
 457     gsl_matrix_set($B->raw, 0,0,2);
 458     gsl_matrix_set($B->raw, 1,0,5);
 459     gsl_matrix_set($B->raw, 0,1,1);
 460     gsl_matrix_set($B->raw, 1,1,3);
 461     my $C = Math::GSL::Matrix->new(2,2);
 462     gsl_matrix_set_zero($C->raw);
 463     ok_status(gsl_blas_dtrmm($CblasLeft, $CblasUpper, $CblasNoTrans, $CblasUnit, 1, $A->raw, $B->raw));
 464     my @got = $B->row(0)->as_list;
 465     ok_similar([@got], [22, 13]);
 466     @got = $B->row(1)->as_list;
 467     ok_similar([@got], [5, 3]);
 468 }
 469
 470 sub GSL_BLAS_ZTRMM : Tests {
 471     my $A = gsl_matrix_complex_alloc(2,2);
 472     my $alpha = gsl_complex_rect(3,1);
 473     gsl_matrix_complex_set($A, 0,1,$alpha);
 474
 475     my $B = gsl_matrix_complex_alloc(2,2);
 476     $alpha = gsl_complex_rect(1,0);
 477     gsl_matrix_complex_set($B, 0,0,$alpha);
 478     $alpha = gsl_complex_rect(2,2);
 479     gsl_matrix_complex_set($B, 0,1,$alpha);
 480     $alpha = gsl_complex_rect(1,-2);
 481     gsl_matrix_complex_set($B, 1,0,$alpha);
 482     $alpha = gsl_complex_rect(4,2);
 483     gsl_matrix_complex_set($B, 1,1,$alpha);
 484
 485     $alpha = gsl_complex_rect(1,0);
 486     ok_status(gsl_blas_ztrmm($CblasLeft, $CblasUpper, $CblasNoTrans, $CblasUnit, $alpha, $A, $B));
 487     ok_similar([gsl_parts(gsl_matrix_complex_get($B, 0, 0))], [6, -5]);
 488     ok_similar([gsl_parts(gsl_matrix_complex_get($B, 0, 1))], [12, 12]);
 489     ok_similar([gsl_parts(gsl_matrix_complex_get($B, 1, 0))], [1, -2]);
 490     ok_similar([gsl_parts(gsl_matrix_complex_get($B, 1, 1))], [4, 2]);
 491 }
 492
 493 sub GSL_BLAS_DSYRK : Tests {
 494     my $A = Math::GSL::Matrix->new(2,2);
 495     gsl_matrix_set($A->raw, 0, 0, 1);
 496     gsl_matrix_set($A->raw, 0, 1, 4);
 497     gsl_matrix_set($A->raw, 1, 0, 4);
 498     gsl_matrix_set($A->raw, 1, 1, 3);
 499     my $C = Math::GSL::Matrix->new(2,2);
 500     gsl_matrix_set_zero($C->raw);
 501     ok_status(gsl_blas_dsyrk ($CblasUpper, $CblasNoTrans, 1, $A->raw, 1, $C->raw));
 502     ok_similar([$C->row(0)->as_list], [17,16]);
 503     ok_similar([$C->row(1)->as_list], [0,25]);
 504 }
 505
 506 sub GSL_BLAS_ZSYRK : Tests {
 507     my $A = gsl_matrix_complex_alloc(2,2);
 508     my $alpha = gsl_complex_rect(3,1);
 509     gsl_matrix_complex_set($A, 0,0,$alpha);
 510     $alpha = gsl_complex_rect(2,1);
 511     gsl_matrix_complex_set($A, 0,1,$alpha);
 512     gsl_matrix_complex_set($A, 1,0,$alpha);
 513     $alpha = gsl_complex_rect(1,1);
 514     gsl_matrix_complex_set($A, 1,1,$alpha);
 515
 516     my $C = gsl_matrix_complex_alloc(2,2);
 517     $alpha = gsl_complex_rect(0,0);
 518     map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 519     map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 520
 521     $alpha = gsl_complex_rect(1,0);
 522     my $beta = gsl_complex_rect(1,0);
 523     ok_status(gsl_blas_zsyrk($CblasUpper, $CblasNoTrans, $alpha, $A, $beta, $C));
 524     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [11, 10]);
 525     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [6, 8]);
 526     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 527     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [3, 6]);
 528 }
 529
 530 sub GSL_BLAS_ZHERK : Tests {
 531     my $A = gsl_matrix_complex_alloc(2,2);
 532     my $alpha = gsl_complex_rect(3,0);
 533     gsl_matrix_complex_set($A, 0,0,$alpha);
 534     $alpha = gsl_complex_rect(2,1);
 535     gsl_matrix_complex_set($A, 0,1,$alpha);
 536     $alpha = gsl_complex_rect(2,-1);
 537     gsl_matrix_complex_set($A, 1,0,$alpha);
 538     $alpha = gsl_complex_rect(1,0);
 539     gsl_matrix_complex_set($A, 1,1,$alpha);
 540
 541     my $C = gsl_matrix_complex_alloc(2,2);
 542     $alpha = gsl_complex_rect(0,0);
 543     map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 544     map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 545
 546     ok_status(gsl_blas_zherk ($CblasUpper, $CblasNoTrans, 1, $A, 1, $C));
 547     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [14, 0]);
 548     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [8, 4]);
 549     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 550     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [6, 0]);
 551 }
 552
 553 sub GSL_BLAS_ZHER2K : Tests {
 554     my $A = gsl_matrix_complex_alloc(2,2);
 555     my $alpha = gsl_complex_rect(3,0);
 556     gsl_matrix_complex_set($A, 0,0,$alpha);
 557     $alpha = gsl_complex_rect(2,1);
 558     gsl_matrix_complex_set($A, 0,1,$alpha);
 559     $alpha = gsl_complex_rect(2,-1);
 560     gsl_matrix_complex_set($A, 1,0,$alpha);
 561     $alpha = gsl_complex_rect(1,0);
 562     gsl_matrix_complex_set($A, 1,1,$alpha);
 563
 564     my $B = gsl_matrix_complex_alloc(2,2);
 565     $alpha = gsl_complex_rect(6,0);
 566     gsl_matrix_complex_set($A, 0,0,$alpha);
 567     $alpha = gsl_complex_rect(3,1);
 568     gsl_matrix_complex_set($A, 0,1,$alpha);
 569     $alpha = gsl_complex_rect(3,-1);
 570     gsl_matrix_complex_set($A, 1,0,$alpha);
 571     $alpha = gsl_complex_rect(5,0);
 572     gsl_matrix_complex_set($A, 1,1,$alpha);
 573
 574     my $C = gsl_matrix_complex_alloc(2,2);
 575     $alpha = gsl_complex_rect(0,0);
 576     map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 577     map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 578
 579     $alpha = gsl_complex_rect(1,0);
 580
 581     ok_status(gsl_blas_zher2k($CblasUpper, $CblasNoTrans, $alpha, $A, $B, 1, $C));
 582     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 583     local $TODO = "These results follow the formula given by the documentation, don't know why it fails";
 584     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [50, 0]);
 585     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [34, 15]);
 586     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [24, 0]);
 587 }
 588
 589 sub GSL_BLAS_DSYR2K : Tests {
 590     my $A = Math::GSL::Matrix->new(2,2);
 591     gsl_matrix_set($A->raw, 0, 0, 1);
 592     gsl_matrix_set($A->raw, 0, 1, 4);
 593     gsl_matrix_set($A->raw, 1, 0, 4);
 594     gsl_matrix_set($A->raw, 1, 1, 3);
 595     my $B = Math::GSL::Matrix->new(2,2);
 596     gsl_matrix_set($B->raw, 0, 0, 2);
 597     gsl_matrix_set($B->raw, 0, 1, 5);
 598     gsl_matrix_set($B->raw, 1, 0, 5);
 599     gsl_matrix_set($B->raw, 1, 1, 1);
 600     my $C = Math::GSL::Matrix->new(2,2);
 601     gsl_matrix_set_zero($C->raw);
 602     ok_status(gsl_blas_dsyr2k ($CblasUpper, $CblasNoTrans, 1, $A->raw, $B->raw, 1, $C->raw, ));
 603     ok_similar([$C->row(0)->as_list], [44,32]);
 604     ok_similar([$C->row(1)->as_list], [0,46]);
 605 }
 606
 607 sub GSL_BLAS_ZSYR2K : Tests {
 608     my $A = gsl_matrix_complex_alloc(2,2);
 609     my $alpha = gsl_complex_rect(3,0);
 610     gsl_matrix_complex_set($A, 0,0,$alpha);
 611     $alpha = gsl_complex_rect(2,1);
 612     gsl_matrix_complex_set($A, 0,1,$alpha);
 613     $alpha = gsl_complex_rect(2,1);
 614     gsl_matrix_complex_set($A, 1,0,$alpha);
 615     $alpha = gsl_complex_rect(1,0);
 616     gsl_matrix_complex_set($A, 1,1,$alpha);
 617
 618     my $B = gsl_matrix_complex_alloc(2,2);
 619     $alpha = gsl_complex_rect(6,0);
 620     gsl_matrix_complex_set($A, 0,0,$alpha);
 621     $alpha = gsl_complex_rect(3,1);
 622     gsl_matrix_complex_set($A, 0,1,$alpha);
 623     $alpha = gsl_complex_rect(3,1);
 624     gsl_matrix_complex_set($A, 1,0,$alpha);
 625     $alpha = gsl_complex_rect(5,0);
 626     gsl_matrix_complex_set($A, 1,1,$alpha);
 627
 628     my $C = gsl_matrix_complex_alloc(2,2);
 629     $alpha = gsl_complex_rect(0,0);
 630     map { gsl_matrix_complex_set($C, 0,$_,$alpha) } (0..1);
 631     map { gsl_matrix_complex_set($C, 1,$_,$alpha) } (0..1);
 632
 633     $alpha = gsl_complex_rect(1,0);
 634     my $beta = gsl_complex_rect(1,0);
 635
 636     ok_status(gsl_blas_zsyr2k($CblasUpper, $CblasNoTrans, $alpha, $A, $B, $beta, $C));
 637     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 0))], [0, 0]);
 638     local $TODO = "These results follow the formula given by the documentation, don't know why it fails";
 639     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 0))], [46, 10]);
 640     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 0, 1))], [34, 15]);
 641     ok_similar([gsl_parts(gsl_matrix_complex_get($C, 1, 1))], [24, 4]);
 642 }
 643
 644 Test::Class->runtests;