#!/usr/bin/env perl6

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

subtest {
  my Math::Libgsl::Vector::Complex32 $x .= new(10);
  my Math::Libgsl::Vector::Complex32 $y .= new(10);
  $x.setall(1+2i);
  $y.setall(2+3i);
  ok cdotu($x, $y) == -40+70i, 'xy';
  ok cdotc($x, $y) == 80-10i, 'H(x)y';
  ok scnrm2($x) == 7.071067810058594, 'euclidean norm ||x||₂';
  ok scasum($x) == 30, 'absolute sum';
  ok icamax($x) == 0, 'index of the largest element';
  ok cswap($x, $y) == GSL_SUCCESS, 'swap vectors';
  is-deeply $x[^10], (2+3i xx 10), 'swap ok';
  ok ccopy($y, $x) == GSL_SUCCESS, 'copy vectors';
  is-deeply $y[^10], (1+2i xx 10), 'copy ok';
  ok caxpy(2+3i, $x, $y) == GSL_SUCCESS, 'sum y = αx + y';
  is-deeply $y[^10], (-3+9i xx 10), 'result: ok';
  cscal(2+3i, $x);
  is-deeply $x[^10], (-4+7i xx 10), 'scale a vector by a complex32';
  csscal(3, $x);
  is-deeply $x[^10], (-12+21i xx 10), 'scale a vector by a num32';
}, 'level 1 functions';

