#!/usr/bin/env perl6

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

subtest {
  my Math::Libgsl::Vector::Num32 $x .= new(10);
  my Math::Libgsl::Vector::Num32 $y .= new(10);
  $x.setall(2);
  $y.setall(3);
  ok sdsdot(5, $x, $y) == 65, 'α + x · y';

  ok sdot($x, $y) == 60, 'scalar product of vectors';
  ok dsdot($x, $y) == 60, 'scalar product of vectors';
  ok snrm2($x) == 6.324555397033691, 'euclidean norm ||x||₂';
  ok sasum($x) == 20, 'absolute sum';
  ok isamax($x) == 0, 'index of the largest element';
  ok sswap($x, $y) == GSL_SUCCESS, 'swap vectors';
  is-deeply $x[^10], (3e0 xx 10), 'swap ok';
  ok scopy($x, $y) == GSL_SUCCESS, 'copy vectors';
  is-deeply $y[^10], (3e0 xx 10), 'copy ok';
  $x.setall(2);
  ok saxpy(2, $x, $y) == GSL_SUCCESS, 'sum y = αx + y';
  is-deeply $y[^10], (7e0 xx 10), 'result: ok';
  sscal(2, $x);
  is-deeply $x[^10], (4e0 xx 10), 'scale a vector';
  is-deeply srotg(1, 2), (2.2360680103302, 0.4472135901451111, 0.8944271802902222)».Num, 'Givens rotation';
  is-deeply srotmg(1, 1, 3, 4), (1e0, 0.75e0, 0e0, 0e0, 0.75e0), 'modified Givens transformation';
  $x[$_] = $_ for ^10;
  srot($x, $y, 1, 2);
  is-deeply $x[^10], (14e0, 15e0, 16e0, 17e0, 18e0, 19e0, 20e0, 21e0, 22e0, 23e0), 'Givens rotation: x';
  is-deeply $y[^10], (7e0, 5e0, 3e0, 1e0, -1e0, -3e0, -5e0, -7e0, -9e0, -11e0), 'Givens rotation: y';
  my Math::Libgsl::Vector::Num32 $rx .= new(1);
  my Math::Libgsl::Vector::Num32 $ry .= new(1);
  $rx[0] = .84;
  $ry[0] = -.711;
  my @param = -1, -8.00850735044, 0.0204647351647, 1.898461360078e-04, -4.32701487194;
  srotm($rx, $ry, @param);
  ok $rx[0] == -6.727281093597412e0, 'modified Givens transformation: rotated vectors: x';
  ok $ry[0] == 3.093698263168335e0, 'modified Givens transformation: rotated vectors: y';
}, 'level 1 functions';

