### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: lib/Math/Libgsl/LinearAlgebra.rakumod
### ----------------------------------------------------
### -- Chunk 2 of 5
);
  fail X::Libgsl.new: errno => $ret, error => "Error in SV-decomp" if $ret ≠ GSL_SUCCESS;
  return $V, $S;
}

sub SV-decomp-mod(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2 --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Matrix $V    .= new: $N, $N;
  my Math::Libgsl::Matrix $X    .= new: $N, $N;
  my Math::Libgsl::Vector $S    .= new: $N;
  my Math::Libgsl::Vector $work .= new: $N;
  my $ret = gsl_linalg_SV_decomp_mod($A.matrix, $X.matrix, $V.matrix, $S.vector, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in SV-decomp-mod" if $ret ≠ GSL_SUCCESS;
  return $V, $S;
}

sub SV-decomp-jacobi(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2 --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Matrix $V    .= new: $N, $N;
  my Math::Libgsl::Vector $S    .= new: $N;
  my $ret = gsl_linalg_SV_decomp_jacobi($A.matrix, $V.matrix, $S.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in SV-decomp-jacobi" if $ret ≠ GSL_SUCCESS;
  return $V, $S;
}

sub SV-solve(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2,
             Math::Libgsl::Matrix $V where { $V.matrix.size1 == $V.matrix.size2 && $V.matrix.size1 == $A.matrix.size2 },
             Math::Libgsl::Vector $S where *.vector.size == $A.matrix.size2,
             Math::Libgsl::Vector $b where *.vector.size == $A.matrix.size1
             --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $A.matrix.size2;
  my $ret = gsl_linalg_SV_solve($A.matrix, $V.matrix, $S.vector, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in SV-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub SV-leverage(Math::Libgsl::Matrix $A --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $h .= new: $A.matrix.size1;
  my $ret = gsl_linalg_SV_leverage($A.matrix, $h.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in SV-leverage" if $ret ≠ GSL_SUCCESS;
  return $h;
}

# Cholesky Decomposition
sub cholesky-decomp1(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Int) is export
{
  gsl_linalg_cholesky_decomp1($A.matrix);
}

sub cholesky-solve(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                   Math::Libgsl::Vector $b where *.vector.size == $A.matrix.size1 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $A.matrix.size2;
  my $ret = gsl_linalg_cholesky_solve($A.matrix, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub cholesky-svx(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                 Math::Libgsl::Vector $x where *.vector.size == $A.matrix.size1 --> Int) is export
{
  gsl_linalg_cholesky_svx($A.matrix, $x.vector);
}

sub cholesky-solve-mat(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                       Math::Libgsl::Matrix $B where *.matrix.size1 == $A.matrix.size1
                       --> Math::Libgsl::Matrix) is export
{
  my Math::Libgsl::Matrix $X .= new: $B.matrix.size1, $B.matrix.size2;
  my $ret = gsl_linalg_cholesky_solve_mat($A.matrix, $B.matrix, $X.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-solve-mat" if $ret ≠ GSL_SUCCESS;
  return $X;
}

sub cholesky-svx-mat(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                     Math::Libgsl::Matrix $X where *.matrix.size1 == $A.matrix.size1 --> Int) is export
{
  gsl_linalg_cholesky_svx_mat($A.matrix, $X.matrix);
}

sub cholesky-invert(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Int) is export
{
  gsl_linalg_cholesky_invert($A.matrix);
}

sub cholesky-decomp2(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $s .= new: $A.matrix.size1;
  my $ret = gsl_linalg_cholesky_decomp2($A.matrix, $s.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-decomp2" if $ret ≠ GSL_SUCCESS;
  return $s;
}

sub cholesky-solve2(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                    Math::Libgsl::Vector $s where *.vector.size == $A.matrix.size1,
                    Math::Libgsl::Vector $b where *.vector.size == $A.matrix.size1 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $A.matrix.size2;
  my $ret = gsl_linalg_cholesky_solve2($A.matrix, $s.vector, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-solve2" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub cholesky-svx2(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                  Math::Libgsl::Vector $s where *.size == $A.matrix.size2,
                  Math::Libgsl::Vector $x where *.vector.size == $A.matrix.size2 --> Int) is export
{
  gsl_linalg_cholesky_svx2($A.matrix, $s.vector, $x.vector);
}

sub cholesky-decomp-unit(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $d .= new: $A.matrix.size1;
  my $ret = gsl_linalg_cholesky_decomp_unit($A.matrix, $d.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-decomp-unit" if $ret ≠ GSL_SUCCESS;
  return $d;
}

sub cholesky-scale(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $s .= new: $A.matrix.size2;
  my $ret = gsl_linalg_cholesky_scale($A.matrix, $s.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-scale" if $ret ≠ GSL_SUCCESS;
  return $s;
}

sub cholesky-scale-apply(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                         Math::Libgsl::Vector $s where *.vector.size == $A.matrix.size2 --> Int) is export
{
  gsl_linalg_cholesky_scale_apply($A.matrix, $s.vector);
}

sub cholesky-rcond(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Num) is export
{
  my Math::Libgsl::Vector $work .= new: 3 * $A.matrix.size2;
  my num64 $rcond;
  my $ret = gsl_linalg_cholesky_rcond($A.matrix, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

# Pivoted Cholesky Decomposition
sub pcholesky-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2
                     --> Math::Libgsl::Permutation) is export
{
  my Math::Libgsl::Permutation $p .= new: $A.matrix.size1;
  my $ret = gsl_linalg_pcholesky_decomp($A.matrix, $p.p);
  fail X::Libgsl.new: errno => $ret, error => "Error in pcholesky-decomp" if $ret ≠ GSL_SUCCESS;
  return $p;
}

sub pcholesky-solve(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                    Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1,
                    Math::Libgsl::Vector $b where *.vector.size == $LDLT.matrix.size1
                    --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $LDLT.matrix.size2;
  my $ret = gsl_linalg_pcholesky_solve($LDLT.matrix, $p.p, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in pcholesky-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub pcholesky-svx(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                  Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1,
                  Math::Libgsl::Vector $x where *.vector.size == $LDLT.matrix.size1 --> Int) is export
{
  gsl_linalg_pcholesky_svx($LDLT.matrix, $p.p, $x.vector);
}

sub pcholesky-decomp2(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> List) is export
{
  my Math::Libgsl::Permutation $p .= new: $A.matrix.size1;
  my Math::Libgsl::Vector $s .= new: $A.matrix.size2;
  my $ret = gsl_linalg_pcholesky_decomp2($A.matrix, $p.p, $s.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in pcholesky-decomp2" if $ret ≠ GSL_SUCCESS;
  return $p, $s;
}

sub pcholesky-solve2(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                     Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1,
                     Math::Libgsl::Vector $s where *.vector.size == $LDLT.matrix.size1,
                     Math::Libgsl::Vector $b where *.vector.size == $LDLT.matrix.size1
                     --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $LDLT.matrix.size2;
  my $ret = gsl_linalg_pcholesky_solve2($LDLT.matrix, $p.p, $s.vector, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in pcholesky-solve2" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub pcholesky-svx2(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                   Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1,
                   Math::Libgsl::Vector $s where *.vector.size == $LDLT.matrix.size1,
                   Math::Libgsl::Vector $x where *.vector.size == $LDLT.matrix.size1 --> Int) is export
{
  gsl_linalg_pcholesky_svx2($LDLT.matrix, $p.p, $s.vector, $x.vector);
}

sub pcholesky-invert(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                     Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1
                     --> Math::Libgsl::Matrix) is export
{
  my Math::Libgsl::Matrix $Ainv .= new: $LDLT.matrix.size1, $LDLT.matrix.size2;
  my $ret = gsl_linalg_pcholesky_invert($LDLT.matrix, $p.p, $Ainv.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in pcholesky-invert" if $ret ≠ GSL_SUCCESS;
  return $Ainv;
}

sub pcholesky-rcond(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                    Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1
                    --> Num) is export
{
  my Math::Libgsl::Vector $work .= new: 3 * $LDLT.matrix.size2;
  my num64 $rcond;
  my $ret = gsl_linalg_pcholesky_rcond($LDLT.matrix, $p.p, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in pcholesky-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

# Modified Cholesky Decomposition
sub mcholesky-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                     Bool :$perturbation = True
                     --> List) is export
{
  my Math::Libgsl::Permutation $p .= new: $A.matrix.size1;
  my Math::Libgsl::Vector $E;
  $E .= new($A.matrix.size1) if $perturbation;
  my $ret = gsl_linalg_mcholesky_decomp($A.matrix, $p.p, $perturbation ?? $E.vector !! Math::Libgsl::Vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in mcholesky-decomp" if $ret ≠ GSL_SUCCESS;
  return $p, $E;
}

sub mcholesky-solve(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                    Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1,
                    Math::Libgsl::Vector $b where *.vector.size == $LDLT.matrix.size1
                    --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $LDLT.matrix.size2;
  my $ret = gsl_linalg_mcholesky_solve($LDLT.matrix, $p.p, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in mcholesky-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub mcholesky-svx(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                  Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1,
                  Math::Libgsl::Vector $x where *.vector.size == $LDLT.matrix.size1 --> Int) is export
{
  gsl_linalg_mcholesky_svx($LDLT.matrix, $p.p, $x.vector);
}

sub mcholesky-rcond(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                    Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1
                    --> Num) is export
{
  my Math::Libgsl::Vector $work .= new: 3 * $LDLT.matrix.size2;
  my num64 $rcond;
  my $ret = gsl_linalg_mcholesky_rcond($LDLT.matrix, $p.p, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in mcholesky-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

sub mcholesky-invert(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
                     Math::Libgsl::Permutation $p where *.p.size == $LDLT.matrix.size1
                     --> Math::Libgsl::Matrix) is export
{
  my Math::Libgsl::Matrix $Ainv .= new: $LDLT.matrix.size1, $LDLT.matrix.size2;
  my $ret = gsl_linalg_mcholesky_invert($LDLT.matrix, $p.p, $Ainv.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in mcholesky-invert" if $ret ≠ GSL_SUCCESS;
  return $Ainv;
}

# LDLT Decomposition
sub ldlt-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Int) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in ldlt-decomp: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_ldlt_decomp($A.matrix);
}

sub ldlt-solve(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
               Math::Libgsl::Vector $b where *.vector.size == $LDLT.matrix.size1
               --> Math::Libgsl::Vector) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in ldlt-solve: version < v2.6" if $gsl-version < v2.6;
  my Math::Libgsl::Vector $x .= new: $LDLT.matrix.size2;
  my $ret = gsl_linalg_ldlt_solve($LDLT.matrix, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in ldlt-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub ldlt-svx(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2,
             Math::Libgsl::Vector $x where *.vector.size == $LDLT.matrix.size1 --> Int) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in ldlt-svx: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_ldlt_svx($LDLT.matrix, $x.vector);
}

sub ldlt-rcond(Math::Libgsl::Matrix $LDLT where *.matrix.size1 == $LDLT.matrix.size2 --> Num) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in ldlt-rcond: version < v2.6" if $gsl-version < v2.6;
  my Math::Libgsl::Vector $work .= new: 3 * $LDLT.matrix.size2;
  my num64 $rcond;
  my $ret = gsl_linalg_ldlt_rcond($LDLT.matrix, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in ldlt-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

# Tridiagonal Decomposition of Real Symmetric Matrices
sub symmtd-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $tau .= new: $A.matrix.size1 - 1;
  my $ret = gsl_linalg_symmtd_decomp($A.matrix, $tau.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in symmtd-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau;
}

sub symmtd-unpack(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                  Math::Libgsl::Vector $tau where *.vector.size == $A.matrix.size1 - 1 --> List) is export
{
  my $M = $A.matrix.size1;
  my Math::Libgsl::Matrix $Q     .= new: $M, $M;
  my Math::Libgsl::Vector $diag  .= new: $A.matrix.size1;
  my Math::Libgsl::Vector $sdiag .= new: $A.matrix.size1 - 1;
  my $ret = gsl_linalg_symmtd_unpack($A.matrix, $tau.vector, $Q.matrix, $diag.vector, $sdiag.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in symmtd-unpack" if $ret ≠ GSL_SUCCESS;
  return $Q, $diag, $sdiag;
}

sub symmtd-unpack-T(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> List) is export
{
  my $M = $A.matrix.size1;
  my Math::Libgsl::Vector $diag  .= new: $A.matrix.size1;
  my Math::Libgsl::Vector $sdiag .= new: $A.matrix.size1 - 1;
  my $ret = gsl_linalg_symmtd_unpack_T($A.matrix, $diag.vector, $sdiag.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in symmtd-unpack-T" if $ret ≠ GSL_SUCCESS;
  return $diag, $sdiag;
}

# Hessenberg Decomposition of Real Matrices
sub hessenberg-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $tau .= new: $A.matrix.size1;
  my $ret = gsl_linalg_hessenberg_decomp($A.matrix, $tau.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in hessenberg-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau;
}

sub hessenberg-unpack(Math::Libgsl::Matrix $H where *.matrix.size1 == $H.matrix.size2,
                      Math::Libgsl::Vector $tau where *.vector.size == $H.matrix.size1
                      --> Math::Libgsl::Matrix) is export
{
  my $M = $H.matrix.size1;
  my Math::Libgsl::Matrix $U .= new: $M, $M;
  my $ret = gsl_linalg_hessenberg_unpack($H.matrix, $tau.vector, $U.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in hessenberg-unpack" if $ret ≠ GSL_SUCCESS;
  return $U;
}

sub hessenberg-unpack-accum(Math::Libgsl::Matrix $H where *.matrix.size1 == $H.matrix.size2,
                            Math::Libgsl::Vector $tau where *.vector.size == $H.matrix.size1
                            --> Math::Libgsl::Matrix) is export
{
  my $M = $H.matrix.size1;
  my Math::Libgsl::Matrix $V .= new: $M, $M;
  my $ret = gsl_linalg_hessenberg_unpack_accum($H.matrix, $tau.vector, $V.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in hessenberg-unpack-accum" if $ret ≠ GSL_SUCCESS;
  return $V;
}

sub hessenberg-set-zero(Math::Libgsl::Matrix $H where *.matrix.size1 == $H.matrix.size2 --> Int) is export
{
  gsl_linalg_hessenberg_set_zero($H.matrix);
}

# Hessenberg-Triangular Decomposition of Real Matrices
sub hesstri-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
                   Math::Libgsl::Matrix $B where { $B.matrix.size1 == $A.matrix.size1 && $B.matrix.size2 == $A.matrix.size2},
                   Bool :$similarity = True
                   --> List) is export
{
  my $M = $A.matrix.size1;
  my Math::Libgsl::Vector $work .= new: $M;
  my Math::Libgsl::Matrix $U;
  my Math::Libgsl::Matrix $V;
  if $similarity {
    $U .= new: $M, $M;
    $V .= new: $M, $M;
  }
  my $ret = gsl_linalg_hesstri_decomp($A.matrix,
                                      $B.matrix,
                                      $similarity ?? $U.matrix !! Math::Libgsl::Matrix,
                                      $similarity ?? $V.matrix !! Math::Libgsl::Matrix,
                                      $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in hesstri-decomp" if $ret ≠ GSL_SUCCESS;
  return $U, $V;
}

# Bidiagonalization
sub bidiag-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2 --> List) is export
{
  my Math::Libgsl::Vector $tau_U .= new: $A.matrix.size2;
  my Math::Libgsl::Vector $tau_V .= new: $A.matrix.size2 - 1;
  my $ret = gsl_linalg_bidiag_decomp($A.matrix, $tau_U.vector, $tau_V.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in bidiag-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau_U, $tau_V;
}

sub bidiag-unpack(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2,
                  Math::Libgsl::Vector $tau_U where *.vector.size == $A.matrix.size2,
                  Math::Libgsl::Vector $tau_V where *.vector.size == $A.matrix.size2 - 1
                  --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Matrix $U     .= new: $M, $N;
  my Math::Libgsl::Matrix $V     .= new: $N, $N;
  my Math::Libgsl::Vector $diag  .= new: $N;
  my Math::Libgsl::Vector $sdiag .= new: $N - 1;
  my $ret = gsl_linalg_bidiag_unpack($A.matrix, $tau_U.vector, $U.matrix, $tau_V.vector, $V.matrix, $diag.vector, $sdiag.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in bidiag-unpack" if $ret ≠ GSL_SUCCESS;
  return $U, $V, $diag, $sdiag;
}

sub bidiag-unpack2(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2,
                   Math::Libgsl::Vector $tau_U where *.vector.size == $A.matrix.size2,
                   Math::Libgsl::Vector $tau_V where *.vector.size == $A.matrix.size2 - 1
                   --> Math::Libgsl::Matrix) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Matrix $V .= new: $N, $N;
  my $ret = gsl_linalg_bidiag_unpack2($A.matrix, $tau_U.vector, $tau_V.vector, $V.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in bidiag-unpack2" if $ret ≠ GSL_SUCCESS;
  return $V;
}

sub bidiag-unpack-B(Math::Libgsl::Matrix $A where *.matrix.size1 ≥ $A.matrix.size2 --> List) is export
{
  my $N = $A.matrix.size2;
  my Math::Libgsl::Vector $diag  .= new: $N;
  my Math::Libgsl::Vector $sdiag .= new: $N - 1;
  my $ret = gsl_linalg_bidiag_unpack_B($A.matrix, $diag.vector, $sdiag.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in bidiag-unpack-B" if $ret ≠ GSL_SUCCESS;
  return $diag, $sdiag;
}

# Givens Rotations
sub givens(Num() $a, Num() $b --> List) is export
{
  my num64 $A = $a;
  my num64 $B = $b;
  my num64 ($c, $s);
  gsl_linalg_givens($A, $B, $c, $s);
  return $c, $s;
}

sub givens-gv(Math::Libgsl::Vector $v, Int, $i, Int $j, Num() $c, Num() $s) is export
{
  my num64 $C = $c;
  my num64 $S = $s;
  gsl_linalg_givens_gv($v.vector, $i, $j, $c, $s);
}

# Householder Transformations
sub householder-transform(Math::Libgsl::Vector $w --> Num) is export
{
  gsl_linalg_householder_transform($w.vector);
}

sub householder-hm(Num() $tau, Math::Libgsl::Vector $v, Math::Libgsl::Matrix $A --> Int) is export
{
  my num64 $τ = $tau;
  gsl_linalg_householder_hm($τ, $v.vector, $A.matrix);
}

sub householder-mh(Num() $tau, Math::Libgsl::Vector $v, Math::Libgsl::Matrix $A --> Int) is export
{
  my num64 $τ = $tau;
  gsl_linalg_householder_mh($τ, $v.vector, $A.matrix);
}

sub householder-hv(Num() $tau, Math::Libgsl::Vector $v, Math::Libgsl::Vector $w --> Int) is export
{
  my num64 $τ = $tau;
  gsl_linalg_householder_hv($τ, $v.vector, $w.vector);
}

# Householder solver for linear systems
sub HH-solve(Math::Libgsl::Matrix $A where *.matrix.size1 ≤ $A.matrix.size2,
             Math::Libgsl::Vector $b where *.vector.size == $A.matrix.size1
             --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $A.matrix.size2;
  my $ret = gsl_linalg_HH_solve($A.matrix, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in HH-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub HH-svx(Math::Libgsl::Matrix $A where *.matrix.size1 ≤ $A.matrix.size2,
           Math::Libgsl::Vector $x where *.vector.size == $A.matrix.size2 --> Int) is export
{
  gsl_linalg_HH_svx($A.matrix, $x.vector);
}

# Tridiagonal Systems
sub tridiag-solve(Math::Libgsl::Vector $diag,
                  Math::Libgsl::Vector $abovediag where *.vector.size == $diag.vector.size - 1,
                  Math::Libgsl::Vector $belowdiag where *.vector.size == $diag.vector.size - 1,
                  Math::Libgsl::Vector $rhs where *.vector.size == $diag.vector.size
                  --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $diag.vector.size;
  my $ret = gsl_linalg_solve_tridiag($diag.vector, $abovediag.vector, $belowdiag.vector, $rhs.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in tridiag-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub tridiag-symm-solve(Math::Libgsl::Vector $diag,
                       Math::Libgsl::Vector $offdiag where *.vector.size == $diag.vector.size - 1,
                       Math::Libgsl::Vector $rhs where *.vector.size == $diag.vector.size
                       --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $diag.vector.size;
  my $ret = gsl_linalg_solve_symm_tridiag($diag.vector, $offdiag.vector, $rhs.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in tridiag-symm-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub tridiag-cyc-solve(Math::Libgsl::Vector $diag where *.vector.size ≥ 3,
                      Math::Libgsl::Vector $abovediag where *.vector.size == $diag.vector.size,
                      Math::Libgsl::Vector $belowdiag where *.vector.size == $diag.vector.size,
                      Math::Libgsl::Vector $rhs where *.vector.size == $diag.vector.size
                      --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $diag.vector.size;
  my $ret = gsl_linalg_solve_cyc_tridiag($diag.vector, $abovediag.vector, $belowdiag.vector, $rhs.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in tridiag-cyc-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub tridiag-symm-cyc-solve(Math::Libgsl::Vector $diag where *.vector.size ≥ 3,
                           Math::Libgsl::Vector $offdiag where *.vector.size == $diag.vector.size,
                           Math::Libgsl::Vector $rhs where *.vector.size == $diag.vector.size
                           --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $diag.vector.size;
  my $ret = gsl_linalg_solve_symm_cyc_tridiag($diag.vector, $offdiag.vector, $rhs.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in tridiag-symm-cyc-solve" if $ret ≠ GSL_SUCCESS;
  return $x