### ----------------------------------------------------
### -- Math::Libgsl::BLAS
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/05-complex64.rakutest
### ----------------------------------------------------
#!/usr/bin/env raku

use Test;
use lib 'lib';
use Math::Libgsl::Constants;
use Math::Libgsl::Vector::Complex64;
use Math::Libgsl::Matrix::Complex64;
use Math::Libgsl::BLAS::Complex64;

subtest {
  my Math::Libgsl::Vector::Complex64 $x .= new(10);
  my Math::Libgsl::Vector::Complex64 $y .= new(10);
  $x.setall(1+2i);
  $y.setall(2+3i);
  ok zdotu($x, $y) == -40+70i, 'xy';
  ok zdotc($x, $y) == 80-10i, 'H(x)y';
  ok dznrm2($x) == 7.0710678118654755, 'euclidean norm ||x||₂';
  ok dzasum($x) == 30, 'absolute sum';
  ok izamax($x) == 0, 'index of the largest element';
  ok zswap($x, $y) == GSL_SUCCESS, 'swap vectors';
  is-deeply $x[^10], (2+3i xx 10), 'swap ok';
  ok zcopy($y, $x) == GSL_SUCCESS, 'copy vectors';
  is-deeply $y[^10], (1+2i xx 10), 'copy ok';
  ok zaxpy(2+3i, $x, $y) == GSL_SUCCESS, 'sum y = αx + y';
  is-deeply $y[^10], (-3+9i xx 10), 'result: ok';
  zscal(2+3i, $x);
  is-deeply $x[^10], (-4+7i xx 10), 'scale a vector by a complex64';
  zdscal(3, $x);
  is-deeply $x[^10], (-12+21i xx 10), 'scale a vector by a num64';
}, 'level 1 functions';

