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

use Test;
use lib 'lib';
use Math::Libgsl::Raw::LinearAlgebra :ALL;
use Math::Libgsl::Raw::Complex :ALL;
use Math::Libgsl::Raw::Matrix :ALL;
use Math::Libgsl::Raw::Matrix::Complex64 :ALL;
use Math::Libgsl::Raw::Permutation :ALL;
use Math::Libgsl::Raw::BLAS :ALL;
use Math::Libgsl::Raw::Random;
use Math::Libgsl::Constants;
use NativeCall;

sub create-random-matrix(gsl_matrix $m, gsl_rng $r)
{
  for ^$m.size1 X ^$m.size2 -> ($i, $j) {
    gsl_matrix_set($m, $i, $j, gsl_rng_uniform_pos($r));
  }
  return GSL_SUCCESS;
}

sub create-random-complex-matrix(gsl_matrix_complex $m, gsl_rng $r)
{
  my gsl_complex $mij = alloc_gsl_complex;
  for ^$m.size1 X ^$m.size2 -> ($i, $j) {
    mgsl_complex_rect(gsl_rng_uniform($r), gsl_rng_uniform($r), $mij);
    mgsl_matrix_complex_set($m, $i, $j, $mij);
  }
  free_gsl_complex($mij);
  return GSL_SUCCESS;
}

sub create-random-vector(gsl_vector $v, gsl_rng $r)
{
  for ^$v.size -> $i {
    my num64 $vi = gsl_rng_uniform($r);
    gsl_vector_set($v, $i, $vi);
  }
  return GSL_SUCCESS;
}

sub create-random-complex-vector(gsl_vector_complex $v, gsl_rng $r)
{
  my gsl_complex $vi = alloc_gsl_complex;
  for ^$v.size -> $i {
    mgsl_complex_rect(gsl_rng_uniform($r), gsl_rng_uniform($r), $vi);
    mgsl_vector_complex_set($v, $i, $vi);
  }
  free_gsl_complex($vi);
  return GSL_SUCCESS;
}

sub create-general-matrix(Int $size1, Int $size2 --> gsl_matrix)
{
  my gsl_matrix $m = gsl_matrix_calloc($size1, $size2);
  for ^$size1 X ^$size2 -> ($i, $j) {
    gsl_matrix_set($m, $i, $j, 1e0 / ($i + $j + 1e0));
  }
  return $m;
}

sub create-hilbert-matrix(Int $size)
{
  my gsl_matrix $m = gsl_matrix_calloc($size, $size);
  for ^$size X ^$size -> ($i, $j) {
    gsl_matrix_set($m, $i, $j, 1e0/($i + $j + 1e0));
  }
  return $m;
}

sub create-rank-matrix(size_t $rank, gsl_matrix $m, gsl_rng $r)
{
  my size_t $M = $m.size1;
  my size_t $N = $m.size2;
  my gsl_vector $u = gsl_vector_calloc($M);
  my gsl_vector $v = gsl_vector_calloc($N);
  gsl_matrix_set_zero($m);

  for ^$rank -> $i {
    create-random-vector($u, $r);
    create-random-vector($v, $r);
    gsl_blas_dger(1e0, $u, $v, $m);
  }

  gsl_vector_free($u);
  gsl_vector_free($v);
  return GSL_SUCCESS;
}

sub create-symm-matrix(gsl_matrix $m)
{
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  for ^$m.size1 -> $i {
    for 0..$i -> $j {
      gsl_matrix_set($m, $i, $j, gsl_rng_uniform($r));
    }
  }
  if $gsl-version > v2.5 {
    gsl_matrix_transpose_tricpy(CblasLower, CblasNonUnit, $m, $m);
  } else {
    gsl_matrix_transpose_tricpy('L'.ord, 0, $m, $m);
  }
  gsl_rng_free($r);
  return GSL_SUCCESS;
}

sub create-posdef-band-matrix(gsl_matrix $m)
{
  create-symm-matrix($m);
  my $size = $m.size1;
  my $α = $size * 10e0;
  gsl_matrix_set($m, $_, $_, gsl_matrix_get($m, $_, $_) + $α) for ^$size;
}

sub test-num64-matrix(gsl_matrix $m1, gsl_matrix $m2 where $m1.size1 == $m2.size1 && $m1.size2 == $m2.size2 --> Bool)
{
  my Bool $result = True;
  my $*TOLERANCE = 1e-12;
  for ^$m1.size1 X ^$m1.size2 -> ($i, $j) {
    $result &&= gsl_matrix_get($m1, $i, $j) ≅ gsl_matrix_get($m2, $i, $j);
  }
  $result
}

