### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 2 of 4
 my gsl_vector $work = gsl_vector_calloc(3 * $N);
  my num64 $rcond;
  gsl_linalg_QRPT_rcond($QR, $rcond, $work);
  my $*TOLERANCE = 1e-10;
  @res.push: $rcond ≅ 2.915362697820e-03;

  gsl_vector_free($work);
  gsl_permutation_free($p);
  gsl_vector_free($norm);
  gsl_vector_free($tau);
  gsl_matrix_free($QR);
  gsl_matrix_free($A);
  gsl_matrix_free($Q);
  gsl_matrix_free($R);
  return @res;
}

sub QRPT-decomp2(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $QR   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $Q    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $R    = gsl_matrix_calloc($M, $M);
  my gsl_vector $rhs  = gsl_vector_calloc($M);
  my gsl_vector $d    = gsl_vector_calloc($M);
  my gsl_vector $x    = gsl_vector_calloc($M);
  my gsl_vector $norm = gsl_vector_calloc($M);
  my gsl_vector $sol  = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my int32 $signum;
  my Bool @res;
  gsl_matrix_memcpy($QR, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$M;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_QRPT_decomp2($QR, $Q, $R, $d, $p, $signum, $norm) == GSL_SUCCESS;
  @res.push: gsl_linalg_QRPT_QRsolve($Q, $R, $p, $rhs, $x) == GSL_SUCCESS;

  @res.push: test-num64-vector($x, $sol);

  gsl_permutation_free($p);
  gsl_vector_free($norm);
  gsl_vector_free($x);
  gsl_vector_free($d);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_matrix_free($QR);
  gsl_matrix_free($Q);
  gsl_matrix_free($R);
  return @res;
}

sub QRPT-solve(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_vector $sol  = gsl_vector_calloc($M);
  my gsl_vector $rhs  = gsl_vector_calloc($M);
  my gsl_vector $d    = gsl_vector_calloc($M);
  my gsl_vector $x    = gsl_vector_calloc($M);
  my gsl_vector $norm = gsl_vector_calloc($M);
  my gsl_matrix $QR   = gsl_matrix_calloc($M, $M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my int32 $signum;
  my Bool @res;
  gsl_matrix_memcpy($QR, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$M;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_QRPT_decomp($QR, $d, $p, $signum, $norm) == GSL_SUCCESS;
  @res.push: gsl_linalg_QRPT_solve($QR, $d, $p, $rhs, $x) == GSL_SUCCESS;

  @res.push: test-num64-vector($x, $sol);

  gsl_permutation_free($p);
  gsl_vector_free($norm);
  gsl_vector_free($x);
  gsl_vector_free($d);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_matrix_free($QR);
  return @res;
}

sub QRPT-lssolve(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_vector $sol  = gsl_vector_calloc($M);
  my gsl_vector $rhs  = gsl_vector_calloc($M);
  my gsl_vector $tau  = gsl_vector_calloc($M);
  my gsl_vector $x    = gsl_vector_calloc($M);
  my gsl_vector $r    = gsl_vector_calloc($M);
  my gsl_vector $res  = gsl_vector_calloc($M);
  my gsl_vector $norm = gsl_vector_calloc($M);
  my gsl_vector $work = gsl_vector_calloc($M);
  my gsl_matrix $QRPT = gsl_matrix_calloc($M, $M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my int32 $signum;
  my Bool @res;
  gsl_matrix_memcpy($QRPT, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$M;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_QRPT_decomp($QRPT, $tau, $p, $signum, $work) == GSL_SUCCESS;
  @res.push: gsl_linalg_QRPT_lssolve($QRPT, $tau, $p, $rhs, $x, $res) == GSL_SUCCESS;

  @res.push: test-num64-vector($x, $sol);

  gsl_vector_set_zero($r);

  @res.push: test-num64-vector($r, $res);

  gsl_permutation_free($p);
  gsl_vector_free($work);
  gsl_vector_free($norm);
  gsl_vector_free($r);
  gsl_vector_free($x);
  gsl_vector_free($tau);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_matrix_free($QRPT);
  return @res;
}

sub QRPT-lssolve2(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_vector $sol  = gsl_vector_calloc($M);
  my gsl_vector $rhs  = gsl_vector_calloc($M);
  my gsl_vector $tau  = gsl_vector_calloc($M);
  my gsl_vector $x    = gsl_vector_calloc($M);
  my gsl_vector $r    = gsl_vector_calloc($M);
  my gsl_vector $res  = gsl_vector_calloc($M);
  my gsl_vector $work = gsl_vector_calloc($M);
  my gsl_matrix $QRPT = gsl_matrix_calloc($M, $M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my int32 $signum;
  my Bool @res;
  gsl_matrix_memcpy($QRPT, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$M;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_QRPT_decomp($QRPT, $tau, $p, $signum, $work) == GSL_SUCCESS;
  my size_t $rank = gsl_linalg_QRPT_rank($QRPT, -1e0);
  @res.push: gsl_linalg_QRPT_lssolve2($QRPT, $tau, $p, $rhs, $rank, $x, $res) == GSL_SUCCESS;

  @res.push: test-num64-vector($x, $sol);

  gsl_vector_set_zero($r);

  @res.push: test-num64-vector($r, $res);

  gsl_permutation_free($p);
  gsl_vector_free($work);
  gsl_vector_free($r);
  gsl_vector_free($x);
  gsl_vector_free($tau);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_matrix_free($QRPT);
  return @res;
}

sub QRPT-update(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $QR1  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $QR2  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $Q1   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $Q2   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $R1   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $R2   = gsl_matrix_calloc($M, $M);
  my gsl_vector $d    = gsl_vector_calloc($M);
  my gsl_vector $u    = gsl_vector_calloc($M);
  my gsl_vector $v    = gsl_vector_calloc($M);
  my gsl_vector $w    = gsl_vector_calloc($M);
  my gsl_vector $norm = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my int32 $signum;
  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 ^$M;
  for ^$M -> $i {
    my num64 $ui = gsl_vector_get($u, $i);
    for ^$M -> $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_QRPT_decomp($QR2, $d, $p, $signum, $norm) == 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_QRPT_update($Q2, $R2, $p, $w, $v) == GSL_SUCCESS;

  for ^$M -> $i {
    for ^$M -> $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);
    }
  }

  my gsl_vector_view $r_i = alloc_gsl_vector_view();
  for ^$M -> $i {
    my gsl_vector $v = mgsl_matrix_row($r_i, $QR2, $i);
    gsl_permute_vector_inverse($p, $v);
  }
  free_gsl_vector_view($r_i);

  @res.push: test-num64-matrix($QR1, $QR2);

  gsl_permutation_free($p);
  gsl_vector_free($norm);
  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($Q2);
  gsl_matrix_free($R1);
  gsl_matrix_free($R2);
  return @res;
}

sub LQ-solve(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_matrix $lq  = gsl_matrix_calloc($dim, $dim);
  my gsl_vector $rhs = gsl_vector_calloc($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);
  gsl_matrix_transpose_memcpy($lq, $m);
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$dim;
  my Bool @res;

  @res.push: gsl_linalg_LQ_decomp($lq, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_LQ_solve_T($lq, $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($lq);
  gsl_vector_free($rhs);
  return @res;
}

sub LQ-LQsolve(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_matrix $lq  = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $q   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $l   = gsl_matrix_calloc($dim, $dim);
  my gsl_vector $rhs = gsl_vector_calloc($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);
  gsl_matrix_transpose_memcpy($lq, $m);
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$dim;
  my Bool @res;

  @res.push: gsl_linalg_LQ_decomp($lq, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_LQ_unpack($lq, $d, $q, $l) == GSL_SUCCESS;
  @res.push: gsl_linalg_LQ_LQsolve($q, $l, $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($lq);
  gsl_matrix_free($q);
  gsl_matrix_free($l);
  gsl_vector_free($rhs);
  return @res;
}

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

  @res.push: gsl_linalg_LQ_decomp($lq, $d) == GSL_SUCCESS;
  @res.push: gsl_linalg_LQ_lssolve_T($lq, $d, $rhs, $x, $res) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_set_zero($r);
  @res.push: test-num64-vector($res, $r);

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

sub LQ-decomp(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_matrix $lq  = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $a   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $q   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $l   = gsl_matrix_calloc($dim, $dim);
  my gsl_vector $d   = gsl_vector_calloc($dim);
  gsl_matrix_memcpy($lq, $m);
  my Bool @res;

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

  gsl_matrix_free($a);
  gsl_matrix_free($q);
  gsl_matrix_free($l);
  gsl_matrix_free($lq);
  gsl_vector_free($d);
  return @res;
}

sub COD-decomp(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $QRTZ   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $Q      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $R      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $QR     = gsl_matrix_calloc($M, $M);
  my gsl_matrix $Z      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $lhs    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $rhs    = gsl_matrix_calloc($M, $M);
  my gsl_vector $tau_Q  = gsl_vector_calloc($M);
  my gsl_vector $tau_Z  = gsl_vector_calloc($M);
  my gsl_vector $work   = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  gsl_matrix_memcpy($QRTZ, $m);
  my size_t $rank;
  my Bool @res;

  @res.push: gsl_linalg_COD_decomp($QRTZ, $tau_Q, $tau_Z, $p, $rank, $work) == GSL_SUCCESS;
  @res.push: gsl_linalg_COD_unpack($QRTZ, $tau_Q, $tau_Z, $rank, $Q, $R, $Z) == GSL_SUCCESS;

  gsl_matrix_memcpy($lhs, $m);
  gsl_permute_matrix($p, $lhs);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $Q, $R, 0e0, $QR);
  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1e0, $QR, $Z, 0e0, $rhs);

  @res.push: test-num64-matrix($rhs, $lhs);

  gsl_permutation_free($p);
  gsl_vector_free($work);
  gsl_vector_free($tau_Q);
  gsl_vector_free($tau_Z);
  gsl_matrix_free($QRTZ);
  gsl_matrix_free($lhs);
  gsl_matrix_free($rhs);
  gsl_matrix_free($QR);
  gsl_matrix_free($Q);
  gsl_matrix_free($R);
  gsl_matrix_free($Z);
  return @res;
}

sub COD-lssolve(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $QRTZ   = gsl_matrix_calloc($M, $M);
  my gsl_vector $lhs    = gsl_vector_calloc($M);
  my gsl_vector $rhs    = gsl_vector_calloc($M);
  my gsl_vector $tau_Q  = gsl_vector_calloc($M);
  my gsl_vector $tau_Z  = gsl_vector_calloc($M);
  my gsl_vector $work   = gsl_vector_calloc($M);
  my gsl_vector $x      = gsl_vector_calloc($M);
  my gsl_vector $r      = gsl_vector_calloc($M);
  my gsl_vector $res    = gsl_vector_calloc($M);
  my gsl_vector $sol    = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my size_t $rank;
  my Bool @res;
  gsl_matrix_memcpy($QRTZ, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$M;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_COD_decomp($QRTZ, $tau_Q, $tau_Z, $p, $rank, $work) == GSL_SUCCESS;
  @res.push: gsl_linalg_COD_lssolve($QRTZ, $tau_Q, $tau_Z, $p, $rank, $rhs, $x, $res) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_set_zero($r);

  @res.push: test-num64-vector($res, $r);

  gsl_permutation_free($p);
  gsl_vector_free($sol);
  gsl_vector_free($res);
  gsl_vector_free($r);
  gsl_vector_free($x);
  gsl_vector_free($work);
  gsl_vector_free($tau_Q);
  gsl_vector_free($tau_Z);
  gsl_matrix_free($QRTZ);
  gsl_vector_free($lhs);
  gsl_vector_free($rhs);
  return @res;
}

sub COD-lssolve2
{
  my $M = 100;
  my $N = 50;
  my num64 $lambda = 2.3e0;
  my gsl_rng $rng       = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $A      = gsl_matrix_calloc($M, $N);
  my gsl_vector $b      = gsl_vector_calloc($M);
  create-rank-matrix(26, $A, $rng);
  create-random-vector($b, $rng);
  my gsl_vector $lhs    = gsl_vector_calloc($M);
  my gsl_matrix $QRZT   = gsl_matrix_calloc($M, $N);
  my gsl_vector $tau_Q  = gsl_vector_calloc(min($M, $N));
  my gsl_vector $tau_Z  = gsl_vector_calloc(min($M, $N));
  my gsl_vector $work   = gsl_vector_calloc($N);
  my gsl_vector $x      = gsl_vector_calloc($N);
  my gsl_vector $x_aug  = gsl_vector_calloc($N);
  my gsl_vector $r      = gsl_vector_calloc($M);
  my gsl_vector $res    = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($N);
  my size_t $rank;
  my Bool @res;

  {
    my gsl_vector_view $v = alloc_gsl_vector_view();
    my gsl_matrix_view $m = alloc_gsl_matrix_view();
    my gsl_vector_view $d = alloc_gsl_vector_view();
    my gsl_permutation $p = gsl_permutation_alloc($N);
    my gsl_matrix $B = gsl_matrix_calloc($M + $N, $N);
    my gsl_vector $f = gsl_vector_calloc($M + $N);
    my gsl_vector $tau = gsl_vector_calloc($N);
    my gsl_vector $residual = gsl_vector_calloc($M + $N);
    my int32 $signum;

    my gsl_matrix $sm = mgsl_matrix_submatrix($m, $B, 0, 0, $M, $N);
    gsl_matrix_memcpy($sm, $A);
    $sm = mgsl_matrix_submatrix($m, $B, $M, 0, $N, $N);
    my gsl_vector $dv = mgsl_matrix_diagonal($d, $sm);
    gsl_vector_set_all($dv, $lambda);

    my gsl_vector $sv = mgsl_vector_subvector($v, $f, 0, $M);
    gsl_vector_memcpy($sv, $b);
    gsl_linalg_QRPT_decomp($B, $tau, $p, $signum, $work);
    gsl_linalg_QRPT_lssolve($B, $tau, $p, $f, $x_aug, $residual);

    free_gsl_vector_view($d);
    free_gsl_vector_view($v);
    free_gsl_matrix_view($m);
    gsl_permutation_free($p);
    gsl_matrix_free($B);
    gsl_vector_free($f);
    gsl_vector_free($tau);
    gsl_vector_free($residual);
  }

  gsl_matrix_memcpy($QRZT, $A);

  @res.push: gsl_linalg_COD_decomp($QRZT, $tau_Q, $tau_Z, $p, $rank, $work) == GSL_SUCCESS;

  {
    my gsl_matrix $S = gsl_matrix_calloc($rank, $rank);
    my gsl_vector $workr = gsl_vector_calloc($rank);

    @res.push: gsl_linalg_COD_lssolve2($lambda, $QRZT, $tau_Q, $tau_Z, $p, $rank, $b, $x, $res, $S, $workr) == GSL_SUCCESS;

    gsl_matrix_free($S);
    gsl_vector_free($workr);
  }

  @res.push: test-num64-vector($x, $x_aug);

  if $M == $N {
    gsl_vector_set_zero($r);
  } else {
    gsl_vector_memcpy($r, $b);
    gsl_blas_dgemv(CblasNoTrans, -1e0, $A, $x, 1e0, $r);
  }

  @res.push: test-num64-vector($r, $res);

  gsl_rng_free($rng);
  gsl_vector_free($r);
  gsl_vector_free($res);
  gsl_vector_free($x);
  gsl_vector_free($x_aug);
  gsl_vector_free($tau_Q);
  gsl_vector_free($tau_Z);
  gsl_matrix_free($QRZT);
  gsl_vector_free($lhs);
  gsl_vector_free($work);
  gsl_permutation_free($p);
  return @res;
}

sub SV-decomp(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $v   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $a   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $q   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $dqt = gsl_matrix_calloc($M, $M);
  my gsl_vector $d   = gsl_vector_calloc($M);
  my gsl_vector $w   = gsl_vector_calloc($M);
  my num64 $di1 = 0e0;
  my Bool @res;
  gsl_matrix_memcpy($v, $m);

  @res.push: gsl_linalg_SV_decomp($v, $q, $d, $w) == GSL_SUCCESS;

  for ^$M -> $i {
    my num64 $di = gsl_vector_get($d, $i);
    for ^$M -> $j {
      my num64 $qji = gsl_matrix_get($q, $j, $i);
      gsl_matrix_set($dqt, $i, $j, $qji * $di);
    }
  }

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $v, $dqt, 0e0, $a);

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

  gsl_vector_free($w);
  gsl_vector_free($d);
  gsl_matrix_free($v);
  gsl_matrix_free($a);
  gsl_matrix_free($q);
  gsl_matrix_free($dqt);
  return @res;
}

sub SV-decomp-jacobi(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $v   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $a   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $q   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $dqt = gsl_matrix_calloc($M, $M);
  my gsl_vector $d   = gsl_vector_calloc($M);
  my num64 $di1 = 0e0;
  my Bool @res;
  gsl_matrix_memcpy($v, $m);

  @res.push: gsl_linalg_SV_decomp_jacobi($v, $q, $d) == GSL_SUCCESS;

  for ^$M -> $i {
    my num64 $di = gsl_vector_get($d, $i);
    for ^$M -> $j {
      my num64 $qji = gsl_matrix_get($q, $j, $i);
      gsl_matrix_set($dqt, $i, $j, $qji * $di);
    }
  }

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $v, $dqt, 0e0, $a);

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

  gsl_vector_free($d);
  gsl_matrix_free($v);
  gsl_matrix_free($a);
  gsl_matrix_free($q);
  gsl_matrix_free($dqt);
  return @res;
}

sub SV-decomp-mod(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $v   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $a   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $q   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $x   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $dqt = gsl_matrix_calloc($M, $M);
  my gsl_vector $d   = gsl_vector_calloc($M);
  my gsl_vector $w   = gsl_vector_calloc($M);
  my num64 $di1 = 0e0;
  my Bool @res;
  gsl_matrix_memcpy($v, $m);

  @res.push: gsl_linalg_SV_decomp_mod($v, $x, $q, $d, $w) == GSL_SUCCESS;

  for ^$M -> $i {
    my num64 $di = gsl_vector_get($d, $i);
    for ^$M -> $j {
      my num64 $qji = gsl_matrix_get($q, $j, $i);
      gsl_matrix_set($dqt, $i, $j, $qji * $di);
    }
  }

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $v, $dqt, 0e0, $a);

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

  gsl_vector_free($w);
  gsl_vector_free($d);
  gsl_matrix_free($x);
  gsl_matrix_free($v);
  gsl_matrix_free($a);
  gsl_matrix_free($q);
  gsl_matrix_free($dqt);
  return @res;
}

sub SV-solve(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_vector $rhs = gsl_vector_calloc($dim);
  my gsl_matrix $u   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $q   = 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($u, $m);
  gsl_vector_set($rhs, $_, $_ + 1e0) for ^$dim;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_SV_decomp($u, $q, $d, $x) == GSL_SUCCESS;
  @res.push: gsl_linalg_SV_solve($u, $q, $d, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

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

sub Cholesky-decomp(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_matrix $V    = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $A    = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $L    = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $LT   = gsl_matrix_calloc($dim, $dim);
  my Bool @res;
  gsl_matrix_memcpy($V, $m);

  @res.push: gsl_linalg_cholesky_decomp1($V) == GSL_SUCCESS;

  if $gsl-version > v2.5 {
    gsl_matrix_tricpy(CblasLower, CblasNonUnit, $L, $V);
    gsl_matrix_transpose_tricpy(CblasLower, CblasNonUnit, $LT, $L);
  } else {
    gsl_matrix_tricpy('L'.ord, 1, $L, $V);
    gsl_matrix_transpose_tricpy('L'.ord, 1, $LT, $L);
  }
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $L, $LT, 0e0, $A);

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

  my gsl_vector $work = gsl_vector_calloc(3 * $dim);
  my num64 $rcond;
  gsl_linalg_cholesky_rcond($V, $rcond, $work);
  my $*TOLERANCE = 1e-12;
  @res.push: $rcond ≅ 0.03703703703703703;
  gsl_vector_free($work);

  gsl_matrix_free($V);
  gsl_matrix_free($A);
  gsl_matrix_free($L);
  gsl_matrix_free($LT);
  return @res;
}

sub Cholesky-invert(gsl_matrix $m)
{
  my $dim = $m.size1;
  my gsl_matrix $v   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $c   = gsl_matrix_calloc($dim, $dim);
  my gsl_matrix $sol = gsl_matrix_calloc($dim, $dim);
  my Bool @res;
  gsl_matrix_memcpy($v, $m);

  @res.push: gsl_linalg_cholesky_decomp1($v) == GSL_SUCCESS;
  @res.push: gsl_linalg_cholesky_invert($v) == GSL_SUCCESS;

  gsl_blas_dsymm(CblasLeft, CblasUpper, 1e0, $m, $v, 0e0, $c);

  gsl_matrix_set_identity($sol);
  @res.push: test-num64-matrix($c, $sol);

  gsl_matrix_free($v);
  gsl_matrix_free($c);
  gsl_matrix_free($sol);
  return @res;
}

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