### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 3 of 4
;
  gsl_vector_set($sol, 0, -8e0);
  gsl_vector_set($sol, 1, 18e0);

  @res.push: gsl_linalg_cholesky_decomp1($u) == GSL_SUCCESS;
  @res.push: gsl_linalg_cholesky_solve($u, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

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

sub Cholesky-solve2(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 $D   = 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_cholesky_decomp2($u, $D) == GSL_SUCCESS;
  @res.push: gsl_linalg_cholesky_solve2($u, $D, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

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

sub Cholesky-decomp-unit(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 $l   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $lt  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $dm  = gsl_matrix_calloc($M, $M);
  my gsl_vector $dv  = gsl_vector_calloc($M);
  my Bool @res;
  gsl_matrix_memcpy($v, $m);

  @res.push: gsl_linalg_cholesky_decomp_unit($v, $dv) == GSL_SUCCESS;

  for ^$M -> $i {
    for ^$M -> $j {
      my num64 $vij = gsl_matrix_get($v, $i, $j);
      gsl_matrix_set($l,  $i, $j, $i ≥ $j ?? $vij !! 0e0);
      gsl_matrix_set($lt, $i, $j, $i ≤ $j ?? $vij !! 0e0);
    }
  }
  gsl_matrix_set_zero($dm);
  gsl_matrix_set($dm, $_, $_, gsl_vector_get($dv, $_)) for ^$M;
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $dm, $lt, 0e0, $v);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $l,   $v, 0e0, $a);

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

  gsl_vector_free($dv);
  gsl_matrix_free($v);
  gsl_matrix_free($a);
  gsl_matrix_free($l);
  gsl_matrix_free($lt);
  gsl_matrix_free($dm);
  return @res;
}

sub Cholesky-complex-solve
{
  my CArray[num64] $data .= new: 66e0, 0e0, 0e0, 64e0, 126e0, 63e0, 124e0, -62e0, 61e0, -61e0, 60e0, 60e0, 0e0, -59e0,
                   0e0, -64e0, 65e0, 0e0, 62e0, -124e0, -61e0, -122e0, -60e0, -60e0, 59e0, -59e0, -58e0, 0e0,
                   126e0, -63e0, 62e0, 124e0, 308e0, 0e0, 180e0, -240e0, 59e0, -177e0, 174e0, 58e0, -57e0, -114e0,
                   124e0, 62e0, -61e0, 122e0, 180e0, 240e0, 299e0, 0e0, 174e0, -58e0, 57e0, 171e0, 56e0, -112e0,
                   61e0, 61e0, -60e0, 60e0, 59e0, 177e0, 174e0, 58e0, 119e0, 0e0, 0e0, 112e0, 55e0, -55e0,
                   60e0, -60e0, 59e0, 59e0, 174e0, -58e0, 57e0, -171e0, 0e0, -112e0, 116e0, 0e0, -54e0, -54e0,
                   0e0, 59e0, -58e0, 0e0, -57e0, 114e0, 56e0, 112e0, 55e0, 55e0, -54e0, 54e0, 60e0, 0e0;
  my CArray[num64] $dsol .= new: -0.524944196428570e0, 0.209123883928571e0,
                   1.052873883928572e0, 0.712444196428571e0,
                   0.117568824404762e0, 0.443191964285714e0,
                   0.412862723214286e0, -0.356696428571429e0,
                   0.815931919642858e0, -0.265820312500000e0,
                   0.777929687500000e0, 0.119484747023810e0,
                   1.058733258928571e0, -0.132087053571429e0;
  my gsl_matrix_complex_view $mv = alloc_gsl_matrix_complex_view;
  my gsl_vector_complex_view $vv = alloc_gsl_vector_complex_view;
  my gsl_matrix_complex $m = mgsl_matrix_complex_view_array($mv, $data, 7, 7);
  my gsl_vector_complex $sol = mgsl_vector_complex_view_array($vv, $dsol, 7);
  my $dim = $m.size1;
  my Bool @res;
  my gsl_vector_complex $rhs = gsl_vector_complex_alloc($dim);
  my gsl_matrix_complex $u   = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_vector_complex $x   = gsl_vector_complex_calloc($dim);
  gsl_matrix_complex_memcpy($u, $m);
  my gsl_complex $z = alloc_gsl_complex;
  for ^$dim -> $i {
    mgsl_complex_rect($i + 1e0, 0e0, $z);
    mgsl_vector_complex_set($rhs, $i, $z);
  }
  free_gsl_complex($z);

  @res.push: gsl_linalg_complex_cholesky_decomp($u) == GSL_SUCCESS;
  @res.push: gsl_linalg_complex_cholesky_solve($u, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-complex64-vector($x, $sol);

  free_gsl_matrix_complex_view($mv);
  free_gsl_vector_complex_view($vv);
  gsl_vector_complex_free($x);
  gsl_matrix_complex_free($u);
  gsl_vector_complex_free($rhs);
  return @res;
}

sub Cholesky-complex-decomp
{
  my CArray[num64] $data .= new: 59.75e0, 0e0, 49.25e0, 172.25e0, 66.75e0, -162.75e0,
                                 49.25e0, -172.25e0, 555.5e0, 0e0, -429e0, -333.5e0,
                                 66.75e0, 162.75e0, -429e0, 333.5e0, 536.5e0, 0e0;
  my gsl_matrix_complex_view $mv = alloc_gsl_matrix_complex_view;
  my gsl_matrix_complex $m = mgsl_matrix_complex_view_array($mv, $data, 3, 3);
  my $dim = $m.size1;
  my Bool @res;
  my gsl_matrix_complex $v  = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_matrix_complex $a  = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_matrix_complex $l  = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_matrix_complex $lh = gsl_matrix_complex_alloc($dim, $dim);
  gsl_matrix_complex_memcpy($v, $m);
  gsl_matrix_complex_set_zero($l);
  gsl_matrix_complex_set_zero($lh);

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

  my gsl_complex $vij = alloc_gsl_complex;
  my gsl_complex $cvij = alloc_gsl_complex;
  for ^$dim -> $i {
    for 0..$i -> $j {
      mgsl_matrix_complex_get($v, $i, $j, $vij);
      mgsl_matrix_complex_set($l, $i, $j, $vij);
      mgsl_complex_conjugate($vij, $cvij);
      mgsl_matrix_complex_set($lh, $j, $i, $cvij);
    }
  }
  free_gsl_complex($cvij);
  free_gsl_complex($vij);

  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_zgemm(CblasNoTrans, CblasNoTrans, $one, $l, $lh, $zero, $a);
  free_gsl_complex($zero);
  free_gsl_complex($one);

  my $*TOLERANCE = 1e-12;
  my gsl_complex $aij = alloc_gsl_complex;
  my gsl_complex $mij = alloc_gsl_complex;
  for ^$dim X ^$dim -> ($i, $j) {
    mgsl_matrix_complex_get($a, $i, $j, $aij);
    mgsl_matrix_complex_get($m, $i, $j, $mij);
    @res.push: $aij.dat[0] ≅ $mij.dat[0];
    @res.push: $aij.dat[1] ≅ $mij.dat[1];
  }
  free_gsl_complex($aij);
  free_gsl_complex($mij);

  free_gsl_matrix_complex_view($mv);
  gsl_matrix_complex_free($v);
  gsl_matrix_complex_free($a);
  gsl_matrix_complex_free($l);
  gsl_matrix_complex_free($lh);
  return @res;
}

sub Cholesky-complex-invert
{
  my CArray[num64] $data .= new: 92.303e0, 0.000e0, 10.858e0, 1.798e0, 10.858e0, -1.798e0, 89.027e0, 0.000e0;
  my gsl_matrix_complex_view $mv = alloc_gsl_matrix_complex_view;
  my gsl_matrix_complex $m = mgsl_matrix_complex_view_array($mv, $data, 2, 2);
  my $dim = $m.size1;
  my Bool @res;
  my gsl_matrix_complex $v   = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_matrix_complex $c   = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_matrix_complex $sol = gsl_matrix_complex_alloc($dim, $dim);
  my gsl_complex $af = alloc_gsl_complex;
  my gsl_complex $bt = alloc_gsl_complex;
  gsl_matrix_complex_set_identity($sol);
  gsl_matrix_complex_memcpy($v, $m);

  @res.push: gsl_linalg_complex_cholesky_decomp($v) == GSL_SUCCESS;
  @res.push: gsl_linalg_complex_cholesky_invert($v) == GSL_SUCCESS;

  mgsl_complex_rect(1e0, 0e0, $af);
  mgsl_complex_rect(0e0, 0e0, $bt);
  mgsl_blas_zhemm(CblasLeft, CblasUpper, $af, $m, $v, $bt, $c);

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

  gsl_matrix_complex_free($v);
  gsl_matrix_complex_free($c);
  gsl_matrix_complex_free($sol);
  free_gsl_complex($af);
  free_gsl_complex($bt);
  return @res;
}

sub mCholesky-decomp
{
  my $M = 2;
  my gsl_matrix $m      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $LDLT   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $V      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $A      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $L      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $LT     = gsl_matrix_calloc($M, $M);
  my gsl_vector $S      = gsl_vector_calloc($M);
  my gsl_vector $E      = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my Bool @res;
  create-symm-matrix($m);
  my gsl_vector_view $vv = alloc_gsl_vector_view;
  my gsl_vector $D = mgsl_matrix_diagonal($vv, $LDLT);
  gsl_matrix_memcpy($LDLT, $m);

  @res.push: gsl_linalg_mcholesky_decomp($LDLT, $p, $E) == GSL_SUCCESS;

  @res[1] = True;
  my $*TOLERANCE = 1e-12;
  for ^$M -> $i {
    for ($i + 1)..^$M -> $j {
      @res[1] &&= gsl_matrix_get($m, $i, $j) ≅ gsl_matrix_get($LDLT, $i, $j);
    }
  }

  gsl_matrix_set_identity($L);
  gsl_matrix_set_identity($LT);
  if $gsl-version > v2.5 {
    gsl_matrix_tricpy(CblasLower, CblasUnit, $L, $LDLT);
    gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, $LT, $L);
  } else {
    gsl_matrix_tricpy('L'.ord, 0, $L, $LDLT);
    gsl_matrix_transpose_tricpy('L'.ord, 0, $LT, $L);
  }

  my gsl_vector_view $vv1 = alloc_gsl_vector_view;
  my gsl_vector_view $vv2 = alloc_gsl_vector_view;
  for ^$M -> $i {
    my $v = mgsl_matrix_column($vv1, $L, $i);
    my $w = mgsl_matrix_row($vv2, $LT, $i);
    my num64 $di = sqrt(gsl_vector_get($D, $i));
    gsl_vector_scale($v, $di);
    gsl_vector_scale($w, $di);
  }
  free_gsl_vector_view($vv1);
  free_gsl_vector_view($vv2);

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $L, $LT, 0e0, $A);
  gsl_matrix_memcpy($V, $m);
  $D = mgsl_matrix_diagonal($vv, $V);
  gsl_vector_add($D, $E);

  my gsl_vector_view $vv3 = alloc_gsl_vector_view;
  for ^$M -> $i {
    my gsl_vector $v = mgsl_matrix_row($vv3, $V, $i);
    gsl_permute_vector($p, $v);
  }
  for ^$M -> $i {
    my gsl_vector $v = mgsl_matrix_column($vv3, $V, $i);
    gsl_permute_vector($p, $v);
  }
  free_gsl_vector_view($vv3);

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

  gsl_matrix_free($m);
  gsl_matrix_free($LDLT);
  gsl_matrix_free($V);
  gsl_matrix_free($A);
  gsl_matrix_free($L);
  gsl_matrix_free($LT);
  gsl_vector_free($S);
  gsl_vector_free($E);
  gsl_permutation_free($p);
  free_gsl_vector_view($vv);
  return @res;
}