sub test-complex64-matrix(gsl_matrix_complex $m1, gsl_matrix_complex $m2 where $m1.size1 == $m2.size1 && $m1.size2 == $m2.size2 --> Bool)
{
  my Bool $result = True;
  my $*TOLERANCE = 1e-13;
  my gsl_complex $zres1 = alloc_gsl_complex;
  my gsl_complex $zres2 = alloc_gsl_complex;
  for ^$m1.size1 X ^$m1.size2 -> ($i, $j) {
    mgsl_matrix_complex_get($m1, $i, $j, $zres1);
    mgsl_matrix_complex_get($m2, $i, $j, $zres2);
    $result &&= $zres1.dat[0] ≅ $zres2.dat[0];
    $result &&= $zres1.dat[1] ≅ $zres2.dat[1];
  }
  free_gsl_complex($zres1);
  free_gsl_complex($zres2);
  $result
}

sub test-num64-vector(gsl_vector $v1, gsl_vector $v2 where $v1.size == $v2.size --> Bool)
{
  my $*TOLERANCE = 1e-10;
  #my Bool $result = gsl_vector_get($v1, $_) ≅ gsl_vector_get($v2, $_) for ^$v1.size;
  my Bool $result = True;
  for ^$v1.size -> $i {
    $result &&= gsl_vector_get($v1, $i) ≅ gsl_vector_get($v2, $i);
  }
  $result
}

sub test-complex64-vector(gsl_vector_complex $v1, gsl_vector_complex $v2 where $v1.size == $v2.size --> Bool)
{
  my Bool $result = True;
  my $*TOLERANCE = 1e-13;
  my gsl_complex $zres1 = alloc_gsl_complex;
  my gsl_complex $zres2 = alloc_gsl_complex;
  for ^$v1.size -> $i {
    mgsl_vector_complex_get($v1, $i, $zres1);
    mgsl_vector_complex_get($v2, $i, $zres2);
    $result &&= $zres1.dat[0] ≅ $zres2.dat[0];
    $result &&= $zres1.dat[1] ≅ $zres2.dat[1]
  }
  free_gsl_complex($zres1);
  free_gsl_complex($zres2);
  $result
}

sub LU-decomp-num64(gsl_matrix $m --> Bool)
{
  my int32 $signum;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix $A = gsl_matrix_calloc(10, 10);
  my gsl_matrix $PLU = gsl_matrix_calloc(10, 10);
  gsl_matrix_memcpy($A, $m);
  gsl_linalg_LU_decomp($A, $p, $signum);

  my gsl_vector_view $vv1 = alloc_gsl_vector_view();
  my gsl_vector_view $vv2 = alloc_gsl_vector_view();
  for ^10 -> $j {
    my gsl_vector $v1 = mgsl_matrix_subcolumn($vv1, $A, $j, $j, 10 - $j);
    my gsl_vector $v2 = mgsl_matrix_subcolumn($vv2, $PLU, $j, $j, 10 - $j);
    gsl_vector_memcpy($v2, $v1);
    gsl_matrix_set($PLU, $j, $j, 1e0);
  }
  free_gsl_vector_view($vv1);
  free_gsl_vector_view($vv2);
  gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1e0, $A, $PLU);
  my gsl_vector_view $vv = alloc_gsl_vector_view();
  for ^10 -> $j {
    my gsl_vector $v = mgsl_matrix_column($vv, $PLU, $j);
    gsl_permute_vector_inverse($p, $v);
  }
  free_gsl_vector_view($vv);
  my Bool $res = test-num64-matrix($m, $PLU);
  gsl_matrix_free($A);
  gsl_matrix_free($PLU);
  gsl_permutation_free($p);
  return $res;
}

sub LU-solve-num64(gsl_matrix $m, gsl_rng $r --> Array[Bool])
{
  my int32 $signum;
  my Bool @res;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_vector $rhs = gsl_vector_calloc(10);
  my gsl_vector $sol = gsl_vector_calloc(10);
  @res.push: create-random-vector($sol, $r) == GSL_SUCCESS;
  gsl_blas_dgemv(CblasNoTrans, 1e0, $m, $sol, 0e0, $rhs);
  my gsl_matrix $lu  = gsl_matrix_calloc(10, 10);
  my gsl_vector $x = gsl_vector_calloc(10);
  my gsl_vector $residual = gsl_vector_calloc(10);
  gsl_matrix_memcpy($lu, $m);

  gsl_linalg_LU_decomp($lu, $p, $signum);
  @res.push: gsl_linalg_LU_solve($lu, $p, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);
  @res.push: gsl_linalg_LU_refine($m, $lu, $p, $rhs, $x, $residual) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_matrix_free($lu);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_vector_free($residual);
  gsl_vector_free($x);
  gsl_permutation_free($p);
  return @res;
}

