### ----------------------------------------------------
### -- Math::Libgsl::BLAS
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw.rakutest
### ----------------------------------------------------
### -- Chunk 1 of 4
#!/usr/bin/env raku

use Test;
use lib 'lib';
use Math::Libgsl::Constants;
use Math::Libgsl::Raw::Complex;
use Math::Libgsl::Raw::Matrix :ALL;
use Math::Libgsl::Raw::Matrix::Num32 :ALL;
use Math::Libgsl::Raw::Matrix::Complex64 :ALL;
use Math::Libgsl::Raw::Matrix::Complex32 :ALL;
use Math::Libgsl::Raw::BLAS :ALL;
use NativeCall;

subtest {
  my gsl_vector_float $x = gsl_vector_float_calloc(10);
  gsl_vector_float_set_all($x, 2e0);
  my gsl_vector_float $y = gsl_vector_float_calloc(10);
  gsl_vector_float_set_all($y, 3e0);
  my num32 $alpha = 5e0;
  my num32 $result;
  ok gsl_blas_sdsdot($alpha, $x, $y, $result) == GSL_SUCCESS, 'gsl_blas_sdsdot returns success';
  ok $result == 65, 'α + x · y (vectors: num32, result: num32)';

  ok gsl_blas_sdot($x, $y, $result) == GSL_SUCCESS, 'gsl_blas_sdot returns success';
  ok $result == 60, 'x · y (vectors: num32, result: num32)';

  my num64 $dresult;
  ok gsl_blas_dsdot($x, $y, $dresult) == GSL_SUCCESS, 'gsl_blas_dsdot returns success';
  ok $dresult == 60, 'x · y (vectors: num32, result: num64)';

  my gsl_vector $dx = gsl_vector_calloc(10);
  gsl_vector_set_all($dx, 2e0);
  my gsl_vector $dy = gsl_vector_calloc(10);
  gsl_vector_set_all($dy, 3e0);
  ok gsl_blas_ddot($dx, $dy, $dresult) == GSL_SUCCESS, 'gsl_blas_ddot returns success';
  ok $dresult == 60, 'x · y (vectors: num64, result: num64)';

  my gsl_vector_complex_float $cx =  gsl_vector_complex_float_calloc(10);
  my gsl_complex_float $cf = alloc_gsl_complex_float;
  mgsl_complex_float_rect(1e0, 2e0, $cf);
  mgsl_vector_complex_float_set_all($cx, $cf);
  my gsl_vector_complex_float $cy =  gsl_vector_complex_float_calloc(10);
  mgsl_complex_float_rect(2e0, 3e0, $cf);
  mgsl_vector_complex_float_set_all($cy, $cf);
  my gsl_complex_float $cresult = alloc_gsl_complex_float;
  ok gsl_blas_cdotu($cx, $cy, $cresult) == GSL_SUCCESS, 'gsl_blas_cdotu returns success';
  ok $cresult.dat[0] == -40, 'result: real part';
  ok $cresult.dat[1] ==  70, 'result: imaginary part';

  my gsl_vector_complex $zx =  gsl_vector_complex_calloc(10);
  my gsl_complex $zf = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 2e0, $zf);
  mgsl_vector_complex_set_all($zx, $zf);
  my gsl_vector_complex $zy =  gsl_vector_complex_calloc(10);
  mgsl_complex_rect(2e0, 3e0, $zf);
  mgsl_vector_complex_set_all($zy, $zf);
  my gsl_complex $zresult = alloc_gsl_complex;
  ok gsl_blas_zdotu($zx, $zy, $zresult) == GSL_SUCCESS, 'gsl_blas_zdotu returns success';
  ok $zresult.dat[0] == -40, 'result: real part';
  ok $zresult.dat[1] ==  70, 'result: imaginary part';

  ok gsl_blas_cdotc($cx, $cy, $cresult) == GSL_SUCCESS, 'gsl_blas_cdotc returns success';
  ok $cresult.dat[0] ==  80, 'result: real part';
  ok $cresult.dat[1] == -10, 'result: imaginary part';

  ok gsl_blas_zdotc($zx, $zy, $zresult) == GSL_SUCCESS, 'gsl_blas_zdotc returns success';
  ok $zresult.dat[0] ==  80, 'result: real part';
  ok $zresult.dat[1] == -10, 'result: imaginary part';

  ok gsl_blas_snrm2($x) == 6.324555397033691, 'gsl_blas_snrm2: num32 euclidean norm ||x||₂';
  ok gsl_blas_dnrm2($dx) == 6.324555320336759, 'gsl_blas_dnrm2: num64 euclidean norm ||x||₂';
  ok gsl_blas_scnrm2($cx) == 7.071067810058594, 'gsl_blas_scnrm2: complex32 euclidean norm ||x||₂';
  ok gsl_blas_dznrm2($zx) == 7.0710678118654755, 'gsl_blas_dznrm2: complex64 euclidean norm ||x||₂';
  ok gsl_blas_sasum($x) == 20, 'gsl_blas_sasum: num32 absolute sum';
  ok gsl_blas_dasum($dx) == 20, 'gsl_blas_dasum: num64 absolute sum';
  ok gsl_blas_scasum($cx) == 30, 'gsl_blas_scasum: complex32 absolute sum';
  ok gsl_blas_dzasum($zx) == 30, 'gsl_blas_dzasum: complex64 absolute sum';
  ok gsl_blas_isamax($x) == 0, 'index of the largest element: num32';
  ok gsl_blas_idamax($dx) == 0, 'index of the largest element: num64';
  ok gsl_blas_icamax($cx) == 0, 'index of the largest element: complex32';
  ok gsl_blas_izamax($zx) == 0, 'index of the largest element: complex64';
  ok gsl_blas_sswap($x, $y) == GSL_SUCCESS, 'swap vectors: num32';
  is-deeply (gather take gsl_vector_float_get($x, $_) for ^10), (3e0 xx 10), 'swapped vectors';
  ok gsl_blas_dswap($dx, $dy) == GSL_SUCCESS, 'swap vectors: num64';
  is-deeply (gather take gsl_vector_get($dx, $_) for ^10), (3e0 xx 10), 'swapped vectors';
  ok gsl_blas_cswap($cx, $cy) == GSL_SUCCESS, 'swap vectors: complex32';
  my gsl_complex_float $cres = alloc_gsl_complex_float;
  is-deeply (gather for ^10 { mgsl_vector_complex_float_get($cx, $_, $cres); take $cres.dat[0], $cres.dat[1] }),
            ((2e0, 3e0) xx 10), 'swapped vectors';
  ok gsl_blas_zswap($zx, $zy) == GSL_SUCCESS, 'swap vectors: complex64';
  my gsl_complex $zres = alloc_gsl_complex;
  is-deeply (gather for ^10 { mgsl_vector_complex_get($zx, $_, $zres); take $zres.dat[0], $zres.dat[1] }),
            ((2e0, 3e0) xx 10), 'swapped vectors';
  ok gsl_blas_scopy($y, $x) == GSL_SUCCESS, 'copy vectors: num32';
  is-deeply (gather take gsl_vector_float_get($x, $_) for ^10), (2e0 xx 10), 'copyped vectors';
  ok gsl_blas_dcopy($dy, $dx) == GSL_SUCCESS, 'copy vectors: num64';
  is-deeply (gather take gsl_vector_get($dx, $_) for ^10), (2e0 xx 10), 'copyped vectors';
  ok gsl_blas_ccopy($cy, $cx) == GSL_SUCCESS, 'copy vectors: complex32';
  is-deeply (gather for ^10 { mgsl_vector_complex_float_get($cx, $_, $cres); take $cres.dat[0], $cres.dat[1] }),
            ((1e0, 2e0) xx 10), 'copied vectors';
  ok gsl_blas_zcopy($zy, $zx) == GSL_SUCCESS, 'copy vectors: complex64';
  is-deeply (gather for ^10 { mgsl_vector_complex_get($zx, $_, $zres); take $zres.dat[0], $zres.dat[1] }),
            ((1e0, 2e0) xx 10), 'copied vectors';

  gsl_vector_float_set_all($y, 3e0);
  ok gsl_blas_saxpy(2e0, $x, $y) == GSL_SUCCESS, 'sum y = αx + y: num32';
  is-deeply (gather take gsl_vector_float_get($y, $_) for ^10), (7e0 xx 10), 'result: ok';

  gsl_vector_set_all($dy, 3e0);
  ok gsl_blas_daxpy(2e0, $dx, $dy) == GSL_SUCCESS, 'sum y = αx + y: num64';
  is-deeply (gather take gsl_vector_get($dy, $_) for ^10), (7e0 xx 10), 'result: ok';

  ok mgsl_blas_caxpy($cf, $cx, $cy) == GSL_SUCCESS, 'sum y = αx + y: complex32';
  is-deeply (gather for ^10 { mgsl_vector_complex_float_get($cy, $_, $cres); take $cres.dat[0], $cres.dat[1] }),
            ((-3e0, 9e0) xx 10), 'result: ok';

  ok mgsl_blas_zaxpy($zf, $zx, $zy) == GSL_SUCCESS, 'sum y = αx + y: complex64';
  is-deeply (gather for ^10 { mgsl_vector_complex_get($zy, $_, $zres); take $zres.dat[0], $zres.dat[1] }),
            ((-3e0, 9e0) xx 10), 'result: ok';

  gsl_blas_sscal(3e0, $x);
  is-deeply (gather take gsl_vector_float_get($x, $_) for ^10), (6e0 xx 10), 'rescale: num32';

  gsl_blas_dscal(3e0, $dx);
  is-deeply (gather take gsl_vector_get($dx, $_) for ^10), (6e0 xx 10), 'rescale: num64';

  mgsl_vector_complex_float_set_all($cy, $cf);
  mgsl_blas_cscal($cf, $cx);
  is-deeply (gather for ^10 { mgsl_vector_complex_float_get($cx, $_, $cres); take $cres.dat[0], $cres.dat[1] }),
            ((-4e0, 7e0) xx 10), 'rescaled vector: complex32 * complex32 vector';

  mgsl_vector_complex_set_all($zy, $zf);
  mgsl_blas_zscal($zf, $zx);
  is-deeply (gather for ^10 { mgsl_vector_complex_get($zx, $_, $zres); take $zres.dat[0], $zres.dat[1] }),
            ((-4e0, 7e0) xx 10), 'rescaled vector: complex64 * complex64 vector';

  gsl_blas_csscal(3e0, $cx);
  is-deeply (gather for ^10 { mgsl_vector_complex_float_get($cx, $_, $cres); take $cres.dat[0], $cres.dat[1] }),
            ((-12e0, 21e0) xx 10), 'rescaled vector: num32 * complex32';

  gsl_blas_zdscal(3e0, $zx);
  is-deeply (gather for ^10 { mgsl_vector_complex_get($zx, $_, $zres); take $zres.dat[0], $zres.dat[1] }),
            ((-12e0, 21e0) xx 10), 'rescaled vector: num64 * complex64';

  my num32 $a = 1e0; my num32 $b = 2e0; my num32 $c; my num32 $s;
  ok gsl_blas_srotg($a, $b, $c, $s) == GSL_SUCCESS, 'Givens rotation: num32';
  ok $a == 2.2360680103302, 'R: ok';
  ok $c == 0.4472135901451111, 'C: ok';
  ok $s == 0.8944271802902222, 'S: ok';

  my num64 $da = 1e0; my num64 $db = 2e0; my num64 $dc; my num64 $ds;
  ok gsl_blas_drotg($da, $db, $dc, $ds) == GSL_SUCCESS, 'Givens rotation: num64';
  ok $da == 2.23606797749979, 'R: ok';
  ok $dc == 0.44721359549995793, 'C: ok';
  ok $ds == 0.8944271909999159, 'S: ok';

  gsl_vector_float_set($x, $_, $_.Num) for ^10;
  $c = 1e0; $s = 2e0;
  ok gsl_blas_srot($x, $y, $c, $s) == GSL_SUCCESS, 'Givens rotation: vectors of num32-s';
  is-deeply (gather take gsl_vector_float_get($x, $_) for ^10), (14e0, 15e0, 16e0, 17e0, 18e0, 19e0, 20e0, 21e0, 22e0, 23e0), "x': ok";
  is-deeply (gather take gsl_vector_float_get($y, $_) for ^10), (7e0, 5e0, 3e0, 1e0, -1e0, -3e0, -5e0, -7e0, -9e0, -11e0), "y': ok";

  gsl_vector_set($dx, $_, $_.Num) for ^10;
  $dc = 1e0; $ds = 2e0;
  ok gsl_blas_drot($dx, $dy, $dc, $ds) == GSL_SUCCESS, 'Givens rotation: vectors of num64-s';
  is-deeply (gather take gsl_vector_get($dx, $_) for ^10), (14e0, 15e0, 16e0, 17e0, 18e0, 19e0, 20e0, 21e0, 22e0, 23e0), "x': ok";
  is-deeply (gather take gsl_vector_get($dy, $_) for ^10), (7e0, 5e0, 3e0, 1e0, -1e0, -3e0, -5e0, -7e0, -9e0, -11e0), "y': ok";

  my num32 $d1 = 1e0; my num32 $d2 = 1e0; my num32 $b1 = 3e0; my num32 $b2 = 4e0;
  my $P = CArray[num32].allocate(5);
  ok gsl_blas_srotmg($d1, $d2, $b1, $b2, $P) == GSL_SUCCESS, 'Modified Givens transformation: num32';
  ok $b1 == 6.25e0, 'rotated vector before scaling';
  is-deeply $P.list, (1e0, 0.75e0, 0e0, 0e0, 0.75e0), 'Givens transformation matrix';

  my num64 $dd1 = 1e0; my num64 $dd2 = 1e0; my num64 $db1 = 3e0; my num64 $db2 = 4e0;
  my $dP = CArray[num64].allocate(5);
  ok gsl_blas_drotmg($dd1, $dd2, $db1, $db2, $dP) == GSL_SUCCESS, 'Modified Givens transformation: num64';
  ok $db1 == 6.25e0, 'rotated vector before scaling';
  is-deeply $dP.list, (1e0, 0.75e0, 0e0, 0e0, 0.75e0), 'Givens transformation matrix';

  my gsl_vector_float $rx = gsl_vector_float_calloc(1);
  gsl_vector_float_set($rx, 0, -.034e0);
  my gsl_vector_float $ry = gsl_vector_float_calloc(1);
  gsl_vector_float_set($ry, 0, -.56e0);
  $P[0] = -1e0;
  $P[1] = -4.44982e3;
  $P[2] = -15.5826e0;
  $P[3] = 7.091334e4;
  $P[4] = 2.95912e4;
  ok gsl_blas_srotm($rx, $ry, $P) == GSL_SUCCESS, 'Modified Givens transformation: num32 vectors';
  ok gsl_vector_float_get($rx, 0) == -39560.1796875e0, 'rotated vectors: x';
  ok gsl_vector_float_get($ry, 0) == -16570.54296875e0, 'rotated vectors: y';

  my gsl_vector $rdx = gsl_vector_calloc(1);
  gsl_vector_set_all($rdx, .84e0);
  my gsl_vector $rdy = gsl_vector_calloc(1);
  gsl_vector_set_all($rdy, -.711e0);
  $dP[0] = -1e0;
  $dP[1] = -8.00850735044e0;
  $dP[2] = 0.0204647351647e0;
  $dP[3] = 1.898461360078e-04;
  $dP[4] = -4.32701487194e0;
  ok gsl_blas_drotm($rdx, $rdy, $dP) == GSL_SUCCESS, 'Modified Givens transformation: num64 vectors';
  ok gsl_vector_get($rdx, 0) == -6.727281154972301e0, 'rotated vectors: x';
  ok gsl_vector_get($rdy, 0) == 3.093697951487688e0, 'rotated vectors: y';
}, 'level 1 functions';

