### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: lib/Math/Libgsl/LinearAlgebra.rakumod
### ----------------------------------------------------
### -- Chunk 1 of 5
use v6.c;

unit class Math::Libgsl::LinearAlgebra:ver<0.0.4>:auth<zef:FRITH>;

use Math::Libgsl::Raw::LinearAlgebra :ALL;
use NativeCall;
use Math::Libgsl::Vector;
use Math::Libgsl::Matrix;
use Math::Libgsl::Permutation;
use Math::Libgsl::Exception;
use Math::Libgsl::Constants;

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

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

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

sub LU-refine(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2,
              Math::Libgsl::Matrix $LU where { $LU.matrix.size1 == $LU.matrix.size2 && $A.matrix.size1 == $LU.matrix.size2 },
              Math::Libgsl::Permutation $p where *.size == $LU.matrix.size1,
              Math::Libgsl::Vector $b where *.vector.size == $LU.matrix.size1,
              Math::Libgsl::Vector $x where *.vector.size == $LU.matrix.size1 --> Int) is export
{
  my Math::Libgsl::Vector $work .= new($LU.matrix.size1);
  gsl_linalg_LU_refine($A.matrix, $LU.matrix, $p.p, $b.vector, $x.vector, $work.vector);
}

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

sub LU-det(Math::Libgsl::Matrix $LU, Int $signum where * ~~ -1|1 --> Num) is export
{
  gsl_linalg_LU_det($LU.matrix, $signum);
}

sub LU-lndet(Math::Libgsl::Matrix $LU --> Num) is export
{
  gsl_linalg_LU_lndet($LU.matrix);
}

sub LU-sgndet(Math::Libgsl::Matrix $LU, Int $signum where * ~~ -1|1 --> Int) is export
{
  gsl_linalg_LU_sgndet($LU.matrix, $signum);
}

# QR Decomposition
sub QR-decomp(Math::Libgsl::Matrix $A --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $tau .= new: min($A.matrix.size1, $A.matrix.size2);
  my $ret = gsl_linalg_QR_decomp($A.matrix, $tau.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in QR-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau;
}

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

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

sub QR-lssolve(Math::Libgsl::Matrix $QR where *.matrix.size1 ≥ $QR.matrix.size2,
               Math::Libgsl::Vector $tau where *.vector.size == $QR.matrix.size2,
               Math::Libgsl::Vector $b where *.vector.size == $QR.matrix.size1 --> List) is export
{
  my Math::Libgsl::Vector $x        .= new: $QR.matrix.size2;
  my Math::Libgsl::Vector $residual .= new: $QR.matrix.size1;
  my $ret = gsl_linalg_QR_lssolve($QR.matrix, $tau.vector, $b.vector, $x.vector, $residual.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in QR-lssolve" if $ret ≠ GSL_SUCCESS;
  return $x, $residual;
}

sub QR-QTvec(Math::Libgsl::Matrix $QR,
             Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2),
             Math::Libgsl::Vector $v where *.vector.size == $QR.matrix.size1 --> Int) is export
{
  gsl_linalg_QR_QTvec($QR.matrix, $tau.vector, $v.vector);
}

sub QR-Qvec(Math::Libgsl::Matrix $QR,
            Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2),
            Math::Libgsl::Vector $v where *.vector.size == $QR.matrix.size1 --> Int) is export
{
  gsl_linalg_QR_Qvec($QR.matrix, $tau.vector, $v.vector);
}

sub QR-QTmat(Math::Libgsl::Matrix $QR,
             Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2),
             Math::Libgsl::Matrix $B where *.matrix.size1 == $QR.matrix.size1 --> Int) is export
{
  gsl_linalg_QR_QTmat($QR.matrix, $tau.vector, $B.matrix);
}

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

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

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

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

sub QR-update(Math::Libgsl::Matrix $Q,
              Math::Libgsl::Matrix $R where { $Q.matrix.size1 == $R.matrix.size1 && $Q.matrix.size2 == $R.matrix.size1 },
              Math::Libgsl::Vector $w where *.vector.size == $R.matrix.size1,
              Math::Libgsl::Vector $v where *.vector.size == $R.matrix.size2 --> Int) is export
{
  gsl_linalg_QR_update($Q.matrix, $R.matrix, $w.vector, $v.vector);
}

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

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