sub LU-solve-complex64(gsl_matrix_complex $cm, gsl_rng $r --> Array[Bool])
{
  my int32 $signum;
  my Bool @res;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_vector_complex $crhs = gsl_vector_complex_alloc(10);
  my gsl_vector_complex $csol = gsl_vector_complex_alloc(10);
  @res.push: create-random-complex-vector($csol, $r) == GSL_SUCCESS;
  my gsl_complex $zero = alloc_gsl_complex;
  mgsl_complex_rect(0e0, 0e0, $zero);
  my gsl_complex $one = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 0e0, $one);
  mgsl_blas_zgemv(CblasNoTrans, $one, $cm, $csol, $zero, $crhs);
  my gsl_matrix_complex $clu  = gsl_matrix_complex_alloc(10, 10);
  my gsl_vector_complex $cx = gsl_vector_complex_alloc(10);
  my gsl_vector_complex $cresidual = gsl_vector_complex_alloc(10);
  gsl_matrix_complex_memcpy($clu, $cm);
  gsl_linalg_complex_LU_decomp($clu, $p, $signum);
  @res.push: gsl_linalg_complex_LU_solve($clu, $p, $crhs, $cx) == GSL_SUCCESS;
  @res.push: test-complex64-vector($cx, $csol);
  @res.push: gsl_linalg_complex_LU_refine($cm, $clu, $p, $crhs, $cx, $cresidual) == GSL_SUCCESS;
  @res.push: test-complex64-vector($cx, $csol);

  gsl_matrix_complex_free($clu);
  gsl_vector_complex_free($crhs);
  gsl_vector_complex_free($csol);
  gsl_vector_complex_free($cresidual);
  gsl_vector_complex_free($cx);
  free_gsl_complex($zero);
  free_gsl_complex($one);
  gsl_permutation_free($p);
  return @res;
}

sub LU-invert-num64(gsl_matrix $m --> Array[Bool])
{
  my int32 $signum;
  my Bool @res;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix $lu  = gsl_matrix_calloc(10, 10);
  my gsl_matrix $inv = gsl_matrix_calloc(10, 10);
  my gsl_matrix $c   = gsl_matrix_calloc(10, 10);
  gsl_matrix_memcpy($lu, $m);
  gsl_linalg_LU_decomp($lu, $p, $signum);

  @res.push: gsl_linalg_LU_invert($lu, $p, $inv) == GSL_SUCCESS;

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $m, $inv, 0e0, $c);
  @res.push: True;
  for ^10 -> $i {
    for ^10 -> $j {
      @res[1] &&= gsl_matrix_get($c, $i, $j) ≅ ($i == $j ?? 1e0 !! 0e0);
    }
  }

  gsl_matrix_free($c);
  gsl_matrix_free($inv);
  gsl_matrix_free($lu);
  gsl_permutation_free($p);
  return @res;
}