sub mCholesky-solve
{
  my $M = 2;
  my gsl_matrix $m      = create-hilbert-matrix($M);
  my gsl_vector $rhs    = gsl_vector_calloc($M);
  my gsl_vector $sol    = gsl_vector_calloc($M);
  my gsl_vector $x      = gsl_vector_calloc($M);
  my gsl_vector $S      = gsl_vector_calloc($M);
  my gsl_matrix $u      = gsl_matrix_calloc($M, $M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my gsl_rng $rng       = mgsl_rng_setup(DEFAULT);
  create-random-vector($sol, $rng);
  my Bool @res;

  gsl_blas_dsymv(CblasLower, 1e0, $m, $sol, 0e0, $rhs);
  gsl_matrix_memcpy($u, $m);

  @res.push: gsl_linalg_mcholesky_decomp($u, $p, gsl_vector) == GSL_SUCCESS;
  @res.push: gsl_linalg_mcholesky_solve($u, $p, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_matrix_free($u);
  gsl_matrix_free($m);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_vector_free($x);
  gsl_vector_free($S);
  gsl_rng_free($rng);
  gsl_permutation_free($p);
  return @res;
}

sub mCholesky-invert
{
  my $M = 2;
  my gsl_matrix $m      = create-hilbert-matrix(2);
  my gsl_matrix $v      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $c      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $minv   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $sol    = gsl_matrix_calloc($M, $M);
  my gsl_vector $E      = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my Bool @res;
  gsl_matrix_memcpy($v, $m);
  gsl_matrix_set_identity($sol);

  @res.push: gsl_linalg_mcholesky_decomp($v, $p, $E) == GSL_SUCCESS;
  @res.push: gsl_linalg_mcholesky_invert($v, $p, $minv) == GSL_SUCCESS;
  gsl_blas_dsymm(CblasLeft, CblasUpper, 1e0, $m, $minv, 0e0, $c);

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

  gsl_matrix_free($minv);
  gsl_matrix_free($c);
  gsl_matrix_free($v);
  gsl_matrix_free($m);
  gsl_matrix_free($sol);
  gsl_vector_free($E);
  gsl_permutation_free($p);
  return @res;
}

sub pCholesky-decomp(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $LDLT   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $V      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $A      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $L      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $LT     = gsl_matrix_calloc($M, $M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my Bool @res;
  my gsl_vector_view $vv = alloc_gsl_vector_view;
  my gsl_vector $D = mgsl_matrix_diagonal($vv, $LDLT);
  gsl_matrix_memcpy($LDLT, $m);

  @res.push: gsl_linalg_pcholesky_decomp($LDLT, $p) == GSL_SUCCESS;

  my $*TOLERANCE = 1e-12;
  @res[1] = True;
  for ^$M -> $i {
    for ($i + 1)..^$M -> $j {
      @res[1] &&= gsl_matrix_get($m, $i, $j) ≅ gsl_matrix_get($LDLT, $i, $j);
    }
  }

  gsl_matrix_set_identity($L);
  gsl_matrix_set_identity($LT);
  if $gsl-version > v2.5 {
    gsl_matrix_tricpy(CblasLower, CblasUnit, $L, $LDLT);
    gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, $LT, $L);
  } else {
    gsl_matrix_tricpy('L'.ord, 0, $L, $LDLT);
    gsl_matrix_transpose_tricpy('L'.ord, 0, $LT, $L);
  }

  my gsl_vector_view $vv1 = alloc_gsl_vector_view;
  my gsl_vector_view $vv2 = alloc_gsl_vector_view;
  for ^$M -> $i {
    my $v = mgsl_matrix_column($vv1, $L, $i);
    my $w = mgsl_matrix_row($vv2, $LT, $i);
    my num64 $di = gsl_vector_get($D, $i);
    gsl_vector_scale($v, sqrt $di);
    gsl_vector_scale($w, sqrt $di);
  }
  free_gsl_vector_view($vv1);
  free_gsl_vector_view($vv2);

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $L, $LT, 0e0, $A);
  gsl_matrix_memcpy($V, $m);

  my gsl_vector_view $vv3 = alloc_gsl_vector_view;
  for ^$M -> $i {
    my gsl_vector $v = mgsl_matrix_row($vv3, $V, $i);
    gsl_permute_vector($p, $v);
  }
  for ^$M -> $i {
    my gsl_vector $v = mgsl_matrix_column($vv3, $V, $i);
    gsl_permute_vector($p, $v);
  }
  free_gsl_vector_view($vv3);

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

  my gsl_vector $work = gsl_vector_calloc(3 * $M);
  my num64 $rcond;
  gsl_linalg_pcholesky_rcond($LDLT, $p, $rcond, $work);
  @res.push: $rcond ≅ 0.03703703703703703;

  gsl_matrix_free($LDLT);
  gsl_matrix_free($V);
  gsl_matrix_free($A);
  gsl_matrix_free($L);
  gsl_matrix_free($LT);
  gsl_vector_free($work);
  gsl_permutation_free($p);
  free_gsl_vector_view($vv);
  return @res;
}

sub pCholesky-solve(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $u      = gsl_matrix_calloc($M, $M);
  my gsl_vector $x      = gsl_vector_calloc($M);
  my gsl_vector $sol    = gsl_vector_calloc($M);
  my gsl_vector $rhs    = gsl_vector_calloc($M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my Bool @res;
  create-random-vector($sol, $r);
  gsl_blas_dsymv(CblasLower, 1e0, $m, $sol, 0e0, $rhs);
  gsl_matrix_memcpy($u, $m);

  @res.push: gsl_linalg_pcholesky_decomp($u, $p) == GSL_SUCCESS;
  @res.push: gsl_linalg_pcholesky_solve($u, $p, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_rng_free($r);
  gsl_matrix_free($u);
  gsl_vector_free($x);
  gsl_vector_free($sol);
  gsl_permutation_free($p);
  return @res;
}

sub pCholesky-invert(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $v      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $c      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $minv   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $sol    = gsl_matrix_calloc($M, $M);
  my gsl_permutation $p = gsl_permutation_alloc($M);
  my Bool @res;
  gsl_matrix_memcpy($v, $m);
  gsl_matrix_set_identity($sol);

  @res.push: gsl_linalg_pcholesky_decomp($v, $p) == GSL_SUCCESS;
  @res.push: gsl_linalg_pcholesky_invert($v, $p, $minv) == GSL_SUCCESS;

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

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

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

sub bidiag-decomp(gsl_matrix $m)
{
  my $M = $m.size1;
  my gsl_matrix $A      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $a      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $b      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $u      = gsl_matrix_calloc($M, $M);
  my gsl_matrix $v      = gsl_matrix_calloc($M, $M);
  my gsl_vector $tau1   = gsl_vector_calloc($M);
  my gsl_vector $tau2   = gsl_vector_calloc($M - 1);
  my gsl_vector $d      = gsl_vector_calloc($M);
  my gsl_vector $sd     = gsl_vector_calloc($M - 1);
  my Bool @res;
  gsl_matrix_memcpy($A, $m);

  @res.push: gsl_linalg_bidiag_decomp($A, $tau1, $tau2) == GSL_SUCCESS;
  @res.push: gsl_linalg_bidiag_unpack($A, $tau1, $u, $tau2, $v, $d, $sd) == GSL_SUCCESS;

  gsl_matrix_set_zero($b);
  gsl_matrix_set($b, $_, $_, gsl_vector_get($d, $_)) for ^$M;
  gsl_matrix_set($b, $_, $_ + 1, gsl_vector_get($sd, $_)) for ^($M - 1);
  for ^$M X ^$M -> ($i, $j) {
    my num64 $sum = 0e0;
    for ^$M X ^$M -> ($k, $r) {
      $sum += gsl_matrix_get($u, $i, $k) * gsl_matrix_get($b, $k, $r) * gsl_matrix_get($v, $j, $r);
    }
    gsl_matrix_set($a, $i, $j, $sum);
  }

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

  gsl_matrix_free($A);
  gsl_matrix_free($a);
  gsl_matrix_free($b);
  gsl_matrix_free($u);
  gsl_matrix_free($v);
  gsl_vector_free($tau1);
  gsl_vector_free($tau2);
  gsl_vector_free($d);
  gsl_vector_free($sd);
  return @res;
}

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

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

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

sub tridiag-symm-solve
{
  my $M = 2;
  my gsl_vector $offdiag = gsl_vector_calloc($M - 1);
  my gsl_vector $diag    = gsl_vector_calloc($M);
  my gsl_vector $rhs     = gsl_vector_calloc($M);
  my gsl_vector $x       = gsl_vector_calloc($M);
  my gsl_vector $sol     = gsl_vector_calloc($M);
  my Bool @res;
  gsl_vector_set($sol, 0, 0e0);
  gsl_vector_set($sol, 1, 2e0);
  for ^$M -> $i {
    gsl_vector_set($diag, $i, 1e0);
    gsl_vector_set($rhs, $i, $i + 1e0);
  }
  for ^($M - 1) -> $i {
    gsl_vector_set($offdiag, $i, .5e0);
  }

  @res.push: gsl_linalg_solve_symm_tridiag($diag, $offdiag, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($offdiag);
  gsl_vector_free($diag);
  gsl_vector_free($rhs);
  gsl_vector_free($x);
  gsl_vector_free($sol);
  return @res;
}

sub tridiag-symm-cyc-solve
{
  my $M = 3;
  my gsl_vector $offdiag = gsl_vector_calloc($M);
  my gsl_vector $diag    = gsl_vector_calloc($M);
  my gsl_vector $rhs     = gsl_vector_calloc($M);
  my gsl_vector $x       = gsl_vector_calloc($M);
  my gsl_vector $sol     = gsl_vector_calloc($M);
  my Bool @res;
  gsl_vector_set($diag, 0, 1e0);
  gsl_vector_set($diag, 1, 1e0);
  gsl_vector_set($diag, 2, 1e0);
  gsl_vector_set($offdiag, 0, 3e0);
  gsl_vector_set($offdiag, 1, 3e0);
  gsl_vector_set($offdiag, 2, 3e0);
  gsl_vector_set($rhs, 0,  7e0);
  gsl_vector_set($rhs, 1, -7e0);
  gsl_vector_set($rhs, 2,  7e0);
  gsl_vector_set($sol, 0, -2e0);
  gsl_vector_set($sol, 1,  5e0);
  gsl_vector_set($sol, 2, -2e0);

  @res.push: gsl_linalg_solve_symm_cyc_tridiag($diag, $offdiag, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($offdiag);
  gsl_vector_free($diag);
  gsl_vector_free($rhs);
  gsl_vector_free($x);
  gsl_vector_free($sol);
  return @res;
}

sub tridiag-solve
{
  my $M = 3;
  my gsl_vector $abovediag = gsl_vector_calloc($M - 1);
  my gsl_vector $belowdiag = gsl_vector_calloc($M - 1);
  my gsl_vector $diag    = gsl_vector_calloc($M);
  my gsl_vector $rhs     = gsl_vector_calloc($M);
  my gsl_vector $x       = gsl_vector_calloc($M);
  my gsl_vector $sol     = gsl_vector_calloc($M);
  my Bool @res;
  gsl_vector_set($sol, 0, .75e0);
  gsl_vector_set($sol, 1, .75e0);
  gsl_vector_set($sol, 2, 2.625e0);
  for ^$M -> $i {
    gsl_vector_set($diag, $i, 1e0);
    gsl_vector_set($rhs,  $i, $i + 1e0);
  }
  for ^($M - 1) -> $i {
    gsl_vector_set($abovediag, $i, (1/3).Num);
    gsl_vector_set($belowdiag, $i, .5e0);
  }

  @res.push: gsl_linalg_solve_tridiag($diag, $abovediag, $belowdiag, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($abovediag);
  gsl_vector_free($belowdiag);
  gsl_vector_free($diag);
  gsl_vector_free