# QR Decomposition with Column Pivoting
sub QRPT-decomp(Math::Libgsl::Matrix $A --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Vector $tau    .= new: min($M, $N);
  my Math::Libgsl::Vector $norm   .= new: $N;
  my Math::Libgsl::Permutation $p .= new: $N;
  my int32 $signum;
  my $ret = gsl_linalg_QRPT_decomp($A.matrix, $tau.vector, $p.p, $signum, $norm.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in QRPT-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau, $p, $signum;
}

sub QRPT-decomp2(Math::Libgsl::Matrix $A --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Matrix $Q      .= new: $M, $M;
  my Math::Libgsl::Matrix $R      .= new: $M, $N;
  my Math::Libgsl::Vector $tau    .= new: min($M, $N);
  my Math::Libgsl::Vector $norm   .= new: $N;
  my Math::Libgsl::Permutation $p .= new: $N;
  my int32 $signum;
  my $ret = gsl_linalg_QRPT_decomp2($A.matrix, $Q.matrix, $R.matrix, $tau.vector, $p.p, $signum, $norm.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in QRPT-decomp2" if $ret ≠ GSL_SUCCESS;
  return $Q, $R, $p, $signum;
}

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

sub QRPT-svx(Math::Libgsl::Matrix $QR where *.matrix.size1 == $QR.matrix.size2,
             Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2),
             Math::Libgsl::Permutation $p where *.p.size == $QR.matrix.size1,
             Math::Libgsl::Vector $x where *.vector.size == $QR.matrix.size2 --> Int) is export
{
  gsl_linalg_QRPT_svx($QR.matrix, $tau.vector, $p.p, $x.vector);
}

sub QRPT-lssolve(Math::Libgsl::Matrix $QR where *.matrix.size1 ≥ $QR.matrix.size2,
               Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2),
               Math::Libgsl::Permutation $p where *.p.size == $QR.matrix.size2,
               Math::Libgsl::Vector $b where *.vector.size == $QR.matrix.size1 --> List) is export
{
  my Math::Libgsl::Vector $x        .= new: $QR.matrix.size2;
  my Math::Libgsl::Vector $residual .= new: $QR.matrix.size1;
  my $ret = gsl_linalg_QRPT_lssolve($QR.matrix, $tau.vector, $p.p, $b.vector, $x.vector, $residual.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in QRPT-lssolve" if $ret ≠ GSL_SUCCESS;
  return $x, $residual;
}

sub QRPT-lssolve2(Math::Libgsl::Matrix $QR where *.matrix.size1 ≥ $QR.matrix.size2,
               Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2),
               Math::Libgsl::Permutation $p where *.p.size == $QR.matrix.size2,
               Math::Libgsl::Vector $b where *.vector.size == $QR.matrix.size1,
               Int $rank where 0 < * ≤ $QR.matrix.size2 --> List) is export
{
  my Math::Libgsl::Vector $x        .= new: $QR.matrix.size2;
  my Math::Libgsl::Vector $residual .= new: $QR.matrix.size1;
  my $ret = gsl_linalg_QRPT_lssolve2($QR.matrix, $tau.vector, $p.p, $b.vector, $rank, $x.vector, $residual.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in QRPT-lssolve2" if $ret ≠ GSL_SUCCESS;
  return $x, $residual;
}

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

sub QRPT-update(Math::Libgsl::Matrix $Q,
                Math::Libgsl::Matrix $R where { $Q.matrix.size1 == $R.matrix.size1 && $Q.matrix.size2 == $R.matrix.size1 },
                Math::Libgsl::Permutation $p where *.p.size == $R.matrix.size1,
                Math::Libgsl::Vector $w where *.vector.size == $R.matrix.size1,
                Math::Libgsl::Vector $v where *.vector.size == $R.matrix.size2 --> Int) is export
{
  gsl_linalg_QRPT_update($Q.matrix, $R.matrix, $p.p, $w.vector, $v.vector);
}

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

sub QRPT-Rsvx(Math::Libgsl::Matrix $QR where *.matrix.size1 == $QR.matrix.size2,
              Math::Libgsl::Permutation $p where *.p.size == $QR.matrix.size2,
              Math::Libgsl::Vector $x where *.vector.size == $QR.matrix.size2 --> Int) is export
{
  gsl_linalg_QRPT_Rsvx($QR.matrix, $p.p, $x.vector);
}

sub QRPT-rank(Math::Libgsl::Matrix $QR, Num() $tolerance --> Int) is export
{
  my num64 $tol = $tolerance;
  gsl_linalg_QRPT_rank($QR.matrix, $tol);
}

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

# LQ Decomposition
sub LQ-decomp(Math::Libgsl::Matrix $A --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $tau .= new: min($A.matrix.size1, $A.matrix.size2);
  my $ret = gsl_linalg_LQ_decomp($A.matrix, $tau.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in LQ-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau;
}

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

sub LQ-svx-T(Math::Libgsl::Matrix $LQ where *.matrix.size1 == $LQ.matrix.size2,
               Math::Libgsl::Vector $tau where *.vector.size == min($LQ.matrix.size1, $LQ.matrix.size2),
               --> Math::Libgsl::Vector) is export
{
  my Math::Libgsl::Vector $x .= new: $LQ.matrix.size1;
  my $ret = gsl_linalg_LQ_svx_T($LQ.matrix, $tau.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in LQ-svx-T" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub LQ-lssolve-T(Math::Libgsl::Matrix $LQ where *.matrix.size1 ≥ $LQ.matrix.size2,
               Math::Libgsl::Vector $tau where *.vector.size == min($LQ.matrix.size1, $LQ.matrix.size2),
               Math::Libgsl::Vector $b where *.vector.size == $LQ.matrix.size1 --> List) is export
{
  my Math::Libgsl::Vector $x        .= new: $LQ.matrix.size2;
  my Math::Libgsl::Vector $residual .= new: $LQ.matrix.size1;
  my $ret = gsl_linalg_LQ_lssolve_T($LQ.matrix, $tau.vector, $b.vector, $x.vector, $residual.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in LQ-lssolve-T" if $ret ≠ GSL_SUCCESS;
  return $x, $residual;
}

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

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

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

sub LQ-vecQ(Math::Libgsl::Matrix $LQ,
            Math::Libgsl::Vector $tau where *.vector.size == min($LQ.matrix.size1, $LQ.matrix.size2),
            Math::Libgsl::Vector $v where *.vector.size == $LQ.matrix.size1 --> Int) is export
{
  gsl_linalg_LQ_vecQ($LQ.matrix, $tau.vector, $v.vector);
}

sub LQ-vecQT(Math::Libgsl::Matrix $LQ,
             Math::Libgsl::Vector $tau where *.vector.size == min($LQ.matrix.size1, $LQ.matrix.size2),
             Math::Libgsl::Vector $v where *.vector.size == $LQ.matrix.size1 --> Int) is export
{
  gsl_linalg_LQ_vecQT($LQ.matrix, $tau.vector, $v.vector);
}

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

sub LQ-update(Math::Libgsl::Matrix $Q,
              Math::Libgsl::Matrix $L where { $L.matrix.size2 == $Q.matrix.size1 && $L.matrix.size2 == $Q.matrix.size2 },
              Math::Libgsl::Vector $v where *.vector.size == $L.matrix.size1,
              Math::Libgsl::Vector $w where *.vector.size == $L.matrix.size2
              --> Int) is export
{
  gsl_linalg_LQ_update($Q.matrix, $L.matrix, $v.vector, $w.vector);
}

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

# Complete Orthogonal Decomposition
sub COD-decomp(Math::Libgsl::Matrix $A --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Vector $tau-Q .= new: min($M, $N);
  my Math::Libgsl::Vector $tau-Z .= new: min($M, $N);
  my Math::Libgsl::Vector $work   .= new: $N;
  my Math::Libgsl::Permutation $p .= new: $N;
  my size_t $rank;
  my $ret = gsl_linalg_COD_decomp($A.matrix, $tau-Q.vector, $tau-Z.vector, $p.p, $rank, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in COD-decomp" if $ret ≠ GSL_SUCCESS;
  return $tau-Q, $tau-Z, $p, $rank;
}

sub COD-decomp-e(Math::Libgsl::Matrix $A, Num() $tolerance --> List) is export
{
  my $M = $A.matrix.size1;
  my $N = $A.matrix.size2;
  my Math::Libgsl::Vector $tau-Q .= new: min($M, $N);
  my Math::Libgsl::Vector $tau-Z .= new: min($M, $N);
  my Math::Libgsl::Vector $work   .= new: $N;
  my Math::Libgsl::Permutation $p .= new: $N;
  my Int $rank;
  my num64 $tol = $tolerance;
  my $ret = gsl_linalg_COD_decomp_e($A.matrix, $tau-Q.vector, $tau-Z.vector, $p.p, $tol, $rank, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in COD-decomp-e" if $ret ≠ GSL_SUCCESS;
  return $tau-Q, $tau-Z, $p, $rank;
}

sub COD-lssolve(Math::Libgsl::Matrix $QRZT where *.matrix.size1 ≥ $QRZT.matrix.size2,
                Math::Libgsl::Vector $tau-Q where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
                Math::Libgsl::Vector $tau-Z where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
                Math::Libgsl::Permutation $p where *.p.size == $QRZT.matrix.size2,
                Int $rank where * ≤ min($QRZT.matrix.size1, $QRZT.matrix.size2),
                Math::Libgsl::Vector $b where *.vector.size == $QRZT.matrix.size1 --> List) is export
{
  my Math::Libgsl::Vector $x        .= new: $QRZT.matrix.size2;
  my Math::Libgsl::Vector $residual .= new: $QRZT.matrix.size1;
  my $ret = gsl_linalg_COD_lssolve($QRZT.matrix, $tau-Q.vector, $tau-Z.vector, $p.p, $rank,
                                   $b.vector, $x.vector, $residual.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in COD-lssolve" if $ret ≠ GSL_SUCCESS;
  return $x, $residual;
}

sub COD-lssolve2(Math::Libgsl::Matrix $QRZT where *.matrix.size1 ≥ $QRZT.matrix.size2,
                 Math::Libgsl::Vector $tau-Q where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
                 Math::Libgsl::Vector $tau-Z where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
                 Math::Libgsl::Permutation $p where *.p.size == $QRZT.matrix.size2,
                 Int $rank where * ≤ min($QRZT.matrix.size1, $QRZT.matrix.size2),
                 Math::Libgsl::Vector $b where *.vector.size == $QRZT.matrix.size1,
                 Num() $lambda --> List) is export
{
  my Math::Libgsl::Vector $x        .= new: $QRZT.matrix.size2;
  my Math::Libgsl::Vector $residual .= new: $QRZT.matrix.size1;
  my Math::Libgsl::Matrix $S        .= new: $rank, $rank;
  my Math::Libgsl::Vector $work     .= new: $rank;
  my $ret = gsl_linalg_COD_lssolve2($lambda, $QRZT.matrix, $tau-Q.vector, $tau-Z.vector, $p.p, $rank,
                                    $b.vector, $x.vector, $residual.vector, $S.matrix, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in COD-lssolve2" if $ret ≠ GSL_SUCCESS;
  return $x, $residual;
}

sub COD-unpack(Math::Libgsl::Matrix $QRZT,
               Math::Libgsl::Vector $tau-Q where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
               Math::Libgsl::Vector $tau-Z where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
               Int $rank where * ≤ min($QRZT.matrix.size1, $QRZT.matrix.size2) --> List) is export
{
  my $M = $QRZT.matrix.size1;
  my $N = $QRZT.matrix.size2;
  my Math::Libgsl::Matrix $Q .= new: $M, $M;
  my Math::Libgsl::Matrix $R .= new: $M, $N;
  my Math::Libgsl::Matrix $Z .= new: $N, $N;
  my $ret = gsl_linalg_COD_unpack($QRZT.matrix, $tau-Q.vector, $tau-Z.vector, $rank, $Q.matrix, $R.matrix, $Z.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in COD-unpack" if $ret ≠ GSL_SUCCESS;
  return $Q, $R, $Z;
}

sub COD-matZ(Math::Libgsl::Matrix $QRZT,
             Math::Libgsl::Vector $tau-Z where *.vector.size == min($QRZT.matrix.size1, $QRZT.matrix.size2),
             Math::Libgsl::Matrix $A where *.matrix.size2 == $QRZT.matrix.size2,
             Int $rank --> Int) is export
{
  my Math::Libgsl::Vector $work .= new: $A.matrix.size1;
  gsl_linalg_COD_matZ($QRZT.matrix, $tau-Z.vector, $rank, $A.matrix, $work.vector);
}

# Singular Value Decomposition
sub SV-decomp(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 Math::Libgsl::Vector $work   .= new: $N;
  my $ret = gsl_linalg_SV_decomp($A.matrix, $V.matrix, $S.vector, $work.vector