subtest {
  my $α = 1+2i; my $β = 2+3i;
  my Math::Libgsl::Matrix::Complex32 $A .= new: 2, 2;
  $A.setall(1+2i);
  my Math::Libgsl::Vector::Complex32 $x .= new: 2;
  $x.setall(1+2i);
  my Math::Libgsl::Vector::Complex32 $y .= new: 2;
  $y.setall(2+3i);
  ok cgemv(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, 1+1i);
  $x.setall(2+1i);
  ok ctrmv(CblasUpper, CblasNoTrans, CblasNonUnit, $A, $x) == GSL_SUCCESS, 'op(A)x';
  is-deeply $x[^2], (0+10i, 1+3i), 'result: ok';
  my Math::Libgsl::Matrix::Complex32 $A1 .= new: 1, 1;
  $A1.set(0, 0, .529-.348i);
  my Math::Libgsl::Vector::Complex32 $x1 .= new: 1;
  $x1.set(0, -.95+.343i);
  ok ctrsv(CblasUpper, CblasNoTrans, CblasNonUnit, $A1, $x1) == GSL_SUCCESS, 'inv op(A)x:';
  ok $x1[0] == -1.551120638847351-0.37200355529785156i, 'result: ok';
  $α = 1+0i; $β = -.3+.1i;
  $A1.set(0, 0, -.434+.837i);
  $x1.set(0, .209-.935i);
  my Math::Libgsl::Vector::Complex32 $y1 .= new: 1;
  $y1.set(0, .346-.412i);
  ok chemv(CblasUpper, $α, $A1, $x1, $β, $y1) == GSL_SUCCESS, 'sum αAx + βy';
  ok $y1[0] == -0.1533060073852539+0.5639899969100952i, 'result: ok';
  $α = 0+0i;
  $A1.set(0, 0, -.651+.856i);
  $x1.set(0, -.38-.235i);
  $y1.set(0, -.651+.856i);
  ok cgeru($α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.6510000228881836+0.8560000061988831i, 'result: ok';
  $α = 0+0i;
  $A1.set(0, 0, -.651+.856i);
  $x1.set(0, -.38-.235i);
  $y1.set(0, -.651+.856i);
  ok cgerc($α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.6510000228881836+0.8560000061988831i, 'result: ok';
  $A1.set(0, 0, -.651+.856i);
  $x1.set(0, -.38-.235i);
  ok cher(CblasUpper, .1e0, $x1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.6310375332832336+0i, 'result: ok';
  $α = -1+0i;
  $A1.set(0, 0, -.821+.954i);
  $x1.set(0, .532+.802i);
  $y1.set(0, .016-.334i);
  ok cher2(CblasUpper, $α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.3022879958152771+0i, 'result: ok';
}, 'level 2 functions';

subtest {
  my Complex $α = i; my Complex $β = 0i;
  my Math::Libgsl::Matrix::Complex32 $A1 .= new: 1, 4;
  $A1.set(0, 0, -.082-.281i);
  $A1.set(0, 1, -.096+.913i);
  $A1.set(0, 2,  .974-.706i);
  $A1.set(0, 3, -.773+.522i);
  my Math::Libgsl::Matrix::Complex32 $B1 .= new: 4, 2;
  $B1.set(0, 0,  .745-.664i);
  $B1.set(0, 1,  .352-.733i);
  $B1.set(1, 0,  .304-.555i);
  $B1.set(1, 1, -.493-.089i);
  $B1.set(2, 0,  .188+.631i);
  $B1.set(2, 1,  .235+.152i);
  $B1.set(3, 0, -.299-.731i);
  $B1.set(3, 1, -.686-.332i);
  my Math::Libgsl::Matrix::Complex32 $C1 .= new: 1, 2;
  $C1.set(0, 0, -.179-.284i);
  $C1.set(0, 1, -.996-.414i);
  ok cgemm(CblasNoTrans, CblasNoTrans, $α, $A1, $B1, $β, $C1) == GSL_SUCCESS, 'sum αop(A)op(B) + βC';
  is-deeply (gather take $C1[0; $_] for ^2), (-1.0667860507965088+1.4711639881134033i, 0.5996890068054199+0.9335318803787231i), 'result: ok';

  $α = -1+0i; $β = 1+0i;
  my Math::Libgsl::Matrix::Complex32 $A2 .= new: 1, 1;
  $A2.set(0, 0,  .476+.816i);
  my Math::Libgsl::Matrix::Complex32 $B2 .= new: 1, 2;
  $B2.set(0, 0,  .282+.852i);
  $B2.set(0, 1, -.891-.588i);
  my Math::Libgsl::Matrix::Complex32 $C2 .= new: 1, 2;
  $C2.set(0, 0,    .9+.486i);
  $C2.set(0, 1,  -.78-.637i);
  ok csymm(CblasLeft, CblasUpper, $α, $A2, $B2, $β, $C2) == GSL_SUCCESS, 'αAB + βC';
  is-deeply (gather take $C2[0; $_] for ^2), (1.4609999656677246-0.1496639847755432i, -0.8356919884681702+0.36994391679763794i), 'result: ok';

  $α = .1i; $β = .1i;
  $A2.set(0, 0, -.126+.079i);
  $B2.set(0, 0, -.954-.059i);
  $B2.set(0, 1,  .296-.988i);
  $C2.set(0, 0, -.859-.731i);
  $C2.set(0, 1,  .737+.593i);
  ok chemm(CblasLeft, CblasUpper, $α, $A2, $B2, $β, $C2) == GSL_SUCCESS, 'sum αAB + βC';
  is-deeply (gather take $C2[0; $_] for ^2), (0.07235660403966904-0.07387959957122803i, -0.07174880057573318+0.06997040659189224i), 'result: ok';

  $α = -1+0i;
  my Math::Libgsl::Matrix::Complex32 $A3 .= new: 2, 2;
  $A3.set(0, 0, -.023+.762i);
  $A3.set(0, 1, -.687-.039i);
  $A3.set(1, 0, -.459+.047i);
  $A3.set(1, 1,   .189+.33i);
  my Math::Libgsl::Matrix::Complex32 $B3 .= new: 2, 3;
  $B3.set(0, 0,  .827-.561i);
  $B3.set(0, 1,  .641-.229i);
  $B3.set(0, 2, -.884-.533i);
  $B3.set(1, 0, -.624-.138i);
  $B3.set(1, 1,  .073+.924i);
  $B3.set(1, 2, -.501-.164i);
  ok ctrmm(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.8317669630050659-0.7622190117835999i, -0.14564001560211182+0.1439259946346283i, -0.7642689943313599+0.5291420221328735i, 0.07239599525928497+0.23200200498104095i, 0.29112303256988525-0.19872601330280304i, 0.040568992495536804+0.19632600247859955i), 'result: ok';

  $α = .3+0i;
  $A3.set(0, 0,  .491-.317i);
  $A3.set(0, 1,  -.14-.739i);
  $A3.set(1, 0, -.969-.518i);
  $A3.set(1, 1,  .702-.287i);
  $B3.set(0, 0, -.962 -.38i);
  $B3.set(0, 1,  .656+.587i);
  $B3.set(0, 2, -.195-.862i);
  $B3.set(1, 0, -.679+.598i);
  $B3.set(1, 1,  .919+.714i);
  $B3.set(1, 2, -.513+.726i);
  ok ctrsm(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.2850664556026459-0.8916955590248108i, -0.46750408411026+0.5161788463592529i, 0.07422492653131485-0.8711394667625427i, -0.33813342452049255+0.11731580644845963i, 0.22961094975471497+0.3990004360675812i, -0.2965131998062134+0.18903233110904694i), 'result: ok';

  $α = 0+0i; $β = -.3+.1i;
  my Math::Libgsl::Matrix::Complex32 $A4 .= new: 2, 1;
  $A4.set(0, 0,  .447-.507i);
  $A4.set(1, 0, -.425+.701i);
  my Math::Libgsl::Matrix::Complex32 $B4 .= new: 2, 2;
  $B4.set(0, 0,  .16 -.245i);
  $B4.set(0, 1,  .922-.437i);
  $B4.set(1, 0,  .24 +.008i);
  $B4.set(1, 1, -.095+.749i);
  ok csyrk(CblasUpper, CblasNoTrans, $α, $A4, $β, $B4) == GSL_SUCCESS, 'αAT(A) + βB';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take $B4[$i; $j] } }), (-0.023499999195337296+0.08950001001358032i, -0.2328999936580658+0.2233000099658966i, 0.23999999463558197+0.00800000037997961i, -0.046400003135204315-0.23420001566410065i), 'result: ok';

  my $a = 0e0; my $b = .1e0;
  $A4.set(0, 0, -.617+.179i);
  $A4.set(1, 0, -.626+.334i);
  $B4.set(0, 0,  .346-.903i);
  $B4.set(0, 1,  .022-.839i);
  $B4.set(1, 0, -.715+.049i);
  $B4.set(1, 1, -.338+.149i);
  ok cherk(CblasUpper, CblasNoTrans, $a, $A4, $b, $B4) == GSL_SUCCESS, 'αAH(A) + βB';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take $B4[$i; $j] } }), (0.03460000082850456+0i, 0.002199999988079071-0.08389999717473984i, -0.7149999737739563+0.04899999871850014i, -0.03380000218749046+0i), 'result: ok';

  $α = -1+0i; $b = -.3e0;
  my Math::Libgsl::Matrix::Complex32 $A5 .= new: 1, 2;
  $A5.set(0, 0,  .178+.545i);
  $A5.set(0, 1, -.491+.979i);
  my Math::Libgsl::Matrix::Complex32 $B5 .= new: 1, 2;
  $B5.set(0, 0, -.665-.531i);
  $B5.set(0, 1,   -.4+.227i);
  my Math::Libgsl::Matrix::Complex32 $C5 .= new: 1, 1;
  $C5.set(0, 0, .115-.193i);
  ok cher2k(CblasUpper, CblasNoTrans, $α, $A5, $B5, $b, $C5) == GSL_SUCCESS, 'αAH(B) + αBH(A) + βC';
  ok $C5[0; 0] == -0.05623596906661987+0i, 'result: ok';
}, 'level 3 functions';

done-testing;