subtest {
  my $α = 1+2i; my $β = 2+3i;
  my Math::Libgsl::Matrix::Complex64 $A .= new: 2, 2;
  $A.setall(1+2i);
  my Math::Libgsl::Vector::Complex64 $x .= new: 2;
  $x.setall(1+2i);
  my Math::Libgsl::Vector::Complex64 $y .= new: 2;
  $y.setall(2+3i);
  ok zgemv(CblasNoTrans, $α, $A, $x, $β, $y) == GSL_SUCCESS, 'αop(A)x + βy';
  is-deeply $y[^2], (-27+8i, -27+8i), 'result: ok';
  $A.set(1, 0, 2+1i);
  $A.set(1, 1, 2+3i);
  $x.setall(1+1i);
  ok ztrmv(CblasUpper, CblasNoTrans, CblasNonUnit, $A, $x) == GSL_SUCCESS, 'op(A)x';
  is-deeply $x[^2], (-2+6i, -1+5i), 'result: ok';
  my Math::Libgsl::Matrix::Complex64 $A1 .= new: 1, 1;
  $A1.set(0, 0, .977-.955i);
  my Math::Libgsl::Vector::Complex64 $x1 .= new: 1;
  $x1.set(0, -.627+.281i);
  ok ztrsv(CblasUpper, CblasNoTrans, CblasNonUnit, $A1, $x1) == GSL_SUCCESS, 'inv op(A)x:';
  ok $x1[0] == -0.47195741457252255-0.17371477064151375i, 'result: ok';
  $α = 0+0i; $β = 1+0i;
  $A1.set(0, 0, .036-.966i);
  $x1.set(0, -.695+.886i);
  my Math::Libgsl::Vector::Complex64 $y1 .= new: 1;
  $y1.set(0, .486+.629i);
  ok zhemv(CblasUpper, $α, $A1, $x1, $β, $y1) == GSL_SUCCESS, 'sum αAx + βy';
  ok $y1[0] == 0.486+0.629i, 'result: ok';
  $α = -1+0i;
  $A1.set(0, 0, -.426+.757i);
  $x1.set(0, -.579-.155i);
  $y1.set(0, .831+.035i);
  ok zgeru($α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == 0.049723999999999935+0.90607i, 'result: ok';
  $α = -1+0i;
  $A1.set(0, 0, -.426+.757i);
  $x1.set(0, -.579-.155i);
  $y1.set(0, .831+.035i);
  ok zgerc($α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == 0.06057399999999996+0.86554i, 'result: ok';
  $A1.set(0, 0, -.426+.757i);
  $x1.set(0, -.579-.155i);
  ok zher(CblasUpper, -.3e0, $x1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.5337798+0i, 'result: ok';
  $α = -.3+.1i;
  $A1.set(0, 0, -.334+.286i);
  $x1.set(0, -.14-.135i);
  $y1.set(0, .455+.358i);
  ok zher2(CblasUpper, $α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.264521+0i, 'result: ok';
}, 'level 2 functions';

subtest {
  my Complex $α = 1+0i; my Complex $β = -1+0i;
  my Math::Libgsl::Matrix::Complex64 $A1 .= new: 1, 4;
  $A1.set(0, 0,  -.33+.457i);
  $A1.set(0, 1,  .428-.19i);
  $A1.set(0, 2,   .86-.53i);
  $A1.set(0, 3,  .058-.942i);
  my Math::Libgsl::Matrix::Complex64 $B1 .= new: 4, 2;
  $B1.set(0, 0,  .434+.653i);
  $B1.set(0, 1, -.124+.191i);
  $B1.set(1, 0, -.112-.84i);
  $B1.set(1, 1,  -.72+.075i);
  $B1.set(2, 0, -.503-.109i);
  $B1.set(2, 1,    .3-.898i);
  $B1.set(3, 0,  .489+.384i);
  $B1.set(3, 1,  .993-.804i);
  my Math::Libgsl::Matrix::Complex64 $C1 .= new: 1, 2;
  $C1.set(0, 0, -.792-.155i);
  $C1.set(0, 1, -.608-.243i);
  ok zgemm(CblasNoTrans, CblasNoTrans, $α, $A1, $B1, $β, $C1) == GSL_SUCCESS, 'sum αop(A)op(B) + βC';
  is-deeply (gather take $C1[0; $_] for ^2), (0.04256299999999996-0.465908i, -0.6499910000000001-1.621116i), 'result: ok';

  $α = .1i; $β = .3i;
  my Math::Libgsl::Matrix::Complex64 $A2 .= new: 1, 1;
  $A2.set(0, 0,  .511-.486i);
  my Math::Libgsl::Matrix::Complex64 $B2 .= new: 1, 2;
  $B2.set(0, 0,  .985-.923i);
  $B2.set(0, 1, -.234-.756i);
  my Math::Libgsl::Matrix::Complex64 $C2 .= new: 1, 2;
  $C2.set(0, 0,  -.16+.049i);
  $C2.set(0, 1,  .618-.349i);
  ok zsymm(CblasLeft, CblasUpper, $α, $A2, $B2, $β, $C2) == GSL_SUCCESS, 'αAB + βC';
  is-deeply (gather take $C2[0; $_] for ^2), (0.08033630000000001-0.0425243i, 0.1319592+0.136701i), 'result: ok';

  $α = .1i; $β = .1i;
  $A2.set(0, 0, -.359+.089i);
  $B2.set(0, 0, -.451-.337i);
  $B2.set(0, 1, -.901-.871i);
  $C2.set(0, 0,  .729+.631i);
  $C2.set(0, 1,  .364+.246i);
  ok zhemm(CblasLeft, CblasUpper, $α, $A2, $B2, $β, $C2) == GSL_SUCCESS, 'sum αAB + βC';
  is-deeply (gather take $C2[0; $_] for ^2), (-0.0751983+0.0890909i, -0.0558689+0.0687459i), 'result: ok';

  $α = .3+0i;
  my Math::Libgsl::Matrix::Complex64 $A3 .= new: 2, 2;
  $A3.set(0, 0,  .463+.033i);
  $A3.set(0, 1, -.929+.949i);
  $A3.set(1, 0,  .864+.986i);
  $A3.set(1, 1,  .393+.885i);
  my Math::Libgsl::Matrix::Complex64 $B3 .= new: 2, 3;
  $B3.set(0, 0, -.321-.852i);
  $B3.set(0, 1, -.337-.175i);
  $B3.set(0, 2,  .607-.613i);
  $B3.set(1, 0,  .688+.973i);
  $B3.set(1, 1, -.331 -.35i);
  $B3.set(1, 2,  .719-.553i);
  ok ztrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $A3, $B3) == GSL_SUCCESS, 'αop(A)B';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take $B3[$i; $j] } }), (-0.5049107999999999-0.1968222i, 0.14681789999999995-0.0243345i, 0.04743480000000003+0.279684i, -0.17721630000000002+0.2973807i, 0.053900099999999986-0.1291455i, 0.2315916+0.1256958i), 'result: ok';

  $α = .3+0i;
  $A3.set(0, 0,  .189+.519i);
  $A3.set(0, 1, -.455-.444i);
  $A3.set(1, 0,  -.21-.507i);
  $A3.set(1, 1, -.591+.859i);
  $B3.set(0, 0, -.779-.484i);
  $B3.set(0, 1,  .249-.107i);
  $B3.set(0, 2, -.755-.047i);
  $B3.set(1, 0,  .941+.675i);
  $B3.set(1, 1, -.757+.645i);
  $B3.set(1, 2, -.649+.242i);
  ok ztrsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $A3, $B3) == GSL_SUCCESS, 'αop(A)B';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take $B3[$i; $j] } }), (-0.5512360375604164-0.04117743465544226i, 0.31534332655118025-0.20784391376238756i, 0.062041254407206994+0.413773152122249i, 0.006538307998256045-0.333136367901012i, 0.2763457515991177+0.074248732019699i, 0.16320575958320835+0.11437182315055161i), 'result: ok';

  $α = 1+0i; $β = 1+0i;
  my Math::Libgsl::Matrix::Complex64 $A4 .= new: 2, 1;
  $A4.set(0, 0, -.049-.687i);
  $A4.set(1, 0, -.434+.294i);
  my Math::Libgsl::Matrix::Complex64 $B4 .= new: 2, 2;
  $B4.set(0, 0,  .937-.113i);
  $B4.set(0, 1,  .796+.293i);
  $B4.set(1, 0,  .876-.199i);
  $B4.set(1, 1, -.757-.103i);
  ok zsyrk(CblasUpper, CblasNoTrans, $α, $A4, $β, $B4) == GSL_SUCCESS, 'αAT(A) + βB';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take $B4[$i; $j] } }), (0.46743199999999996-0.04567399999999999i, 1.019244+0.576752i, 0.876-0.199i, -0.65508-0.35819199999999995i), 'result: ok';

  my $a = 1e0; my $b = 1e0;
  $A4.set(0, 0, .934+.664i);
  $A4.set(1, 0, .426+.263i);
  $B4.set(0, 0,  .251-.97i);
  $B4.set(0, 1,  .76-.349i);
  $B4.set(1, 0, .152-.899i);
  $B4.set(1, 1, -.17+.707i);
  ok zherk(CblasUpper, CblasNoTrans, $a, $A4, $b, $B4) == GSL_SUCCESS, 'αAH(A) + βB';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take $B4[$i; $j] } }), (1.5642520000000002+0i, 1.332516-0.311778i, 0.152-0.899i, 0.080645+0i), 'result: ok';

  $α = 1+0i; $b = 1e0;
  my Math::Libgsl::Matrix::Complex64 $A5 .= new: 1, 2;
  $A5.set(0, 0, .515-.034i);
  $A5.set(0, 1, .067+ .66i);
  my Math::Libgsl::Matrix::Complex64 $B5 .= new: 1, 2;
  $B5.set(0, 0,  .408 -.85i);
  $B5.set(0, 1, -.945-.799i);
  my Math::Libgsl::Matrix::Complex64 $C5 .= new: 1, 1;
  $C5.set(0, 0, -.918-.985i);
  ok zher2k(CblasUpper, CblasNoTrans, $α, $A5, $B5, $b, $C5) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC';
  ok $C5[0; 0] == -1.62127+0i, 'result: ok';
}, 'level 3 functions';

done-testing;


