#!/usr/bin/env perl6

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

subtest {
  my Math::Libgsl::Vector $x .= new(10);
  my Math::Libgsl::Vector $y .= new(10);
  $x.setall(2);
  $y.setall(3);
  ok ddot($x, $y) == 60, 'scalar product of vectors';
  ok dnrm2($x) == 6.324555320336759, 'euclidean norm ||x||₂';
  ok dasum($x) == 20, 'absolute sum';
  ok idamax($x) == 0, 'index of the largest element';
  ok dswap($x, $y) == GSL_SUCCESS, 'swap vectors';
  is-deeply $x[^10], (3e0 xx 10), 'swap ok';
  ok dcopy($x, $y) == GSL_SUCCESS, 'copy vectors';
  is-deeply $y[^10], (3e0 xx 10), 'copy ok';
  $x.setall(2);
  ok daxpy(2, $x, $y) == GSL_SUCCESS, 'sum y = αx + y';
  is-deeply $y[^10], (7e0 xx 10), 'result: ok';
  dscal(2, $x);
  is-deeply $x[^10], (4e0 xx 10), 'scale a vector';
  is-deeply drotg(1, 2), (2.23606797749979, 0.44721359549995793, 0.8944271909999159)».Num, 'Givens rotation';
  is-deeply drotmg(1, 1, 3, 4), (1e0, 0.75e0, 0e0, 0e0, 0.75e0), 'modified Givens transformation';
  $x[$_] = $_ for ^10;
  drot($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 $rx .= new(1);
  my Math::Libgsl::Vector $ry .= new(1);
  $rx[0] = .84;
  $ry[0] = -.711;
  my @param = -1e0, -8.00850735044e0, 0.0204647351647e0, 1.898461360078e-04, -4.32701487194e0;
  drotm($rx, $ry, @param);
  ok $rx[0] == -6.727281154972301e0, 'modified Givens transformation: rotated vectors: x';
  ok $ry[0] == 3.093697951487688e0, 'modified Givens transformation: rotated vectors: y';
}, 'level 1 functions';

subtest {
  my $α = 3; my $β = 2;
  my Math::Libgsl::Matrix $A .= new: 2, 2;
  $A.setall(2);
  my Math::Libgsl::Vector $x .= new: 2;
  $x.setall(2);
  my Math::Libgsl::Vector $y .= new: 2;
  $y.setall(2);
  ok dgemv(CblasNoTrans, $α, $A, $x, $β, $y) == GSL_SUCCESS, 'αop(A)x + βy';
  is-deeply $y[^2], (28e0, 28e0), 'result: ok';
  $A.set(1, 1, 0);
  ok dtrmv(CblasUpper, CblasNoTrans, CblasNonUnit, $A, $x) == GSL_SUCCESS, 'op(A)x';
  is-deeply $x[^2], (8e0, 0e0), 'result: ok';
  my Math::Libgsl::Matrix $A1 .= new: 1, 1;
  $A1.set(0, 0, -.21);
  my Math::Libgsl::Vector $x1 .= new: 1;
  $x1.set(0, .473);
  ok dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, $A1, $x1) == GSL_SUCCESS, 'inv op(A)x';
  ok $x1[0] == -2.2523809523809524, 'result: ok';
  $α = 0; $β = -.3;
  $A1.set(0, 0, .544);
  $x1.set(0, -.601);
  my Math::Libgsl::Vector $y1 .= new: 1;
  $y1.set(0, -.852);
  ok dsymv(CblasUpper, $α, $A1, $x1, $β, $y1) == GSL_SUCCESS, 'sum αAx + βy';
  ok $y1[0] == 0.2556e0, 'result: ok';
  $α = 1;
  $A1.set(0, 0, -.809);
  $x1.set(0, -.652);
  $y1.set(0, .712);
  ok dger($α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(y) + A';
  ok $A1[0;0] == -1.273224e0, 'result: ok';
  $α = -.3;
  $A1.set(0, 0, -.65);
  $x1.set(0, -.891);
  ok dsyr(CblasUpper, $α, $x1, $A1) == GSL_SUCCESS, 'αxT(x) + A';
  ok $A1[0;0] == -0.8881643e0, 'result: ok';
  $α = 0;
  $A1.set(0, 0, -.824);
  $x1.set(0, .684);
  $y1.set(0, .965);
  ok dsyr2(CblasUpper, $α, $x1, $y1, $A1) == GSL_SUCCESS, 'αxT(x) + αyT(x) + A';
  ok $A1[0;0] == -0.824e0, 'result: ok';
}, 'level 2 functions';

subtest {
  my $α = 0; my $β = 0;
  my Math::Libgsl::Matrix $A1 .= new: 1, 4;
  $A1.set(0, 0,  .017);
  $A1.set(0, 1,  .191);
  $A1.set(0, 2,  .863);
  $A1.set(0, 3,  -.97);
  my Math::Libgsl::Matrix $B1 .= new: 4, 2;
  $B1.set(0, 0, -.207);
  $B1.set(0, 1, -.916);
  $B1.set(1, 0, -.278);
  $B1.set(1, 1,  .403);
  $B1.set(2, 0,  .885);
  $B1.set(2, 1,  .409);
  $B1.set(3, 0, -.772);
  $B1.set(3, 1,  -.27);
  my Math::Libgsl::Matrix $C1 .= new: 1, 2;
  $C1.set(0, 0, -.274);
  $C1.set(0, 0, -.858);
  ok dgemm(CblasNoTrans, CblasNoTrans, $α, $A1, $B1, $β, $C1) == GSL_SUCCESS, 'sum αop(A)op(B) + βC';
  is-deeply (gather take $C1[0; $_] for ^2), (0e0, 0e0), 'result: ok';

  $α = -.3; $β = 1;
  my Math::Libgsl::Matrix $A2 .= new: 1, 1;
  $A2.set(0, 0, -.981);
  my Math::Libgsl::Matrix $B2 .= new: 1, 2;
  $B2.set(0, 0, -.823);
  $B2.set(0, 1,   .83);
  my Math::Libgsl::Matrix $C2 .= new: 1, 2;
  $C2.set(0, 0,  .991);
  $C2.set(0, 1,  .382);
  ok dsymm(CblasLeft, CblasUpper, $α, $A2, $B2, $β, $C2) == GSL_SUCCESS, 'αAB + βC';
  is-deeply (gather take $C2[0; $_] for ^2), (0.7487911e0, 0.626269e0), 'result: ok';

  $α = -.3;
  my Math::Libgsl::Matrix $A3 .= new: 2, 2;
  $A3.set(0, 0,  .174);
  $A3.set(0, 1, -.308);
  $A3.set(1, 0,  .997);
  $A3.set(1, 1, -.484);
  my Math::Libgsl::Matrix $B3 .= new: 2, 3;
  $B3.set(0, 0, -.256);
  $B3.set(0, 1, -.178);
  $B3.set(0, 2,  .098);
  $B3.set(1, 0,  .004);
  $B3.set(1, 1,   .97);
  $B3.set(1, 2, -.408);
  ok dtrmm(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.013732799999999998e0, 0.09891959999999998e0, -0.0428148e0, 0.0005808e0, 0.14084399999999997e0, -0.05924159999999999e0), 'result: ok';

  $α = .1; $β = -.3;
  my Math::Libgsl::Matrix $A4 .= new: 2, 1;
  $A4.set(0, 0,  .169);
  $A4.set(1, 0, -.875);
  my Math::Libgsl::Matrix $B4 .= new: 2, 2;
  $B4.set(0, 0, .159);
  $B4.set(0, 1, .277);
  $B4.set(1, 0, .865);
  $B4.set(1, 1, .346);
  ok dsyrk(CblasUpper, CblasNoTrans, $α, $A4, $β, $B4) == GSL_SUCCESS, 'αAT(A) + βB';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take $B4[$i; $j] } }), (-0.0448439e0, -0.09788750000000002e0, 0.865e0, -0.027237499999999984e0), 'result: ok';

  $α = .1; $β = 0;
  my Math::Libgsl::Matrix $A5 .= new: 1, 2;
  $A5.set(0, 0, -.225);
  $A5.set(0, 1,  .857);
  my Math::Libgsl::Matrix $B5 .= new: 1, 2;
  $B5.set(0, 0, -.933);
  $B5.set(0, 1,  .994);
  my Math::Libgsl::Matrix $C5 .= new: 1, 1;
  $C5.set(0, 0, .177);
  ok dsyr2k(CblasUpper, CblasNoTrans, $α, $A5, $B5, $β, $C5) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC';
  ok $C5[0; 0] == 0.21235660000000003e0, 'result: ok';
}, 'level 3 functions';

done-testing;