sub LU-invert-complex64(gsl_matrix_complex $cm --> Array[Bool])
{
  my int32 $signum;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix_complex $clu  = gsl_matrix_complex_alloc(10, 10);
  my gsl_matrix_complex $cinv = gsl_matrix_complex_alloc(10, 10);
  my gsl_matrix_complex $cc   = gsl_matrix_complex_alloc(10, 10);
  my gsl_complex $zero = alloc_gsl_complex;
  mgsl_complex_rect(0e0, 0e0, $zero);
  my gsl_complex $one = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 0e0, $one);
  gsl_matrix_complex_memcpy($clu, $cm);
  gsl_linalg_complex_LU_decomp($clu, $p, $signum);

  my Bool @res;
  @res.push: gsl_linalg_complex_LU_invert($clu, $p, $cinv) == GSL_SUCCESS;

  mgsl_blas_zgemm(CblasNoTrans, CblasNoTrans, $one, $cm, $cinv, $zero, $cc);
  @res.push: True;
  my $*TOLERANCE = 1e-13;
  my gsl_complex $cij = alloc_gsl_complex;
  for ^10 -> $i {
    for ^10 -> $j {
      mgsl_matrix_complex_get($cc, $i, $j, $cij);
      @res[1] &&= (($cij.dat[0] // 0e0) ≅ ($i == $j ?? 1e0 !! 0e0) && ($cij.dat[1] // 0e0) ≅ 0e0);
    }
  }
  free_gsl_complex($cij);

  gsl_matrix_complex_free($cc);
  gsl_matrix_complex_free($cinv);
  gsl_matrix_complex_free($clu);
  free_gsl_complex($zero);
  free_gsl_complex($one);
  gsl_permutation_free($p);
  return @res;
}

sub LU-det-num64(gsl_matrix $m --> Num)
{
  my int32 $signum;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix $lu  = gsl_matrix_calloc(10, 10);
  gsl_matrix_memcpy($lu, $m);
  gsl_linalg_LU_decomp($lu, $p, $signum);
  my $ret = gsl_linalg_LU_det($lu, $signum);
  gsl_matrix_free($lu);
  gsl_permutation_free($p);
  $ret;
}

sub LU-det-complex64(gsl_matrix_complex $cm, gsl_complex $res)
{
  my int32 $signum;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix_complex $clu  = gsl_matrix_complex_alloc(10, 10);
  gsl_matrix_complex_memcpy($clu, $cm);
  gsl_linalg_complex_LU_decomp($clu, $p, $signum);
  mgsl_linalg_complex_LU_det($clu, $signum, $res);
  gsl_matrix_complex_free($clu);
  gsl_permutation_free($p);
  $res;
}

sub LU-sgndet-num64(gsl_matrix $m --> Int)
{
  my int32 $signum;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix $lu  = gsl_matrix_calloc(10, 10);
  gsl_matrix_memcpy($lu, $m);
  gsl_linalg_LU_decomp($lu, $p, $signum);
  my $ret = gsl_linalg_LU_sgndet($lu, $signum);
  gsl_matrix_free($lu);
  gsl_permutation_free($p);
  $ret;
}

sub LU-sgndet-complex64(gsl_matrix_complex $cm --> Complex)
{
  my int32 $signum;
  my gsl_permutation $p = gsl_permutation_alloc(10);
  my gsl_matrix_complex $clu  = gsl_matrix_complex_alloc(10, 10);
  gsl_matrix_complex_memcpy($clu, $cm);
  gsl_linalg_complex_LU_decomp($clu, $p, $signum);
  my gsl_complex $ret = alloc_gsl_complex;
  mgsl_linalg_complex_LU_sgndet($clu, $signum, $ret);
  gsl_matrix_complex_free($clu);
  my Complex $c = $ret.dat[0] + $ret.dat[1] * i;
  free_gsl_complex($ret);
  gsl_permutation_free($p);
  $c;
}

sub QR-decomp(gsl_matrix $m)
{
  my $M = $m.size1;
  my $N = $m.size2;
  my gsl_matrix $qr = gsl_matrix_calloc($M, $N);
  my gsl_matrix $a  = gsl_matrix_calloc($M, $N);
  my gsl_matrix $q  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $r  = gsl_matrix_calloc($M, $N);
  my gsl_vector $d  = gsl_vector_calloc(min($M, $N));
  my Bool @res;
  gsl_matrix_memcpy($qr, $m);

  @res.push: gsl_linalg_QR_decomp($qr, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_unpack($qr, $d, $q, $r) == GSL_SUCCESS;

  @res.push: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $q, $r, 0e0, $a) == GSL_SUCCESS;
  @res.push: test-num64-matrix($a, $m);

  gsl_vector_free($d);
  gsl_matrix_free($qr);
  gsl_matrix_free($a);
  gsl_matrix_free($q);
  gsl_matrix_free($r);
  return @res;
}

sub QR-solve(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_vector $rhs = gsl_vector_calloc($dim);
  my gsl_matrix $qr  = gsl_matrix_calloc($dim, $dim);
  my gsl_vector $d   = gsl_vector_calloc($dim);
  my gsl_vector $x   = gsl_vector_calloc($dim);
  my gsl_vector $sol = gsl_vector_calloc($dim);
  my Bool @res;
  gsl_matrix_memcpy($qr, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$dim;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_QR_decomp($qr, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_solve($qr, $d, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($sol);
  gsl_vector_free($x);
  gsl_vector_free($d);
  gsl_matrix_free($qr);
  gsl_vector_free($rhs);
  return @res;
}

sub QR-lssolve(gsl_matrix $m)
{
  my $M = $m.size1;
  my $N = $m.size2;
  my gsl_vector $rhs = gsl_vector_calloc($M);
  my gsl_matrix $qr  = gsl_matrix_calloc($M, $N);
  my gsl_vector $d   = gsl_vector_calloc($N);
  my gsl_vector $x   = gsl_vector_calloc($N);
  my gsl_vector $r   = gsl_vector_calloc($M);
  my gsl_vector $res = gsl_vector_calloc($M);
  my gsl_vector $sol = gsl_vector_calloc($N);
  my Bool @res;
  gsl_matrix_memcpy($qr, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$M;
  gsl_vector_set($sol, 0, 52.5992295702070e0);
  gsl_vector_set($sol, 1, -337.7263113752073e0);
  gsl_vector_set($sol, 2, 351.8823436427604e0);

  @res.push: gsl_linalg_QR_decomp($qr, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_lssolve($qr, $d, $rhs, $x, $res) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($sol);
  gsl_vector_free($res);
  gsl_vector_free($r);
  gsl_vector_free($x);
  gsl_vector_free($d);
  gsl_matrix_free($qr);
  gsl_vector_free($rhs);
  return @res;
}

sub QR-QTmat(gsl_matrix $A, gsl_matrix $B)
{
  my $M = $A.size1;
  my $N = $A.size2;
  my $K = $B.size2;
  my gsl_matrix $QR   = gsl_matrix_calloc($M, $N);
  my gsl_matrix $work = gsl_matrix_calloc($N, $K);
  my gsl_matrix $B1   = gsl_matrix_calloc($M, $K);
  my gsl_vector $tau  = gsl_vector_calloc($N);
  my @sol[$M; $K] = [
    [-1.208620986625094, -1.0546472240861027, -0.7791902451463131, -1.168773675878457, -1.8023918740364615,],
    [-0.7255813094089609, -0.6599481245008292, -0.6132935912751079, -0.6001656202344267, -0.9656995184226206,],
    [-0.2706435681806497, -0.37683693893675707, 0.1454012389605544, -0.2524856382131687, -0.1832861493887159,],
    [-0.25346494022471294, -0.034976909820802315, -0.40523345579661035, -0.3042534812505212, -0.7456694840765314,],
    [0.3232602101781864, 0.2187167384986574, 0.39770492818806313, 0.18484732428275386, 0.6840232171558648,],
    [-0.12627206077967576, 0.22810152457644187, -0.08358600426058016, 0.3855951837101906, 0.18596814875099849,],
    [-0.5648471145968404, -0.12230132226051842, -0.6255138866668446, -0.6991112162503683, -0.22180105879372478,],
    [-0.2514851779890921, 0.26868394934112794, 0.11129620314236724, 0.7323300527937038, -0.2037361807555147,],
    [0.5294793846838799, 0.12853319270334268, 0.255770526547394, 0.34313072822437585, -0.25510668175740325,],
    [-0.30702000923201295, 0.4330172700274453, 0.4160212523423857, -0.1256545168104251, 0.025424204736825898,],
  ];

  my Bool @res;
  gsl_matrix_memcpy($QR, $A);
  gsl_matrix_memcpy($B1, $B);

  @res.push: gsl_linalg_QR_decomp($QR, $tau) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_QTmat($QR, $tau, $B1) == GSL_SUCCESS;

  my $*TOLERANCE = 1e-13;
  my Bool $pres = True;
  for ^$M -> $i {
    for ^$K -> $j {
      $pres &&= gsl_matrix_get($B1, $i, $j) ≅ @sol[$i; $j];
    }
  }
  @res.push: $pres;

  gsl_vector_free($tau);
  gsl_matrix_free($QR);
  gsl_matrix_free($B1);
  gsl_matrix_free($work);
  return @res;
}

sub QR-QRsolve(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_vector $rhs = gsl_vector_calloc($dim);
  my gsl_matrix $qr  = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $q   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $r   = gsl_matrix_calloc($dim, $dim);
  my gsl_vector $d   = gsl_vector_calloc($dim);
  my gsl_vector $x   = gsl_vector_calloc($dim);
  my gsl_vector $sol = gsl_vector_calloc($dim);
  my Bool @res;
  gsl_matrix_memcpy($qr, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$dim;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_QR_decomp($qr, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_unpack($qr, $d, $q, $r) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_QRsolve($q, $r, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($sol);
  gsl_vector_free($x);
  gsl_vector_free($d);
  gsl_matrix_free($qr);
  gsl_matrix_free($q);
  gsl_matrix_free($r);
  gsl_vector_free($rhs);
  return @res;
}

sub QR-update(gsl_matrix $m)
{
  my $M = $m.size1;
  my $N = $m.size2;
  my gsl_matrix $qr1  = gsl_matrix_calloc($M, $N);
  my gsl_matrix $qr2  = gsl_matrix_calloc($M, $N);
  my gsl_matrix $q1   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $r1   = gsl_matrix_calloc($M, $N);
  my gsl_matrix $q2   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $r2   = gsl_matrix_calloc($M, $N);
  my gsl_vector $d    = gsl_vector_calloc(min($M, $N));
  my gsl_vector $u    = gsl_vector_calloc($M);
  my gsl_vector $v    = gsl_vector_calloc($N);
  my gsl_vector $w    = gsl_vector_calloc($M);
  my Bool @res;
  gsl_matrix_memcpy($qr1, $m);
  gsl_matrix_memcpy($qr2, $m);
  gsl_vector_set($u, $_, sin($_ + 1e0)) for ^$M;
  gsl_vector_set($v, $_, cos($_ + 2e0) + sin($_ * $_ + 3e0)) for ^$N;

  for ^$M -> $i {
    my num64 $ui = gsl_vector_get($u, $i);
    for ^$N -> $j {
      my num64 $vj  = gsl_vector_get($v, $j);
      my num64 $qij = gsl_matrix_get($qr1, $i, $j);
      gsl_matrix_set($qr1, $i, $j, $qij + $ui * $vj);
    }
  }

  @res.push: gsl_linalg_QR_decomp($qr2, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_unpack($qr2, $d, $q2, $r2) == GSL_SUCCESS;

  for ^$M -> $j {
    my num64 $sum = 0e0;
    for ^$M -> $i {
      $sum += gsl_matrix_get($q2, $i, $j) * gsl_vector_get($u, $i);
    }
    gsl_vector_set($w, $j, $sum);
  }

  @res.push: gsl_linalg_QR_update($q2, $r2, $w, $v) == GSL_SUCCESS;

  for ^$M -> $i {
    for ^$N -> $j {
      my num64 $sum = 0e0;
      for 0..min($j, $M - 1) -> $k {
        my num64 $qik = gsl_matrix_get($q2, $i, $k);
        my num64 $rkj = gsl_matrix_get($r2, $k, $j);
        $sum += $qik * $rkj;
      }
      gsl_matrix_set($qr2, $i, $j, $sum);
    }
  }

  @res.push: test-num64-matrix($qr1, $qr2);

  gsl_vector_free($d);
  gsl_vector_free($u);
  gsl_vector_free($v);
  gsl_vector_free($w);
  gsl_matrix_free($qr1);
  gsl_matrix_free($qr2);
  gsl_matrix_free($q1);
  gsl_matrix_free($r1);
  gsl_matrix_free($q2);
  gsl_matrix_free($r2);
  return @res;
}

sub QRPT-decomp(gsl_matrix $m)
{
  my $M = $m.size1;
  my $N = $m.size2;
  my gsl_matrix $QR   = gsl_matrix_calloc($M, $N);
  my gsl_matrix $A    = gsl_matrix_calloc($M, $N);
  my gsl_matrix $Q    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $R    = gsl_matrix_calloc($M, $N);
  my gsl_vector $tau  = gsl_vector_calloc(min($M, $N));
  my gsl_vector $norm = gsl_vector_calloc($N);
  my gsl_permutation $p = gsl_permutation_alloc($N);
  my int32 $signum;
  my Bool @res;
  gsl_matrix_memcpy($QR, $m);

  @res.push: gsl_linalg_QRPT_decomp($QR, $tau, $p, $signum, $norm) == GSL_SUCCESS;
  @res.push: gsl_linalg_QR_unpack($QR, $tau, $Q, $R) == GSL_SUCCESS;
  @res.push: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $Q, $R, 0e0, $A) == GSL_SUCCESS;

  my gsl_matrix_view $mv = alloc_gsl_matrix_view();
  for ^$M -> $i {
    my gsl_vector $vv = mgsl_matrix_row($mv, $A, $i);
    gsl_permute_vector_inverse($p, $vv);
  }
  free_gsl_matrix_view($mv);

  @res.push: test-num64-matrix($A, $m);

 