subtest {
  my $α = 3; my $β = 2;
  my Math::Libgsl::Matrix::Num32 $A .= new: 2, 2;
  $A.setall(2);
  my Math::Libgsl::Vector::Num32 $x .= new: 2;
  $x.setall(2);
  my Math::Libgsl::Vector::Num32 $y .= new: 2;
  $y.setall(2);
  ok sgemv(CblasNoTrans, $α, $A, $x, $β, $y) == GSL_SUCCESS, 'αop(A)x + βy';
  is-deeply $y[^2], (28e0, 28e0), 'result: ok';
  $A.set(1, 1, 0);
  ok strmv(CblasUpper, CblasNoTrans, CblasNonUnit, $A, $x) == GSL_SUCCESS, 'op(A)x';
  is-deeply $x[^2], (8e0, 0e0), 'result: ok';
  my Math::Libgsl::Matrix::Num32 $A1 .= new: 1, 1;
  $A1.set(0, 0, -.21);
  my Math::Libgsl::Vector::Num32 $x1 .= new: 1;
  $x1.set(0, .473);
  ok strsv(CblasUpper, CblasNoTrans, CblasNonUnit, $A1, $x1) == GSL_SUCCESS, 'inv op(A)x';
  ok $x1[0] == -2.2523810863494873, 'result: ok';
  $α = 1; $β = -1;
  $A1.set(0, 0, -.428);
  $x1.set(0, -.34);
  my Math::Libgsl::Vector::Num32 $y1 .= new: 1;
  $y1.set(0, -.888);
  ok ssymv(CblasUpper, $α, $A1, $x1, $β, $y1) == GSL_SUCCESS, 'sum αAx + βy';
  ok $y1[0] == 1.033519983291626e0, 'result: ok';
  $α = 1;
  $A1.set(0, 0, -.515);
  $x1.set(0, .611);
  $y1.set(0, -.082);
  ok sger($α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -0.5651019811630249e0, 'result: ok';
  $α = .1;
  $A1.set(0, 0, -.291);
  $x1.set(0, .845);
  ok ssyr(CblasUpper, $α, $x1, $A1) == GSL_SUCCESS, 'αxT(x) + A';
  ok $A1[0;0] == -0.21959750354290009e0, 'result: ok';
  $α = 0;
  $A1.set(0, 0, .862);
  $x1.set(0, .823);
  $y1.set(0, .699);
  ok ssyr2(CblasUpper, $α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(x) + αyT(x) + A';
  ok $A1[0;0] == 0.8619999885559082e0, 'result: ok';
}, 'level 2 functions';

subtest {
  my $α = 1; my $β = 0;
  my Math::Libgsl::Matrix::Num32 $A1 .= new: 1, 4;
  $A1.set(0, 0,  .199);
  $A1.set(0, 1,  .237);
  $A1.set(0, 2,  .456);
  $A1.set(0, 3,  .377);
  my Math::Libgsl::Matrix::Num32 $B1 .= new: 4, 2;
  $B1.set(0, 0,  .842);
  $B1.set(0, 1, -.734);
  $B1.set(1, 0,  .323);
  $B1.set(1, 1, -.957);
  $B1.set(2, 0, -.303);
  $B1.set(2, 1, -.873);
  $B1.set(3, 0, -.871);
  $B1.set(3, 1, -.819);
  my Math::Libgsl::Matrix::Num32 $C1 .= new: 1, 2;
  $C1.set(0, 0,  .498);
  $C1.set(0, 0, -.925);
  ok sgemm(CblasNoTrans, CblasNoTrans, $α, $A1, $B1, $β, $C1) == GSL_SUCCESS, 'sum αop(A)op(B) + βC';
  is-deeply (gather take $C1[0; $_] for ^2), (-0.22242599725723267e0, -1.0797260999679565e0), 'result: ok';

  $α = -.3; $β = -1;
  my Math::Libgsl::Matrix::Num32 $A2 .= new: 1, 1;
  $A2.set(0, 0, -.581);
  my Math::Libgsl::Matrix::Num32 $B2 .= new: 1, 2;
  $B2.set(0, 0,  .157);
  $B2.set(0, 1,  .451);
  my Math::Libgsl::Matrix::Num32 $C2 .= new: 1, 2;
  $C2.set(0, 0, -.869);
  $C2.set(0, 1, -.871);
  ok ssymm(CblasLeft, CblasUpper, $α, $A2, $B2, $β, $C2) == GSL_SUCCESS, 'αAB + βC';
  is-deeply (gather take $C2[0; $_] for ^2), (0.8963651061058044e0, 0.9496092796325684e0), 'result: ok';

  $α = -.3;
  my Math::Libgsl::Matrix::Num32 $A3 .= new: 2, 2;
  $A3.set(0, 0,   .18);
  $A3.set(0, 1,  .199);
  $A3.set(1, 0,  .122);
  $A3.set(1, 1, -.547);
  my Math::Libgsl::Matrix::Num32 $B3 .= new: 2, 3;
  $B3.set(0, 0, -.874);
  $B3.set(0, 1, -.383);
  $B3.set(0, 2,  .458);
  $B3.set(1, 0,  .124);
  $B3.set(1, 1, -.221);
  $B3.set(1, 2, -.107);
  ok strmm(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.03979320451617241e0, 0.03387570381164551e0, -0.018344102427363396e0, 0.020348399877548218e0, -0.03626609966158867e0, -0.01755870133638382e0), 'result: ok';

  $α = -.3;
  $A3.set(0, 0, -.279);
  $A3.set(0, 1,  .058);
  $A3.set(1, 0,  .437);
  $A3.set(1, 1,  .462);
  $B3.set(0, 0,  .578);
  $B3.set(0, 1,  .473);
  $B3.set(0, 2,  -.34);
  $B3.set(1, 0, -.128);
  $B3.set(1, 1,  .503);
  $B3.set(1, 2,    .2);
  ok strsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $A3, $B3) == GSL_SUCCESS, 'αop(inv(A))B';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take $B3[$i; $j] } }), (0.6387841701507568e0, 0.4407019317150116e0, -0.3925895094871521e0, 0.08311688154935837e0, -0.3266233801841736e0, -0.12987013161182404e0), 'result: ok';

  $α = -1e0; $β = .1;
  my Math::Libgsl::Matrix::Num32 $A4 .= new: 2, 1;
  $A4.set(0, 0,  .412);
  $A4.set(1, 0, -.229);
  my Math::Libgsl::Matrix::Num32 $B4 .= new: 2, 2;
  $B4.set(0, 0,  .628);
  $B4.set(0, 1, -.664);
  $B4.set(1, 0, -.268);
  $B4.set(1, 1,  .096);
  ok ssyrk(CblasUpper, CblasNoTrans, $α, $A4, $β, $B4) == GSL_SUCCESS, 'αAT(A) + βB';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take $B4[$i; $j] } }), (-0.1069439947605133e0, 0.02794799953699112e0, -0.2680000066757202e0, -0.042841002345085144e0), 'result: ok';

  $α = .1; $β = 1;
  my Math::Libgsl::Matrix::Num32 $A5 .= new: 1, 2;
  $A5.set(0, 0, -.915);
  $A5.set(0, 1,  .445);
  my Math::Libgsl::Matrix::Num32 $B5 .= new: 1, 2;
  $B5.set(0, 0,  .213);
  $B5.set(0, 1, -.194);
  my Math::Libgsl::Matrix::Num32 $C5 .= new: 1, 1;
  $C5.set(0, 0, -.117);
  ok ssyr2k(CblasUpper, CblasNoTrans, $α, $A5, $B5, $β, $C5) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC';
  ok $C5[0; 0] == -0.17324499785900116e0, 'result: ok';
}, 'level 3 functions';

done-testing;