subtest {
  my num32 $α = 3e0; my num32 $β = 2e0;
  my gsl_matrix_float $A = gsl_matrix_float_calloc(2, 2);
  gsl_matrix_float_set_all($A, 2e0);
  my gsl_vector_float $x = gsl_vector_float_calloc(2);
  gsl_vector_float_set_all($x, 2e0);
  my gsl_vector_float $y = gsl_vector_float_calloc(2);
  gsl_vector_float_set_all($y, 2e0);
  ok gsl_blas_sgemv(CblasNoTrans, $α, $A, $x, $β, $y) == GSL_SUCCESS, 'αop(A)x + βy: num32';
  is-deeply (gather take gsl_vector_float_get($y, $_) for ^2), (28e0, 28e0), 'result: ok';

  my num64 $dα = 3e0; my num64 $dβ = 2e0;
  my gsl_matrix $dA = gsl_matrix_calloc(2, 2);
  gsl_matrix_set_all($dA, 2e0);
  my gsl_vector $dx = gsl_vector_calloc(2);
  gsl_vector_set_all($dx, 2e0);
  my gsl_vector $dy = gsl_vector_calloc(2);
  gsl_vector_set_all($dy, 2e0);
  ok gsl_blas_dgemv(CblasNoTrans, $dα, $dA, $dx, $dβ, $dy) == GSL_SUCCESS, 'αop(A)x + βy: num64';
  is-deeply (gather take gsl_vector_get($dy, $_) for ^2), (28e0, 28e0), 'result: ok';

  my gsl_vector_complex_float $cx = gsl_vector_complex_float_calloc(2);
  my gsl_complex_float $cα = alloc_gsl_complex_float;
  mgsl_complex_float_rect(1e0, 2e0, $cα);
  mgsl_vector_complex_float_set_all($cx, $cα);
  my gsl_vector_complex_float $cy = gsl_vector_complex_float_calloc(2);
  my gsl_complex_float $cβ = alloc_gsl_complex_float;
  mgsl_complex_float_rect(2e0, 3e0, $cβ);
  mgsl_vector_complex_float_set_all($cy, $cβ);
  my gsl_matrix_complex_float $cA = gsl_matrix_complex_float_calloc(2, 2);
  mgsl_matrix_complex_float_set_all($cA, $cα);
  ok mgsl_blas_cgemv(CblasNoTrans, $cα, $cA, $cx, $cβ, $cy) == GSL_SUCCESS, 'αop(A)x + βy: complex32';
  my gsl_complex_float $cres = alloc_gsl_complex_float;
  is-deeply (gather for ^2 { mgsl_vector_complex_float_get($cy, $_, $cres); take $cres.dat[0], $cres.dat[1] }), ((-27e0, 8e0), (-27e0, 8e0)), 'result: ok';

  my gsl_vector_complex $zx = gsl_vector_complex_calloc(2);
  my gsl_complex $zα = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 2e0, $zα);
  mgsl_vector_complex_set_all($zx, $zα);
  my gsl_vector_complex $zy = gsl_vector_complex_calloc(2);
  my gsl_complex $zβ = alloc_gsl_complex;
  mgsl_complex_rect(2e0, 3e0, $zβ);
  mgsl_vector_complex_set_all($zy, $zβ);
  my gsl_matrix_complex $zA = gsl_matrix_complex_calloc(2, 2);
  mgsl_matrix_complex_set_all($zA, $zα);
  ok mgsl_blas_zgemv(CblasNoTrans, $zα, $zA, $zx, $zβ, $zy) == GSL_SUCCESS, 'αop(A)x + βy: complex64';
  my gsl_complex $zres = alloc_gsl_complex;
  is-deeply (gather for ^2 { mgsl_vector_complex_get($zy, $_, $zres); take $zres.dat[0], $zres.dat[1] }), ((-27e0, 8e0), (-27e0, 8e0)), 'result: ok';

  my gsl_matrix_float $A2 = gsl_matrix_float_calloc(2, 2);
  gsl_matrix_float_set_all($A2, 2e0);
  gsl_matrix_float_set($A2, 1, 1, 0e0);
  my gsl_vector_float $x2 = gsl_vector_float_calloc(2);
  gsl_vector_float_set_all($x2, 2e0);
  ok gsl_blas_strmv(CblasUpper, CblasNoTrans, CblasNonUnit, $A2, $x2) == GSL_SUCCESS, 'op(A)x: num32';
  is-deeply (gather take gsl_vector_float_get($x2, $_) for ^2), (8e0, 0e0), 'result: ok';

  my gsl_matrix $dA2 = gsl_matrix_calloc(2, 2);
  gsl_matrix_set_all($dA2, 2e0);
  gsl_matrix_set($dA2, 1, 1, 0e0);
  my gsl_vector $dx2 = gsl_vector_calloc(2);
  gsl_vector_set_all($dx2, 2e0);
  ok gsl_blas_dtrmv(CblasUpper, CblasNoTrans, CblasNonUnit, $dA2, $dx2) == GSL_SUCCESS, 'op(A)x: num64';
  is-deeply (gather take gsl_vector_get($dx2, $_) for ^2), (8e0, 0e0), 'result: ok';

  my gsl_matrix_complex_float $cA2 = gsl_matrix_complex_float_calloc(2, 2);
  my gsl_complex_float $c2 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(1e0, 2e0, $c2);
  mgsl_matrix_complex_float_set($cA2, 0, 0, $c2);
  mgsl_matrix_complex_float_set($cA2, 0, 1, $c2);
  mgsl_complex_float_rect(2e0, 1e0, $c2);
  mgsl_matrix_complex_float_set($cA2, 1, 0, $c2);
  mgsl_complex_float_rect(0e0, 0e0, $c2);
  mgsl_matrix_complex_float_set($cA2, 1, 1, $c2);
  my gsl_vector_complex_float $cx2 = gsl_vector_complex_float_calloc(2);
  mgsl_complex_float_rect(2e0, 1e0, $c2);
  mgsl_vector_complex_float_set_all($cx2, $c2);
  ok gsl_blas_ctrmv(CblasUpper, CblasNoTrans, CblasNonUnit, $cA2, $cx2) == GSL_SUCCESS, 'op(A)x: complex32';
  is-deeply (gather for ^2 { mgsl_vector_complex_float_get($cx2, $_, $cres); take $cres.dat[0]; take $cres.dat[1] }), (0e0, 10e0, 0e0, 0e0), 'result: ok';

  my gsl_matrix_complex $zA2 = gsl_matrix_complex_calloc(2, 2);
  my gsl_complex $z2 = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 2e0, $z2);
  mgsl_matrix_complex_set($zA2, 0, 0, $z2);
  mgsl_matrix_complex_set($zA2, 0, 1, $z2);
  mgsl_complex_rect(2e0, 1e0, $z2);
  mgsl_matrix_complex_set($zA2, 1, 0, $z2);
  mgsl_complex_rect(0e0, 0e0, $z2);
  mgsl_matrix_complex_set($zA2, 1, 1, $z2);
  my gsl_vector_complex $zx2 = gsl_vector_complex_calloc(2);
  mgsl_complex_rect(2e0, 1e0, $z2);
  mgsl_vector_complex_set_all($zx2, $z2);
  ok gsl_blas_ztrmv(CblasUpper, CblasNoTrans, CblasNonUnit, $zA2, $zx2) == GSL_SUCCESS, 'op(A)x: complex64';
  is-deeply (gather for ^2 { mgsl_vector_complex_get($zx2, $_, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0e0, 10e0, 0e0, 0e0), 'result: ok';

  my gsl_matrix_float $A3 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($A3, 0, 0, .995e0);
  my gsl_vector_float $x3 = gsl_vector_float_calloc(1);
  gsl_vector_float_set_all($x3, .348e0);
  ok gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, $A3, $x3) == GSL_SUCCESS, 'inv op(A)x: num32';
  ok gsl_vector_float_get($x3, 0) == 0.34974873065948486e0, 'result: ok';

  my gsl_matrix $dA3 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dA3, 0, 0, -.21e0);
  my gsl_vector $dx3 = gsl_vector_calloc(1);
  gsl_vector_set($dx3, 0, .473e0);
  ok gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, $dA3, $dx3) == GSL_SUCCESS, 'inv op(A)x: num64';
  ok gsl_vector_get($dx3, 0) == -2.2523809523809524, 'result: ok';

  my gsl_matrix_complex_float $cA3 = gsl_matrix_complex_float_calloc(1, 1);
  my gsl_complex_float $c3 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(.529e0, -.348e0, $c3);
  mgsl_matrix_complex_float_set($cA3, 0, 0, $c3);
  my gsl_vector_complex_float $cx3 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(-.95e0, .343e0, $c3);
  mgsl_vector_complex_float_set_all($cx3, $c3);
  ok gsl_blas_ctrsv(CblasUpper, CblasNoTrans, CblasNonUnit, $cA3, $cx3) == GSL_SUCCESS, 'inv op(A)x: complex32';
  is-deeply (gather { mgsl_vector_complex_float_get($cx3, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-1.551120638847351e0, -0.37200355529785156e0), 'result: ok';

  my gsl_matrix_complex $zA3 = gsl_matrix_complex_calloc(1, 1);
  my gsl_complex $z3 = alloc_gsl_complex;
  mgsl_complex_rect(.977e0, -.955e0, $z3);
  mgsl_matrix_complex_set($zA3, 0, 0, $z3);
  my gsl_vector_complex $zx3 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(-.627e0, .281e0, $z3);
  mgsl_vector_complex_set($zx3, 0, $z3);
  ok gsl_blas_ztrsv(CblasUpper, CblasNoTrans, CblasNonUnit, $zA3, $zx3) == GSL_SUCCESS, 'inv op(A)x: complex64';
  is-deeply (gather { mgsl_vector_complex_get($zx3, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (-0.47195741457252255e0, -0.17371477064151375e0), 'result: ok';

  my $α1 = 1e0; my $β1 = -1e0;
  my gsl_matrix_float $A4 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($A4, 0, 0, -.428e0);
  my gsl_vector_float $x4 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($x4, 0, -.34e0);
  my gsl_vector_float $y4 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($y4, 0, -.888e0);
  ok gsl_blas_ssymv(CblasUpper, $α1, $A4, $x4, $β1, $y4) == GSL_SUCCESS, 'sum αAx + βy: num32';
  ok gsl_vector_float_get($y4, 0) == 1.033519983291626e0, 'result: ok';

  $α1 = 0e0; $β1 = -.3e0;
  my gsl_matrix $dA4 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dA4, 0, 0, .544e0);
  my gsl_vector $dx4 = gsl_vector_calloc(1);
  gsl_vector_set($dx4, 0, -.601e0);
  my gsl_vector $dy4 = gsl_vector_calloc(1);
  gsl_vector_set($dy4, 0, -.852e0);
  ok gsl_blas_dsymv(CblasUpper, $α1, $dA4, $dx4, $β1, $dy4) == GSL_SUCCESS, 'sum αAx + βy: num64';
  ok gsl_vector_get($dy4, 0) == 0.2556e0, 'result: ok';

  mgsl_complex_float_rect(1e0, 0e0, $cα);
  my gsl_complex_float $cβ1 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.3e0, .1e0, $cβ1);
  my gsl_matrix_complex_float $cA4 = gsl_matrix_complex_float_calloc(1, 1);
  my gsl_complex_float $c4 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.434e0, .837e0, $c4);
