Help language development. Donate to The Perl Foundation

Math::Libgsl::BLAS cpan:FRITH last updated on 2020-02-22

t/01-raw.rakutest
#!/usr/bin/env perl6

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);
  mgsl_matrix_complex_float_set($cA4, 0, 0, $c4);
  my gsl_vector_complex_float $cx4 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(.209e0, -.935e0, $c4);
  mgsl_vector_complex_float_set($cx4, 0, $c4);
  my gsl_vector_complex_float $cy4 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(.346e0, -.412e0, $c4);
  mgsl_vector_complex_float_set($cy4, 0, $c4);
  ok mgsl_blas_chemv(CblasUpper, $cα, $cA4, $cx4, $cβ1, $cy4) == GSL_SUCCESS, 'sum αAx + βy: complex32';
  is-deeply (gather { mgsl_vector_complex_float_get($cy4, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-0.1533060073852539e0, 0.5639899969100952e0), 'result: ok';

  mgsl_complex_rect(0e0, 0e0, $zα);
  my gsl_complex $zβ1 = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 0e0, $zβ1);
  my gsl_matrix_complex $zA4 = gsl_matrix_complex_calloc(1, 1);
  my gsl_complex $z4 = alloc_gsl_complex;
  mgsl_complex_rect(.036e0, -.966e0, $z4);
  mgsl_matrix_complex_set($zA4, 0, 0, $z4);
  my gsl_vector_complex $zx4 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(-.695e0, .886e0, $z4);
  mgsl_vector_complex_set($zx4, 0, $z4);
  my gsl_vector_complex $zy4 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(.486e0, .629e0, $z4);
  mgsl_vector_complex_set($zy4, 0, $z4);
  ok mgsl_blas_zhemv(CblasUpper, $zα, $zA4, $zx4, $zβ1, $zy4) == GSL_SUCCESS, 'sum αAx + βy: complex64';
  is-deeply (gather { mgsl_vector_complex_get($zy4, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0.486e0, 0.629e0), 'result: ok';

  $α = 1e0;
  my gsl_matrix_float $A5 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($A5, 0, 0, -.515e0);
  my gsl_vector_float $x5 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($x5, 0, .611e0);
  my gsl_vector_float $y5 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($y5, 0, -.082e0);
  ok gsl_blas_sger($α, $x5, $y5, $A5) == GSL_SUCCESS, 'αxT(y) + A: num32';
  ok gsl_matrix_float_get($A5, 0, 0) == -0.5651019811630249e0, 'result: ok';

  $dα = 1e0;
  my gsl_matrix $dA5 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dA5, 0, 0, -.809e0);
  my gsl_vector $dx5 = gsl_vector_calloc(1);
  gsl_vector_set($dx5, 0, -.652e0);
  my gsl_vector $dy5 = gsl_vector_calloc(1);
  gsl_vector_set($dy5, 0, .712e0);
  ok gsl_blas_dger($dα, $dx5, $dy5, $dA5) == GSL_SUCCESS, 'αxT(y) + A: num64';
  ok gsl_matrix_get($dA5, 0, 0) == -1.273224e0, 'result: ok';

  mgsl_complex_float_rect(0e0, 0e0, $cα);
  my gsl_matrix_complex_float $cA5 = gsl_matrix_complex_float_calloc(1, 1);
  my gsl_complex_float $c5 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.651e0, .856e0, $c5);
  mgsl_matrix_complex_float_set($cA5, 0, 0, $c5);
  my gsl_vector_complex_float $cx5 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(-.38e0, -.235e0, $c5);
  mgsl_vector_complex_float_set($cx5, 0, $c5);
  my gsl_vector_complex_float $cy5 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(-.651e0, .856e0, $c5);
  mgsl_vector_complex_float_set($cy5, 0, $c5);
  ok mgsl_blas_cgeru($cα, $cx5, $cy5, $cA5) == GSL_SUCCESS, 'αxT(y) + A: complex32';
  is-deeply (gather { mgsl_matrix_complex_float_get($cA5, 0, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-0.6510000228881836e0, 0.8560000061988831e0), 'result: ok';

  mgsl_complex_rect(-1e0, 0e0, $zα);
  my gsl_matrix_complex $zA5 = gsl_matrix_complex_calloc(1, 1);
  my gsl_complex $z5 = alloc_gsl_complex;
  mgsl_complex_rect(-.426e0, .757e0, $z5);
  mgsl_matrix_complex_set($zA5, 0, 0, $z5);
  my gsl_vector_complex $zx5 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(-.579e0, -.155e0, $z5);
  mgsl_vector_complex_set($zx5, 0, $z5);
  my gsl_vector_complex $zy5 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(.831e0, .035e0, $z5);
  mgsl_vector_complex_set($zy5, 0, $z5);
  ok mgsl_blas_zgeru($zα, $zx5, $zy5, $zA5) == GSL_SUCCESS, 'αxT(y) + A: complex64';
  is-deeply (gather { mgsl_matrix_complex_get($zA5, 0, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0.049723999999999935e0, 0.90607e0), 'result: ok';

  mgsl_complex_float_rect(0e0, 0e0, $cα);
  my gsl_matrix_complex_float $cA6 = gsl_matrix_complex_float_calloc(1, 1);
  my gsl_complex_float $c6 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.651e0, .856e0, $c6);
  mgsl_matrix_complex_float_set($cA6, 0, 0, $c6);
  my gsl_vector_complex_float $cx6 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(-.38e0, -.235e0, $c6);
  mgsl_vector_complex_float_set($cx6, 0, $c6);
  my gsl_vector_complex_float $cy6 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(-.651e0, .856e0, $c6);
  mgsl_vector_complex_float_set($cy6, 0, $c6);
  ok mgsl_blas_cgerc($cα, $cx6, $cy6, $cA6) == GSL_SUCCESS, 'αxH(y) + A: complex32';
  is-deeply (gather { mgsl_matrix_complex_float_get($cA6, 0, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-0.6510000228881836e0, 0.8560000061988831e0), 'result: ok';

  mgsl_complex_rect(-1e0, 0e0, $zα);
  my gsl_matrix_complex $zA6 = gsl_matrix_complex_calloc(1, 1);
  my gsl_complex $z6 = alloc_gsl_complex;
  mgsl_complex_rect(-.426e0, .757e0, $z6);
  mgsl_matrix_complex_set($zA6, 0, 0, $z6);
  my gsl_vector_complex $zx6 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(-.579e0, -.155e0, $z6);
  mgsl_vector_complex_set($zx6, 0, $z6);
  my gsl_vector_complex $zy6 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(.831e0, .035e0, $z6);
  mgsl_vector_complex_set($zy6, 0, $z6);
  ok mgsl_blas_zgerc($zα, $zx6, $zy6, $zA6) == GSL_SUCCESS, 'αxH(y) + A: complex64';
  is-deeply (gather { mgsl_matrix_complex_get($zA6, 0, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0.06057399999999996e0, 0.86554e0), 'result: ok';

  $α = .1e0;
  my gsl_matrix_float $A7 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($A7, 0, 0, -.291e0);
  my gsl_vector_float $x7 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($x7, 0, .845e0);
  ok gsl_blas_ssyr(CblasUpper, $α, $x7, $A7) == GSL_SUCCESS, 'αxT(x) + A: num32';
  ok gsl_matrix_float_get($A7, 0, 0) == -0.21959750354290009e0, 'result: ok';

  $dα = -.3e0;
  my gsl_matrix $dA7 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dA7, 0, 0, -.65e0);
  my gsl_vector $dx7 = gsl_vector_calloc(1);
  gsl_vector_set($dx7, 0, -.891e0);
  ok gsl_blas_dsyr(CblasUpper, $dα, $dx7, $dA7) == GSL_SUCCESS, 'αxT(x) + A: num64';
  ok gsl_matrix_get($dA7, 0, 0) == -0.8881643e0, 'result: ok';

  my gsl_matrix_complex_float $cA7 = gsl_matrix_complex_float_calloc(1, 1);
  my gsl_complex_float $c7 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.651e0, .856e0, $c7);
  mgsl_matrix_complex_float_set($cA7, 0, 0, $c7);
  my gsl_vector_complex_float $cx7 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(-.38e0, -.235e0, $c7);
  mgsl_vector_complex_float_set($cx7, 0, $c7);
  ok gsl_blas_cher(CblasUpper, $α, $cx7, $cA7) == GSL_SUCCESS, 'αxH(y) + A: complex32';
  is-deeply (gather { mgsl_matrix_complex_float_get($cA7, 0, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-0.6310375332832336e0, 0e0), 'result: ok';

  my gsl_matrix_complex $zA7 = gsl_matrix_complex_calloc(1, 1);
  my gsl_complex $z7 = alloc_gsl_complex;
  mgsl_complex_rect(-.426e0, .757e0, $z7);
  mgsl_matrix_complex_set($zA7, 0, 0, $z7);
  my gsl_vector_complex $zx7 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(-.579e0, -.155e0, $z7);
  mgsl_vector_complex_set($zx7, 0, $z7);
  ok gsl_blas_zher(CblasUpper, $dα, $zx7, $zA7) == GSL_SUCCESS, 'αxH(y) + A: complex64';
  is-deeply (gather { mgsl_matrix_complex_get($zA7, 0, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (-0.5337798e0, 0e0), 'result: ok';

  $α = 0e0;
  my gsl_matrix_float $A8 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($A8, 0, 0, .862e0);
  my gsl_vector_float $x8 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($x8, 0, .823e0);
  my gsl_vector_float $y8 = gsl_vector_float_calloc(1);
  gsl_vector_float_set($y8, 0, .699e0);
  ok gsl_blas_ssyr2(CblasUpper, $α, $x8, $y8, $A8) == GSL_SUCCESS, 'αxT(y) + αyT(x) + A: num32';
  ok gsl_matrix_float_get($A8, 0, 0) == 0.8619999885559082e0, 'result: ok';

  $dα = 0e0;
  my gsl_matrix $dA8 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dA8, 0, 0, -.824e0);
  my gsl_vector $dx8 = gsl_vector_calloc(1);
  gsl_vector_set($dx8, 0, .684e0);
  my gsl_vector $dy8 = gsl_vector_calloc(1);
  gsl_vector_set($dy8, 0, .965e0);
  ok gsl_blas_dsyr2(CblasUpper, $dα, $dx8, $dy8, $dA8) == GSL_SUCCESS, 'αxT(x) + αyT(x) + A: num64';
  ok gsl_matrix_get($dA8, 0, 0) == -0.824e0, 'result: ok';

  mgsl_complex_float_rect(-1e0, 0e0, $cα);
  my gsl_matrix_complex_float $cA8 = gsl_matrix_complex_float_calloc(1, 1);
  my gsl_complex_float $c8 = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.821e0, .954e0, $c8);
  mgsl_matrix_complex_float_set($cA8, 0, 0, $c8);
  my gsl_vector_complex_float $cx8 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(.532e0, .802e0, $c8);
  mgsl_vector_complex_float_set($cx8, 0, $c8);
  my gsl_vector_complex_float $cy8 = gsl_vector_complex_float_calloc(1);
  mgsl_complex_float_rect(.016e0, -.334e0, $c8);
  mgsl_vector_complex_float_set($cy8, 0, $c8);
  ok mgsl_blas_cher2(CblasUpper, $cα, $cx8, $cy8, $cA8) == GSL_SUCCESS, 'αxH(y) + αyH(x) + A: complex32';
  is-deeply (gather { mgsl_matrix_complex_float_get($cA8, 0, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-0.3022879958152771e0, 0e0), 'result: ok';

  mgsl_complex_rect(-.3e0, .1e0, $zα);
  my gsl_matrix_complex $zA8 = gsl_matrix_complex_calloc(1, 1);
  my gsl_complex $z8 = alloc_gsl_complex;
  mgsl_complex_rect(-.334e0, .286e0, $z8);
  mgsl_matrix_complex_set($zA8, 0, 0, $z8);
  my gsl_vector_complex $zx8 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(-.14e0, -.135e0, $z8);
  mgsl_vector_complex_set($zx8, 0, $z8);
  my gsl_vector_complex $zy8 = gsl_vector_complex_calloc(1);
  mgsl_complex_rect(.455e0, .358e0, $z8);
  mgsl_vector_complex_set($zy8, 0, $z8);
  ok mgsl_blas_zher2(CblasUpper, $zα, $zx8, $zy8, $zA8) == GSL_SUCCESS, 'αxH(y) + αyH(x) + A: complex64';
  is-deeply (gather { mgsl_matrix_complex_get($zA8, 0, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (-0.264521e0, 0e0), 'result: ok';
}, 'level 2 functions';

subtest {
  my $α = 1e0; my $β = 0e0;
  my gsl_matrix_float $A = gsl_matrix_float_calloc(1, 4);
  gsl_matrix_float_set($A, 0, 0, .199e0);
  gsl_matrix_float_set($A, 0, 1, .237e0);
  gsl_matrix_float_set($A, 0, 2, .456e0);
  gsl_matrix_float_set($A, 0, 3, .377e0);
  my gsl_matrix_float $B = gsl_matrix_float_calloc(4, 2);
  gsl_matrix_float_set($B, 0, 0,  .842e0);
  gsl_matrix_float_set($B, 0, 1, -.734e0);
  gsl_matrix_float_set($B, 1, 0,  .323e0);
  gsl_matrix_float_set($B, 1, 1, -.957e0);
  gsl_matrix_float_set($B, 2, 0, -.303e0);
  gsl_matrix_float_set($B, 2, 1, -.873e0);
  gsl_matrix_float_set($B, 3, 0, -.871e0);
  gsl_matrix_float_set($B, 3, 1, -.819e0);
  my gsl_matrix_float $C = gsl_matrix_float_calloc(1, 2);
  gsl_matrix_float_set($C, 0, 0,  .498e0);
  gsl_matrix_float_set($C, 0, 1, -.925e0);
  ok gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, $α, $A, $B, $β, $C) == GSL_SUCCESS, 'sum αop(A)op(B) + βC: num32';
  is-deeply (gather take gsl_matrix_float_get($C, 0, $_) for ^2), (-0.22242599725723267e0, -1.0797260999679565e0), 'result: ok';

  $α = 0e0; $β = 0e0;
  my gsl_matrix $dA = gsl_matrix_calloc(1, 4);
  gsl_matrix_set($dA, 0, 0, .017e0);
  gsl_matrix_set($dA, 0, 1, .191e0);
  gsl_matrix_set($dA, 0, 2, .863e0);
  gsl_matrix_set($dA, 0, 3, -.97e0);
  my gsl_matrix $dB = gsl_matrix_calloc(4, 2);
  gsl_matrix_set($dB, 0, 0, -.207e0);
  gsl_matrix_set($dB, 0, 1, -.916e0);
  gsl_matrix_set($dB, 1, 0, -.278e0);
  gsl_matrix_set($dB, 1, 1, .403e0);
  gsl_matrix_set($dB, 2, 0, .885e0);
  gsl_matrix_set($dB, 2, 1, .409e0);
  gsl_matrix_set($dB, 3, 0, -.772e0);
  gsl_matrix_set($dB, 3, 1, -.27e0);
  my gsl_matrix $dC = gsl_matrix_calloc(1, 2);
  gsl_matrix_set($dC, 0, 0, -.274e0);
  gsl_matrix_set($dC, 0, 1, -.858e0);
  ok gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, $α, $dA, $dB, $β, $dC) == GSL_SUCCESS, 'sum αop(A)op(B) + βC: num64';
  is-deeply (gather take gsl_matrix_get($dC, 0, $_) for ^2), (0e0, 0e0), 'result: ok';

  my gsl_complex_float $cα = alloc_gsl_complex_float;
  mgsl_complex_float_rect(0e0, 1e0, $cα);
  my gsl_complex_float $cβ = alloc_gsl_complex_float;
  mgsl_complex_float_rect(0e0, 0e0, $cβ);
  my gsl_matrix_complex_float $cA = gsl_matrix_complex_float_calloc(1, 4);
  my gsl_complex_float $c = alloc_gsl_complex_float;
  mgsl_complex_float_rect(-.082e0, -.281e0, $c);
  mgsl_matrix_complex_float_set($cA, 0, 0, $c);
  mgsl_complex_float_rect(-.096e0, .913e0, $c);
  mgsl_matrix_complex_float_set($cA, 0, 1, $c);
  mgsl_complex_float_rect(.974e0, -.706e0, $c);
  mgsl_matrix_complex_float_set($cA, 0, 2, $c);
  mgsl_complex_float_rect(-.773e0, .522e0, $c);
  mgsl_matrix_complex_float_set($cA, 0, 3, $c);
  my gsl_matrix_complex_float $cB = gsl_matrix_complex_float_calloc(4, 2);
  mgsl_complex_float_rect(.745e0, -.664e0, $c);
  mgsl_matrix_complex_float_set($cB, 0, 0, $c);
  mgsl_complex_float_rect(.352e0, -.733e0, $c);
  mgsl_matrix_complex_float_set($cB, 0, 1, $c);
  mgsl_complex_float_rect(.304e0, -.555e0, $c);
  mgsl_matrix_complex_float_set($cB, 1, 0, $c);
  mgsl_complex_float_rect(-.493e0, -.089e0, $c);
  mgsl_matrix_complex_float_set($cB, 1, 1, $c);
  mgsl_complex_float_rect(.188e0, .631e0, $c);
  mgsl_matrix_complex_float_set($cB, 2, 0, $c);
  mgsl_complex_float_rect(.235e0, .152e0, $c);
  mgsl_matrix_complex_float_set($cB, 2, 1, $c);
  mgsl_complex_float_rect(-.299e0, -.731e0, $c);
  mgsl_matrix_complex_float_set($cB, 3, 0, $c);
  mgsl_complex_float_rect(-.686e0, -.332e0, $c);
  mgsl_matrix_complex_float_set($cB, 3, 1, $c);
  my gsl_matrix_complex_float $cC = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(-.179e0, -.284e0, $c);
  mgsl_matrix_complex_float_set($cC, 0, 0, $c);
  mgsl_complex_float_rect(-.996e0, -.414e0, $c);
  mgsl_matrix_complex_float_set($cC, 0, 1, $c);
  ok mgsl_blas_cgemm(CblasNoTrans, CblasNoTrans, $cα, $cA, $cB, $cβ, $cC) == GSL_SUCCESS, 'sum αop(A)op(B) + βC: complex32';
  my gsl_complex_float $cres = alloc_gsl_complex_float;
  is-deeply (gather for ^2 { mgsl_matrix_complex_float_get($cC, 0, $_, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-1.0667860507965088e0, 1.4711639881134033e0, 0.5996890068054199e0, 0.9335318803787231e0), 'result: ok';

  my gsl_complex $zα = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 0e0, $zα);
  my gsl_complex $zβ = alloc_gsl_complex;
  mgsl_complex_rect(-1e0, 0e0, $zβ);
  my gsl_matrix_complex $zA = gsl_matrix_complex_calloc(1, 4);
  my gsl_complex $z = alloc_gsl_complex;
  mgsl_complex_rect(-.33e0, .457e0, $z);
  mgsl_matrix_complex_set($zA, 0, 0, $z);
  mgsl_complex_rect(.428e0, -.19e0, $z);
  mgsl_matrix_complex_set($zA, 0, 1, $z);
  mgsl_complex_rect(.86e0, -.53e0, $z);
  mgsl_matrix_complex_set($zA, 0, 2, $z);
  mgsl_complex_rect(.058e0, -.942e0, $z);
  mgsl_matrix_complex_set($zA, 0, 3, $z);
  my gsl_matrix_complex $zB = gsl_matrix_complex_calloc(4, 2);
  mgsl_complex_rect(.434e0, .653e0, $z);
  mgsl_matrix_complex_set($zB, 0, 0, $z);
  mgsl_complex_rect(-.124e0, .191e0, $z);
  mgsl_matrix_complex_set($zB, 0, 1, $z);
  mgsl_complex_rect(-.112e0, -.84e0, $z);
  mgsl_matrix_complex_set($zB, 1, 0, $z);
  mgsl_complex_rect(-.72e0, .075e0, $z);
  mgsl_matrix_complex_set($zB, 1, 1, $z);
  mgsl_complex_rect(-.503e0, -.109e0, $z);
  mgsl_matrix_complex_set($zB, 2, 0, $z);
  mgsl_complex_rect(.3e0, -.898e0, $z);
  mgsl_matrix_complex_set($zB, 2, 1, $z);
  mgsl_complex_rect(.489e0, .384e0, $z);
  mgsl_matrix_complex_set($zB, 3, 0, $z);
  mgsl_complex_rect(.993e0, -.804e0, $z);
  mgsl_matrix_complex_set($zB, 3, 1, $z);
  my gsl_matrix_complex $zC = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(-.792e0, -.155e0, $z);
  mgsl_matrix_complex_set($zC, 0, 0, $z);
  mgsl_complex_rect(-.608e0, -.243e0, $z);
  mgsl_matrix_complex_set($zC, 0, 1, $z);
  ok mgsl_blas_zgemm(CblasNoTrans, CblasNoTrans, $zα, $zA, $zB, $zβ, $zC) == GSL_SUCCESS, 'sum αop(A)op(B) + βC: complex64';
  my gsl_complex $zres = alloc_gsl_complex;
  is-deeply (gather for ^2 { mgsl_matrix_complex_get($zC, 0, $_, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0.04256299999999996e0, -0.465908e0, -0.6499910000000001e0, -1.621116e0), 'result: ok';

  $α = -.3e0; $β = -1e0;
  my gsl_matrix_float $A1 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($A1, 0, 0, -.581e0);
  my gsl_matrix_float $B1 = gsl_matrix_float_calloc(1, 2);
  gsl_matrix_float_set($B1, 0, 0, .157e0);
  gsl_matrix_float_set($B1, 0, 1, .451e0);
  my gsl_matrix_float $C1 = gsl_matrix_float_calloc(1, 2);
  gsl_matrix_float_set($C1, 0, 0, -.869e0);
  gsl_matrix_float_set($C1, 0, 1, -.871e0);
  ok gsl_blas_ssymm(CblasLeft, CblasUpper, $α, $A1, $B1, $β, $C1) == GSL_SUCCESS, 'αAB + βC: num32';
  is-deeply (gather take gsl_matrix_float_get($C1, 0, $_) for ^2), (0.8963651061058044e0, 0.9496092796325684e0), 'result: ok';

  $α = -.3e0; $β = 1e0;
  my gsl_matrix $dA1 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dA1, 0, 0, -.981e0);
  my gsl_matrix $dB1 = gsl_matrix_calloc(1, 2);
  gsl_matrix_set($dB1, 0, 0, -.823e0);
  gsl_matrix_set($dB1, 0, 1, .83e0);
  my gsl_matrix $dC1 = gsl_matrix_calloc(1, 2);
  gsl_matrix_set($dC1, 0, 0, .991e0);
  gsl_matrix_set($dC1, 0, 1, .382e0);
  ok gsl_blas_dsymm(CblasLeft, CblasUpper, $α, $dA1, $dB1, $β, $dC1) == GSL_SUCCESS, 'αAB + βC: num64';
  is-deeply (gather take gsl_matrix_get($dC1, 0, $_) for ^2), (0.7487911e0, 0.626269e0), 'result: ok';

  mgsl_complex_float_rect(-1e0, 0e0, $cα);
  mgsl_complex_float_rect(1e0, 0e0, $cβ);
  my gsl_matrix_complex_float $cA1 = gsl_matrix_complex_float_calloc(1, 1);
  mgsl_complex_float_rect(.476e0, .816e0, $c);
  mgsl_matrix_complex_float_set($cA1, 0, 0, $c);
  my gsl_matrix_complex_float $cB1 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(.282e0, .852e0, $c);
  mgsl_matrix_complex_float_set($cB1, 0, 0, $c);
  mgsl_complex_float_rect(-.891e0, -.588e0, $c);
  mgsl_matrix_complex_float_set($cB1, 0, 1, $c);
  my gsl_matrix_complex_float $cC1 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(.9e0, .486e0, $c);
  mgsl_matrix_complex_float_set($cC1, 0, 0, $c);
  mgsl_complex_float_rect(-.78e0, -.637e0, $c);
  mgsl_matrix_complex_float_set($cC1, 0, 1, $c);
  ok mgsl_blas_csymm(CblasLeft, CblasUpper, $cα, $cA1, $cB1, $cβ, $cC1) == GSL_SUCCESS, 'αAB + βC: complex32';
  is-deeply (gather for ^2 { mgsl_matrix_complex_float_get($cC1, 0, $_, $cres); take $cres.dat[0]; take $cres.dat[1] }), (1.4609999656677246e0, -0.1496639847755432e0, -0.8356919884681702e0, 0.36994391679763794e0), 'result: ok';

  mgsl_complex_rect(0e0, 0e0, $zα);
  mgsl_complex_rect(0e0, 0e0, $zβ);
  my gsl_matrix_complex $zA1 = gsl_matrix_complex_calloc(1, 1);
  mgsl_complex_rect(.511e0, -.486e0, $z);
  mgsl_matrix_complex_set($zA1, 0, 0, $z);
  my gsl_matrix_complex $zB1 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(.985e0, -.923e0, $z);
  mgsl_matrix_complex_set($zB1, 0, 0, $z);
  mgsl_complex_rect(-.234e0, -.756e0, $z);
  mgsl_matrix_complex_set($zB1, 0, 1, $z);
  my gsl_matrix_complex $zC1 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(-.16e0, .049e0, $z);
  mgsl_matrix_complex_set($zC1, 0, 0, $z);
  mgsl_complex_rect(.618e0, -.349e0, $z);
  mgsl_matrix_complex_set($zC1, 0, 1, $z);
  ok mgsl_blas_zsymm(CblasLeft, CblasUpper, $zα, $zA1, $zB1, $zβ, $zC1) == GSL_SUCCESS, 'αAB + βC: complex64';
  is-deeply (gather for ^2 { mgsl_matrix_complex_get($zC1, 0, $_, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0e0, 0e0, 0e0, 0e0), 'result: ok';

  mgsl_complex_float_rect(0e0, .1e0, $cα);
  mgsl_complex_float_rect(0e0, .1e0, $cβ);
  my gsl_matrix_complex_float $cA2 = gsl_matrix_complex_float_calloc(1, 1);
  mgsl_complex_float_rect(-.126e0, .079e0, $c);
  mgsl_matrix_complex_float_set($cA2, 0, 0, $c);
  my gsl_matrix_complex_float $cB2 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(-.954e0, -.059e0, $c);
  mgsl_matrix_complex_float_set($cB2, 0, 0, $c);
  mgsl_complex_float_rect(.296e0, -.988e0, $c);
  mgsl_matrix_complex_float_set($cB2, 0, 1, $c);
  my gsl_matrix_complex_float $cC2 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(-.859e0, -.731e0, $c);
  mgsl_matrix_complex_float_set($cC2, 0, 0, $c);
  mgsl_complex_float_rect(.737e0, .593e0, $c);
  mgsl_matrix_complex_float_set($cC2, 0, 1, $c);
  ok mgsl_blas_chemm(CblasLeft, CblasUpper, $cα, $cA2, $cB2, $cβ, $cC2) == GSL_SUCCESS, 'sum αAB + βC: complex32';
  is-deeply (gather for ^2 { mgsl_matrix_complex_float_get($cC2, 0, $_, $cres); take $cres.dat[0]; take $cres.dat[1] }), (0.07235660403966904e0, -0.07387959957122803e0, -0.07174880057573318e0, 0.06997040659189224e0), 'result: ok';

  mgsl_complex_rect(0e0, .1e0, $zα);
  mgsl_complex_rect(0e0, .1e0, $zβ);
  my gsl_matrix_complex $zA2 = gsl_matrix_complex_calloc(1, 1);
  mgsl_complex_rect(-.359e0, .089e0, $z);
  mgsl_matrix_complex_set($zA2, 0, 0, $z);
  my gsl_matrix_complex $zB2 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(-.451e0, -.337e0, $z);
  mgsl_matrix_complex_set($zB2, 0, 0, $z);
  mgsl_complex_rect(-.901e0, -.871e0, $z);
  mgsl_matrix_complex_set($zB2, 0, 1, $z);
  my gsl_matrix_complex $zC2 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(.729e0, .631e0, $z);
  mgsl_matrix_complex_set($zC2, 0, 0, $z);
  mgsl_complex_rect(.364e0, .246e0, $z);
  mgsl_matrix_complex_set($zC2, 0, 1, $z);
  ok mgsl_blas_zhemm(CblasLeft, CblasUpper, $zα, $zA2, $zB2, $zβ, $zC2) == GSL_SUCCESS, 'sum αAB + βC: complex64';
  is-deeply (gather for ^2 { mgsl_matrix_complex_get($zC2, 0, $_, $zres); take $zres.dat[0]; take $zres.dat[1] }), (-0.0751983e0, 0.0890909e0, -0.0558689e0, 0.0687459e0), 'result: ok';

  $α = -.3e0;
  my gsl_matrix_float $A3 = gsl_matrix_float_calloc(2, 2);
  gsl_matrix_float_set($A3, 0, 0, .18e0);
  gsl_matrix_float_set($A3, 0, 1, .199e0);
  gsl_matrix_float_set($A3, 1, 0, .122e0);
  gsl_matrix_float_set($A3, 1, 1, -.547e0);
  my gsl_matrix_float $B3 = gsl_matrix_float_calloc(2, 3);
  gsl_matrix_float_set($B3, 0, 0, -.874e0);
  gsl_matrix_float_set($B3, 0, 1, -.383e0);
  gsl_matrix_float_set($B3, 0, 2, .458e0);
  gsl_matrix_float_set($B3, 1, 0, .124e0);
  gsl_matrix_float_set($B3, 1, 1, -.221e0);
  gsl_matrix_float_set($B3, 1, 2, -.107e0);
  ok gsl_blas_strmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $A3, $B3) == GSL_SUCCESS, 'αop(A)B: num32';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take gsl_matrix_float_get($B3, $i, $j) } }), (0.03979320451617241e0, 0.03387570381164551e0, -0.018344102427363396e0, 0.020348399877548218e0, -0.03626609966158867e0, -0.01755870133638382e0), 'result: ok';

  $α = -.3e0;
  my gsl_matrix $dA3 = gsl_matrix_calloc(2, 2);
  gsl_matrix_set($dA3, 0, 0, .174e0);
  gsl_matrix_set($dA3, 0, 1, -.308e0);
  gsl_matrix_set($dA3, 1, 0, .997e0);
  gsl_matrix_set($dA3, 1, 1, -.484e0);
  my gsl_matrix $dB3 = gsl_matrix_calloc(2, 3);
  gsl_matrix_set($dB3, 0, 0, -.256e0);
  gsl_matrix_set($dB3, 0, 1, -.178e0);
  gsl_matrix_set($dB3, 0, 2, .098e0);
  gsl_matrix_set($dB3, 1, 0, .004e0);
  gsl_matrix_set($dB3, 1, 1, .97e0);
  gsl_matrix_set($dB3, 1, 2, -.408e0);
  ok gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $dA3, $dB3) == GSL_SUCCESS, 'αop(A)B: num64';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take gsl_matrix_get($dB3, $i, $j) } }), (0.013732799999999998e0, 0.09891959999999998e0, -0.0428148e0, 0.0005808e0, 0.14084399999999997e0, -0.05924159999999999e0), 'result: ok';

  mgsl_complex_float_rect(-1e0, 0e0, $cα);
  my gsl_matrix_complex_float $cA3 = gsl_matrix_complex_float_calloc(2, 2);
  mgsl_complex_float_rect(-.023e0, .762e0, $c);
  mgsl_matrix_complex_float_set($cA3, 0, 0, $c);
  mgsl_complex_float_rect(-.687e0, -.039e0, $c);
  mgsl_matrix_complex_float_set($cA3, 0, 1, $c);
  mgsl_complex_float_rect(-.459e0, .047e0, $c);
  mgsl_matrix_complex_float_set($cA3, 1, 0, $c);
  mgsl_complex_float_rect(.189e0, .33e0, $c);
  mgsl_matrix_complex_float_set($cA3, 1, 1, $c);
  my gsl_matrix_complex_float $cB3 = gsl_matrix_complex_float_calloc(2, 3);
  mgsl_complex_float_rect(.827e0, -.561e0, $c);
  mgsl_matrix_complex_float_set($cB3, 0, 0, $c);
  mgsl_complex_float_rect(.641e0, -.229e0, $c);
  mgsl_matrix_complex_float_set($cB3, 0, 1, $c);
  mgsl_complex_float_rect(-.884e0, -.533e0, $c);
  mgsl_matrix_complex_float_set($cB3, 0, 2, $c);
  mgsl_complex_float_rect(-.624e0, -.138e0, $c);
  mgsl_matrix_complex_float_set($cB3, 1, 0, $c);
  mgsl_complex_float_rect(.073e0, .924e0, $c);
  mgsl_matrix_complex_float_set($cB3, 1, 1, $c);
  mgsl_complex_float_rect(-.501e0, -.164e0, $c);
  mgsl_matrix_complex_float_set($cB3, 1, 2, $c);
  ok mgsl_blas_ctrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $cα, $cA3, $cB3) == GSL_SUCCESS, 'αop(A)B: complex32';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { mgsl_matrix_complex_float_get($cB3, $i, $j, $cres); take $cres.dat[0]; take $cres.dat[1] } }), (-0.8317669630050659e0, -0.7622190117835999e0, -0.14564001560211182e0, 0.1439259946346283e0, -0.7642689943313599e0, 0.5291420221328735e0, 0.07239599525928497e0, 0.23200200498104095e0, 0.29112303256988525e0, -0.19872601330280304e0, 0.040568992495536804e0, 0.19632600247859955e0), 'result: ok';

  mgsl_complex_rect(0e0, 0e0, $zα);
  my gsl_matrix_complex $zA3 = gsl_matrix_complex_calloc(2, 2);
  mgsl_complex_rect(.463e0, .033e0, $z);
  mgsl_matrix_complex_set($zA3, 0, 0, $z);
  mgsl_complex_rect(-.929e0, .949e0, $z);
  mgsl_matrix_complex_set($zA3, 0, 1, $z);
  mgsl_complex_rect(.864e0, .986e0, $z);
  mgsl_matrix_complex_set($zA3, 1, 0, $z);
  mgsl_complex_rect(.393e0, .885e0, $z);
  mgsl_matrix_complex_set($zA3, 1, 1, $z);
  my gsl_matrix_complex $zB3 = gsl_matrix_complex_calloc(2, 3);
  mgsl_complex_rect(-.321e0, -.852e0, $z);
  mgsl_matrix_complex_set($zB3, 0, 0, $z);
  mgsl_complex_rect(-.337e0, -.175e0, $z);
  mgsl_matrix_complex_set($zB3, 0, 1, $z);
  mgsl_complex_rect(.607e0, -.613e0, $z);
  mgsl_matrix_complex_set($zB3, 0, 2, $z);
  mgsl_complex_rect(.688e0, .973e0, $z);
  mgsl_matrix_complex_set($zB3, 1, 0, $z);
  mgsl_complex_rect(-.331e0, -.35e0, $z);
  mgsl_matrix_complex_set($zB3, 1, 1, $z);
  mgsl_complex_rect(.719e0, -.553e0, $z);
  mgsl_matrix_complex_set($zB3, 1, 2, $z);
  ok mgsl_blas_ztrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $zα, $zA3, $zB3) == GSL_SUCCESS, 'αop(A)B: complex64';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { mgsl_matrix_complex_get($zB3, $i, $j, $zres); take $zres.dat[0]; take $zres.dat[1] } } ), (0e0, -0e0, 0e0, 0e0, 0e0, 0e0, -0e0, 0e0, 0e0, 0e0, 0e0, 0e0), 'result: ok';

  $α = -.3e0;
  my gsl_matrix_float $A4 = gsl_matrix_float_calloc(2, 2);
  gsl_matrix_float_set($A4, 0, 0, -.279e0);
  gsl_matrix_float_set($A4, 0, 1,  .058e0);
  gsl_matrix_float_set($A4, 1, 0,  .437e0);
  gsl_matrix_float_set($A4, 1, 1,  .462e0);
  my gsl_matrix_float $B4 = gsl_matrix_float_calloc(2, 3);
  gsl_matrix_float_set($B4, 0, 0,  .578e0);
  gsl_matrix_float_set($B4, 0, 1,  .473e0);
  gsl_matrix_float_set($B4, 0, 2,  -.34e0);
  gsl_matrix_float_set($B4, 1, 0, -.128e0);
  gsl_matrix_float_set($B4, 1, 1,  .503e0);
  gsl_matrix_float_set($B4, 1, 2,    .2e0);
  ok gsl_blas_strsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $A4, $B4) == GSL_SUCCESS, 'αop(inv(A))B: num32';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take gsl_matrix_float_get($B4, $i, $j) } }), (0.6387841701507568e0, 0.4407019317150116e0, -0.3925895094871521e0, 0.08311688154935837e0, -0.3266233801841736e0, -0.12987013161182404e0), 'result: ok';

  $α = -1e0;
  my gsl_matrix $dA4 = gsl_matrix_calloc(2, 2);
  gsl_matrix_set($dA4, 0, 0, -.876e0);
  gsl_matrix_set($dA4, 0, 1, -.503e0);
  gsl_matrix_set($dA4, 1, 0, -.062e0);
  gsl_matrix_set($dA4, 1, 1, -.987e0);
  my gsl_matrix $dB4 = gsl_matrix_calloc(2, 3);
  gsl_matrix_set($dB4, 0, 0,  .219e0);
  gsl_matrix_set($dB4, 0, 1, -.986e0);
  gsl_matrix_set($dB4, 0, 2,    -0e0);
  gsl_matrix_set($dB4, 1, 0, -.605e0);
  gsl_matrix_set($dB4, 1, 1,  .289e0);
  gsl_matrix_set($dB4, 1, 2,  .641e0);
  ok gsl_blas_dtrsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $α, $dA4, $dB4) == GSL_SUCCESS, 'αop(inv(A))B: num64';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { take gsl_matrix_get($dB4, $i, $j) } }), (0.6019671251382123e0, -1.2937005269415647e0, -0.3729106234935439e0, -0.6129685916919959e0, 0.292806484295846e0, 0.6494427558257345e0), 'result: ok';

  mgsl_complex_float_rect(0e0, 0e0, $cα);
  my gsl_matrix_complex_float $cA4 = gsl_matrix_complex_float_calloc(2, 2);
  mgsl_complex_float_rect(.491e0, -.317e0, $c);
  mgsl_matrix_complex_float_set($cA4, 0, 0, $c);
  mgsl_complex_float_rect(-.14e0, -.739e0, $c);
  mgsl_matrix_complex_float_set($cA4, 0, 1, $c);
  mgsl_complex_float_rect(-.969e0, -.518e0, $c);
  mgsl_matrix_complex_float_set($cA4, 1, 0, $c);
  mgsl_complex_float_rect(.702e0, -.287e0, $c);
  mgsl_matrix_complex_float_set($cA4, 1, 1, $c);
  my gsl_matrix_complex_float $cB4 = gsl_matrix_complex_float_calloc(2, 3);
  mgsl_complex_float_rect(-.962e0, -.38e0, $c);
  mgsl_matrix_complex_float_set($cB4, 0, 0, $c);
  mgsl_complex_float_rect(.656e0, .587e0, $c);
  mgsl_matrix_complex_float_set($cB4, 0, 1, $c);
  mgsl_complex_float_rect(-.195e0, -.862e0, $c);
  mgsl_matrix_complex_float_set($cB4, 0, 2, $c);
  mgsl_complex_float_rect(-.679e0, .598e0, $c);
  mgsl_matrix_complex_float_set($cB4, 1, 0, $c);
  mgsl_complex_float_rect(.919e0, .714e0, $c);
  mgsl_matrix_complex_float_set($cB4, 1, 1, $c);
  mgsl_complex_float_rect(-.513e0, .726e0, $c);
  mgsl_matrix_complex_float_set($cB4, 1, 2, $c);
  ok mgsl_blas_ctrsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $cα, $cA4, $cB4) == GSL_SUCCESS, 'αop(inv(A))B: complex32';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { mgsl_matrix_complex_float_get($cB4, $i, $j, $cres); take $cres.dat[0]; take $cres.dat[1] } }), (0e0, 0e0, 0e0, 0e0, 0e0, 0e0, -0e0, 0e0, 0e0, 0e0, -0e0, 0e0), 'result: ok';

  mgsl_complex_rect(0e0, 0e0, $zα);
  my gsl_matrix_complex $zA4 = gsl_matrix_complex_calloc(2, 2);
  mgsl_complex_rect(.189e0, .519e0, $z);
  mgsl_matrix_complex_set($zA4, 0, 0, $z);
  mgsl_complex_rect(-.455e0, -.444e0, $z);
  mgsl_matrix_complex_set($zA4, 0, 1, $z);
  mgsl_complex_rect(-.21e0, -.507e0, $z);
  mgsl_matrix_complex_set($zA4, 1, 0, $z);
  mgsl_complex_rect(-.591e0, .859e0, $z);
  mgsl_matrix_complex_set($zA4, 1, 1, $z);
  my gsl_matrix_complex $zB4 = gsl_matrix_complex_calloc(2, 3);
  mgsl_complex_rect(-.779e0, -.484e0, $z);
  mgsl_matrix_complex_set($zB4, 0, 0, $z);
  mgsl_complex_rect(.249e0, -.107e0, $z);
  mgsl_matrix_complex_set($zB4, 0, 1, $z);
  mgsl_complex_rect(-.755e0, -.047e0, $z);
  mgsl_matrix_complex_set($zB4, 0, 2, $z);
  mgsl_complex_rect(.941e0, .675e0, $z);
  mgsl_matrix_complex_set($zB4, 1, 0, $z);
  mgsl_complex_rect(-.757e0, .645e0, $z);
  mgsl_matrix_complex_set($zB4, 1, 1, $z);
  mgsl_complex_rect(-.649e0, .242e0, $z);
  mgsl_matrix_complex_set($zB4, 1, 2, $z);
  ok mgsl_blas_ztrsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, $zα, $zA4, $zB4) == GSL_SUCCESS, 'αop(inv(A))B: complex64';
  is-deeply (gather for ^2 -> $i { for ^3 -> $j { mgsl_matrix_complex_get($zB4, $i, $j, $zres); take $zres.dat[0]; take $zres.dat[1] } } ), (0e0, -0e0, 0e0, 0e0, 0e0, 0e0, 0e0, -0e0, 0e0, 0e0, 0e0, 0e0), 'result: ok';

  $α = -1e0; $β = .1e0;
  my gsl_matrix_float $A5 = gsl_matrix_float_calloc(2, 1);
  gsl_matrix_float_set($A5, 0, 0,  .412e0);
  gsl_matrix_float_set($A5, 1, 0, -.229e0);
  my gsl_matrix_float $B5 = gsl_matrix_float_calloc(2, 2);
  gsl_matrix_float_set($B5, 0, 0,  .628e0);
  gsl_matrix_float_set($B5, 0, 1, -.664e0);
  gsl_matrix_float_set($B5, 1, 0, -.268e0);
  gsl_matrix_float_set($B5, 1, 1,  .096e0);
  ok gsl_blas_ssyrk(CblasUpper, CblasNoTrans, $α, $A5, $β, $B5) == GSL_SUCCESS, 'αAT(A) + βB: num32';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take gsl_matrix_float_get($B5, $i, $j) } }), (-0.1069439947605133e0, 0.02794799953699112e0, -0.2680000066757202e0, -0.042841002345085144e0), 'result: ok';

  $α = .1e0; $β = -.3e0;
  my gsl_matrix $dA5 = gsl_matrix_calloc(2, 1);
  gsl_matrix_set($dA5, 0, 0,  .169e0);
  gsl_matrix_set($dA5, 1, 0, -.875e0);
  my gsl_matrix $dB5 = gsl_matrix_calloc(2, 2);
  gsl_matrix_set($dB5, 0, 0, .159e0);
  gsl_matrix_set($dB5, 0, 1, .277e0);
  gsl_matrix_set($dB5, 1, 0, .865e0);
  gsl_matrix_set($dB5, 1, 1, .346e0);
  ok gsl_blas_dsyrk(CblasUpper, CblasNoTrans, $α, $dA5, $β, $dB5) == GSL_SUCCESS, 'αAT(A) + βB: num64';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { take gsl_matrix_get($dB5, $i, $j) } }), (-0.0448439e0, -0.09788750000000002e0, 0.865e0, -0.027237499999999984e0), 'result: ok';

  mgsl_complex_float_rect(0e0, 0e0, $cα);
  mgsl_complex_float_rect(-.3e0, .1e0, $cβ);
  my gsl_matrix_complex_float $cA5 = gsl_matrix_complex_float_calloc(2, 1);
  mgsl_complex_float_rect(.447e0, -.507e0, $c);
  mgsl_matrix_complex_float_set($cA5, 0, 0, $c);
  mgsl_complex_float_rect(-.425e0, .701e0, $c);
  mgsl_matrix_complex_float_set($cA5, 1, 0, $c);
  my gsl_matrix_complex_float $cB5 = gsl_matrix_complex_float_calloc(2, 2);
  mgsl_complex_float_rect(.16e0, -.245e0, $c);
  mgsl_matrix_complex_float_set($cB5, 0, 0, $c);
  mgsl_complex_float_rect(.922e0, -.437e0, $c);
  mgsl_matrix_complex_float_set($cB5, 0, 1, $c);
  mgsl_complex_float_rect(.24e0, .008e0, $c);
  mgsl_matrix_complex_float_set($cB5, 1, 0, $c);
  mgsl_complex_float_rect(-.095e0, .749e0, $c);
  mgsl_matrix_complex_float_set($cB5, 1, 1, $c);
  ok mgsl_blas_csyrk(CblasUpper, CblasNoTrans, $cα, $cA5, $cβ, $cB5) == GSL_SUCCESS, 'αAT(A) + βB: complex32';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { mgsl_matrix_complex_float_get($cB5, $i, $j, $cres); take $cres.dat[0]; take $cres.dat[1] } }), (-0.023499999195337296e0, 0.08950001001358032e0, -0.2328999936580658e0, 0.2233000099658966e0, 0.23999999463558197e0, 0.00800000037997961e0, -0.046400003135204315e0, -0.23420001566410065e0), 'result: ok';

  mgsl_complex_rect(1e0, 0e0, $zα);
  mgsl_complex_rect(1e0, 0e0, $zβ);
  my gsl_matrix_complex $zA5 = gsl_matrix_complex_calloc(2, 1);
  mgsl_complex_rect(-.049e0, -.687e0, $z);
  mgsl_matrix_complex_set($zA5, 0, 0, $z);
  mgsl_complex_rect(-.434e0, .294e0, $z);
  mgsl_matrix_complex_set($zA5, 1, 0, $z);
  my gsl_matrix_complex $zB5 = gsl_matrix_complex_calloc(2, 2);
  mgsl_complex_rect(.937e0, -.113e0, $z);
  mgsl_matrix_complex_set($zB5, 0, 0, $z);
  mgsl_complex_rect(.796e0, .293e0, $z);
  mgsl_matrix_complex_set($zB5, 0, 1, $z);
  mgsl_complex_rect(.876e0, -.199e0, $z);
  mgsl_matrix_complex_set($zB5, 1, 0, $z);
  mgsl_complex_rect(-.757e0, -.103e0, $z);
  mgsl_matrix_complex_set($zB5, 1, 1, $z);
  ok mgsl_blas_zsyrk(CblasUpper, CblasNoTrans, $zα, $zA5, $zβ, $zB5) == GSL_SUCCESS, 'αAT(A) + βB: complex64';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { mgsl_matrix_complex_get($zB5, $i, $j, $zres); take $zres.dat[0]; take $zres.dat[1] } } ), (0.46743199999999996e0, -0.04567399999999999e0, 1.019244e0, 0.576752e0, 0.876e0, -0.199e0, -0.65508e0, -0.35819199999999995e0), 'result: ok';

  $α = 0e0; $β = .1e0;
  my gsl_matrix_complex_float $cA6 = gsl_matrix_complex_float_calloc(2, 1);
  mgsl_complex_float_rect(-.617e0, .179e0, $c);
  mgsl_matrix_complex_float_set($cA6, 0, 0, $c);
  mgsl_complex_float_rect(-.626e0, .334e0, $c);
  mgsl_matrix_complex_float_set($cA6, 1, 0, $c);
  my gsl_matrix_complex_float $cB6 = gsl_matrix_complex_float_calloc(2, 2);
  mgsl_complex_float_rect(.346e0, -.903e0, $c);
  mgsl_matrix_complex_float_set($cB6, 0, 0, $c);
  mgsl_complex_float_rect(.022e0, -.839e0, $c);
  mgsl_matrix_complex_float_set($cB6, 0, 1, $c);
  mgsl_complex_float_rect(-.715e0, .049e0, $c);
  mgsl_matrix_complex_float_set($cB6, 1, 0, $c);
  mgsl_complex_float_rect(-.338e0, .149e0, $c);
  mgsl_matrix_complex_float_set($cB6, 1, 1, $c);
  ok gsl_blas_cherk(CblasUpper, CblasNoTrans, $α, $cA6, $β, $cB6) == GSL_SUCCESS, 'αAH(A) + βB: complex32';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { mgsl_matrix_complex_float_get($cB6, $i, $j, $cres); take $cres.dat[0]; take $cres.dat[1] } }), (0.03460000082850456e0, 0e0, 0.002199999988079071e0, -0.08389999717473984e0, -0.7149999737739563e0, 0.04899999871850014e0, -0.03380000218749046e0, 0e0), 'result: ok';

  $α = 1e0; $β = 1e0;
  my gsl_matrix_complex $zA6 = gsl_matrix_complex_calloc(2, 1);
  mgsl_complex_rect(.934e0, .664e0, $z);
  mgsl_matrix_complex_set($zA6, 0, 0, $z);
  mgsl_complex_rect(.426e0, .263e0, $z);
  mgsl_matrix_complex_set($zA6, 1, 0, $z);
  my gsl_matrix_complex $zB6 = gsl_matrix_complex_calloc(2, 2);
  mgsl_complex_rect(.251e0, -.97e0, $z);
  mgsl_matrix_complex_set($zB6, 0, 0, $z);
  mgsl_complex_rect(.76e0, -.349e0, $z);
  mgsl_matrix_complex_set($zB6, 0, 1, $z);
  mgsl_complex_rect(.152e0, -.899e0, $z);
  mgsl_matrix_complex_set($zB6, 1, 0, $z);
  mgsl_complex_rect(-.17e0, .707e0, $z);
  mgsl_matrix_complex_set($zB6, 1, 1, $z);
  ok gsl_blas_zherk(CblasUpper, CblasNoTrans, $α, $zA6, $β, $zB6) == GSL_SUCCESS, 'αAH(A) + βB: complex64';
  is-deeply (gather for ^2 -> $i { for ^2 -> $j { mgsl_matrix_complex_get($zB6, $i, $j, $zres); take $zres.dat[0]; take $zres.dat[1] } } ), (1.5642520000000002e0, 0e0, 1.332516e0, -0.311778e0, 0.152e0, -0.899e0, 0.080645e0, 0e0), 'result: ok';

  $α = .1e0; $β = 1e0;
  my gsl_matrix_float $A7 = gsl_matrix_float_calloc(1, 2);
  gsl_matrix_float_set($A7, 0, 0, -.915e0);
  gsl_matrix_float_set($A7, 0, 1,  .445e0);
  my gsl_matrix_float $B7 = gsl_matrix_float_calloc(1, 2);
  gsl_matrix_float_set($B7, 0, 0,  .213e0);
  gsl_matrix_float_set($B7, 0, 1, -.194e0);
  my gsl_matrix_float $C7 = gsl_matrix_float_calloc(1, 1);
  gsl_matrix_float_set($C7, 0, 0, -.117e0);
  ok gsl_blas_ssyr2k(CblasUpper, CblasNoTrans, $α, $A7, $B7, $β, $C7) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC: num32';
  ok gsl_matrix_float_get($C7, 0, 0) == -0.17324499785900116e0, 'result: ok';

  $α = .1e0; $β = 0e0;
  my gsl_matrix $dA7 = gsl_matrix_calloc(1, 2);
  gsl_matrix_set($dA7, 0, 0, -.225e0);
  gsl_matrix_set($dA7, 0, 1,  .857e0);
  my gsl_matrix $dB7 = gsl_matrix_calloc(1, 2);
  gsl_matrix_set($dB7, 0, 0, -.933e0);
  gsl_matrix_set($dB7, 0, 1,  .994e0);
  my gsl_matrix $dC7 = gsl_matrix_calloc(1, 1);
  gsl_matrix_set($dC7, 0, 0, .177e0);
  ok gsl_blas_dsyr2k(CblasUpper, CblasNoTrans, $α, $dA7, $dB7, $β, $dC7) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC: num64';
  ok gsl_matrix_get($dC7, 0, 0) == 0.21235660000000003e0, 'result: ok';

  mgsl_complex_float_rect(0e0, .1e0, $cα);
  mgsl_complex_float_rect(0e0, 0e0, $cβ);
  my gsl_matrix_complex_float $cA7 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(-.248e0, -.037e0, $c);
  mgsl_matrix_complex_float_set($cA7, 0, 0, $c);
  mgsl_complex_float_rect(-.124e0, .998e0, $c);
  mgsl_matrix_complex_float_set($cA7, 0, 1, $c);
  my gsl_matrix_complex_float $cB7 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(-.608e0, -.115e0, $c);
  mgsl_matrix_complex_float_set($cB7, 0, 0, $c);
  mgsl_complex_float_rect(-.718e0, -.551e0, $c);
  mgsl_matrix_complex_float_set($cB7, 0, 1, $c);
  my gsl_matrix_complex_float $cC7 = gsl_matrix_complex_float_calloc(1, 1);
  mgsl_complex_float_rect(.187e0, -.329e0, $c);
  mgsl_matrix_complex_float_set($cC7, 0, 0, $c);
  ok mgsl_blas_csyr2k(CblasUpper, CblasNoTrans, $cα, $cA7, $cB7, $cβ, $cC7) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC: complex32';
  is-deeply (gather { mgsl_matrix_complex_float_get($cC7, 0, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (0.11944480240345001e0, 0.15709181129932404e0), 'result: ok';

  mgsl_complex_rect(0e0, 0e0, $zα);
  mgsl_complex_rect(-.3e0, .1e0, $zβ);
  my gsl_matrix_complex $zA7 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(-.315e0, .03e0, $z);
  mgsl_matrix_complex_set($zA7, 0, 0, $z);
  mgsl_complex_rect(.218e0, .175e0, $z);
  mgsl_matrix_complex_set($zA7, 0, 1, $z);
  my gsl_matrix_complex $zB7 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(-.832e0, -.964e0, $z);
  mgsl_matrix_complex_set($zB7, 0, 0, $z);
  mgsl_complex_rect(.291e0, .476e0, $z);
  mgsl_matrix_complex_set($zB7, 0, 1, $z);
  my gsl_matrix_complex $zC7 = gsl_matrix_complex_calloc(1, 1);
  mgsl_complex_rect(-.341e0, .743e0, $z);
  mgsl_matrix_complex_set($zC7, 0, 0, $z);
  ok mgsl_blas_zsyr2k(CblasUpper, CblasNoTrans, $zα, $zA7, $zB7, $zβ, $zC7) == GSL_SUCCESS, 'αAT(B) + αBT(A) + βC: complex64';
  is-deeply (gather { mgsl_matrix_complex_get($zC7, 0, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (0.027999999999999997e0, -0.257e0), 'result: ok';

  mgsl_complex_float_rect(-1e0, 0e0, $cα);
  $β = -.3e0;
  my gsl_matrix_complex_float $cA8 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(.178e0, .545e0, $c);
  mgsl_matrix_complex_float_set($cA8, 0, 0, $c);
  mgsl_complex_float_rect(-.491e0, .979e0, $c);
  mgsl_matrix_complex_float_set($cA8, 0, 1, $c);
  my gsl_matrix_complex_float $cB8 = gsl_matrix_complex_float_calloc(1, 2);
  mgsl_complex_float_rect(-.665e0, -.531e0, $c);
  mgsl_matrix_complex_float_set($cB8, 0, 0, $c);
  mgsl_complex_float_rect(-.4e0, .227e0, $c);
  mgsl_matrix_complex_float_set($cB8, 0, 1, $c);
  my gsl_matrix_complex_float $cC8 = gsl_matrix_complex_float_calloc(1, 1);
  mgsl_complex_float_rect(.115e0, -.193e0, $c);
  mgsl_matrix_complex_float_set($cC8, 0, 0, $c);
  ok mgsl_blas_cher2k(CblasUpper, CblasNoTrans, $cα, $cA8, $cB8, $β, $cC8) == GSL_SUCCESS, 'αAH(B) + αBH(A) + βC: complex32';
  is-deeply (gather { mgsl_matrix_complex_float_get($cC8, 0, 0, $cres); take $cres.dat[0]; take $cres.dat[1] }), (-0.05623596906661987e0, 0e0), 'result: ok';

  mgsl_complex_rect(1e0, 0e0, $zα);
  $β = 1e0;
  my gsl_matrix_complex $zA8 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(.515e0, -.034e0, $z);
  mgsl_matrix_complex_set($zA8, 0, 0, $z);
  mgsl_complex_rect(.067e0, .66e0, $z);
  mgsl_matrix_complex_set($zA8, 0, 1, $z);
  my gsl_matrix_complex $zB8 = gsl_matrix_complex_calloc(1, 2);
  mgsl_complex_rect(.408e0, -.85e0, $z);
  mgsl_matrix_complex_set($zB8, 0, 0, $z);
  mgsl_complex_rect(-.945e0, -.799e0, $z);
  mgsl_matrix_complex_set($zB8, 0, 1, $z);
  my gsl_matrix_complex $zC8 = gsl_matrix_complex_calloc(1, 1);
  mgsl_complex_rect(-.918e0, -.985e0, $z);
  mgsl_matrix_complex_set($zC8, 0, 0, $z);
  ok mgsl_blas_zher2k(CblasUpper, CblasNoTrans, $zα, $zA8, $zB8, $β, $zC8) == GSL_SUCCESS, 'αAH(B) + αBH(A) + βC: complex64';
  is-deeply (gather { mgsl_matrix_complex_get($zC8, 0, 0, $zres); take $zres.dat[0]; take $zres.dat[1] }), (-1.62127e0, 0e0), 'result: ok';
}, 'level 3 functions';

done-testing;