### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: lib/Math/Libgsl/LinearAlgebra.rakumod
### ----------------------------------------------------
### -- Chunk 3 of 5
;
}

# Triangular Systems
# v2.5 & v2.6
sub tri-upper-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_tri_upper_rcond($A.matrix, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in tri-upper-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

sub tri-lower-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_tri_lower_rcond($A.matrix, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in tri-lower-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

sub tri-upper-invert(Math::Libgsl::Matrix $T where *.matrix.size1 == $T.matrix.size2 --> Int) is export
{
  gsl_linalg_tri_upper_invert($T.matrix);
}

sub tri-lower-invert(Math::Libgsl::Matrix $T where *.matrix.size1 == $T.matrix.size2 --> Int) is export
{
  gsl_linalg_tri_lower_invert($T.matrix);
}

sub tri-upper-unit-invert(Math::Libgsl::Matrix $T where *.matrix.size1 == $T.matrix.size2 --> Int) is export
{
  gsl_linalg_tri_upper_unit_invert($T.matrix);
}

sub tri-lower-unit-invert(Math::Libgsl::Matrix $T where *.matrix.size1 == $T.matrix.size2 --> Int) is export
{
  gsl_linalg_tri_lower_unit_invert($T.matrix);
}

# v2.6
sub tri-invert(Int $Uplo, Int $Diag, Math::Libgsl::Matrix $T where *.matrix.size1 == $T.matrix.size2 --> Int) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in tri-invert: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_tri_invert($Uplo, $Diag, $T.matrix);
}

sub tri-LTL(Math::Libgsl::Matrix $L where *.matrix.size1 == $L.matrix.size2 --> Int) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in tri-LTL: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_tri_LTL($L.matrix);
}

sub tri-UL(Math::Libgsl::Matrix $LU where *.matrix.size1 == $LU.matrix.size2 --> Int) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in tri-UL: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_tri_UL($LU.matrix);
}

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

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

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

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

sub cholesky-band-invert(Math::Libgsl::Matrix $LLT where *.matrix.size1 == $LLT.matrix.size2
                         --> Math::Libgsl::Matrix) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in cholesky-band-invert: version < v2.6" if $gsl-version < v2.6;
  my Math::Libgsl::Matrix $Ainv .= new: $LLT.matrix.size1, $LLT.matrix.size2;
  my $ret = gsl_linalg_cholesky_band_invert($LLT.matrix, $Ainv.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-band-invert" if $ret ≠ GSL_SUCCESS;
  return $Ainv;
}

sub cholesky-band-unpack(Math::Libgsl::Matrix $LLT where *.matrix.size1 == $LLT.matrix.size2
                         --> Math::Libgsl::Matrix) is export
{
  fail X::Libgsl.new: errno => GSL_FAILURE, error => "Error in cholesky-band-unpack: version < v2.6" if $gsl-version < v2.6;
  my Math::Libgsl::Matrix $L .= new: $LLT.matrix.size1, $LLT.matrix.size2;
  my $ret = gsl_linalg_cholesky_band_unpack($LLT.matrix, $L.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in cholesky-band-unpack" if $ret ≠ GSL_SUCCESS;
  return $L;
}

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

# Banded LDLT Decomposition v2.6
sub ldlt-band-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-band-decomp: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_ldlt_band_decomp($A.matrix);
}

sub ldlt-band-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-band-solve: version < v2.6" if $gsl-version < v2.6;
  my Math::Libgsl::Vector $x .= new: $b.vector.size;
  my $ret = gsl_linalg_ldlt_band_solve($LDLT.matrix, $b.vector, $x.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in ldlt-band-solve" if $ret ≠ GSL_SUCCESS;
  return $x;
}

sub ldlt-band-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-band-svx: version < v2.6" if $gsl-version < v2.6;
  gsl_linalg_ldlt_band_svx($LDLT.matrix, $x.vector);
}

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

sub ldlt-band-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-band-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_band_rcond($LDLT.matrix, $rcond, $work.vector);
  fail X::Libgsl.new: errno => $ret, error => "Error in ldlt-band-rcond" if $ret ≠ GSL_SUCCESS;
  return $rcond;
}

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

=begin pod

=head1 NAME

Math::Libgsl::LinearAlgebra - An interface to libgsl, the Gnu Scientific Library - Linear Algebra.

=head1 SYNOPSIS

=begin code :lang<raku>

use Math::Libgsl::LinearAlgebra;

=end code

=head1 DESCRIPTION

Math::Libgsl::LinearAlgebra is an interface to the linear algebra functions of libgsl, the GNU Scientific Library.
This package provides both the low-level interface to the C library (Raw) and a more comfortable interface layer for the Raku programmer.

This module provides functions for Num and Complex data types.

=head2 Num

=head3 LU-decomp(Math::Libgsl::Matrix $A where *.matrix.size1 == $A.matrix.size2 --> List)

This function factorizes the matrix A into the LU decomposition PA = LU.
The factorization is done in place, so on output the diagonal and upper triangular (or trapezoidal) part of the input matrix A contain the matrix U. The lower triangular (or trapezoidal) part of the input matrix (excluding the diagonal) contains L. The diagonal elements of L are unity, and are not stored.
The return value is a List: the sign of the permutation and a permutation object, which encodes the permutation matrix P.
In case of error a failure object is returned.

=head3 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)

This function solves the square system Ax = b using the LU decomposition of A into (LU, p) given by the output of LU-decomp.
In case of error a failure object is returned.

=head3 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)

This function solves the square system Ax = b in-place using the precomputed LU decomposition of A into (LU, p). On input $x should contain the right-hand side b, which is replaced by the solution on output.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function applies an iterative improvement to x, the solution of Ax = b, from the precomputed LU decomposition of A into (LU, p).
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 LU-invert(Math::Libgsl::Matrix $LU, Math::Libgsl::Permutation $p where *.size == $LU.matrix.size1 --> Math::Libgsl::Matrix)

This function computes the inverse of a matrix A from its LU decomposition (LU, p), returning the matrix inverse.
In case of error a failure object is returned.

=head3 LU-det(Math::Libgsl::Matrix $LU, Int $signum where * ~~ -1|1 --> Num)

This function computes the determinant of a matrix A from its LU decomposition, $LU, and the sign of the permutation, $signum.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 LU-lndet(Math::Libgsl::Matrix $LU --> Num)

This function computes the determinant the logarithm of the absolute value of the determinant of a matrix A, ln |det(A)|, from its LU decomposition, $LU.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 LU-sgndet(Math::Libgsl::Matrix $LU, Int $signum where * ~~ -1|1 --> Int)

This function computes the sign or phase factor of the determinant of a matrix A, det(A)/|det(A)| from its LU decomposition, $LU.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 QR-decomp(Math::Libgsl::Matrix $A --> Math::Libgsl::Vector)

This function factorizes the M-by-N matrix $A into the QR decomposition A = QR. On output the diagonal and upper triangular part of the input matrix contain the matrix R. The returned vector and the columns of the lower triangular part of the matrix A contain the Householder coefficients and Householder vectors which encode the orthogonal matrix Q.
In case of error a failure object is returned.

=head3 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)

This function solves the square system Ax = b using the QR decomposition of A held in ($QR, $tau) which must have been computed previously with QR-decomp().
In case of error a failure object is returned.

=head3 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)

This function solves the square system Ax = b in-place using the QR decomposition of A held in ($QR, $tau) which must have been computed previously by QR-decomp(). On input $x should contain the right-hand side b, which is replaced by the solution on output.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function finds the least squares solution to the overdetermined system Ax = b where the matrix A has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual, ||Ax − b||. The routine requires as input the QR decomposition of A into ($QR, $tau) given by QR-decomp().
The function returns a List of two Math::Libgsl::Vector objects: the solution x and the residual.
In case of error a failure object is returned.

=head3 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)

These function applies the matrix T(Q) encoded in the decomposition (QR, tau) to the vector $v, storing the result T(Q) v in $v. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix T(Q).
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function applies the matrix Q encoded in the decomposition ($QR, $tau) to the vector $v, storing the result Qv in $v. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix Q.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function applies the matrix T(Q) encoded in the decomposition ($QR, $tau) to the M-by-K matrix $B, storing the result T(Q) B in $B. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix T(Q).
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function solves the triangular system Rx = b and returns the Math::Libgsl::Vector object $x.
In case of error a failure object is returned.

=head3 QR-Rsvx(Math::Libgsl::Matrix $QR where *.matrix.size1 == $QR.matrix.size2, Math::Libgsl::Vector $x where *.vector.size == $QR.matrix.size2 --> Int)

This function solves the triangular system Rx = b for x in-place. On input $x should contain the right-hand side b and is replaced by the solution on output.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 QR-unpack(Math::Libgsl::Matrix $QR, Math::Libgsl::Vector $tau where *.vector.size == min($QR.matrix.size1, $QR.matrix.size2) --> List)

This function unpacks the encoded QR decomposition ($QR, $tau) into the matrices Q and R.
The function returns a List of two Math::Libgsl::Matrix objects: $Q which is M-by-M and $R which is M-by-N.
In case of error a failure object is returned.

=head3 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)

This function solves the system Rx = T(Q) b for x. It can be used when the QR decomposition of a matrix is available in unpacked form as ($Q, $R).
The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 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)

This function performs a rank-1 update wT(v) of the QR decomposition ($Q, $R). The update is given by Q'R' = Q(R+wT(v)) where the output matrices $Q and $R are also orthogonal and right triangular. Note that $w is destroyed by the update.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function solves the triangular system Rx = b for the N-by-N matrix $R.
The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 R-svx(Math::Libgsl::Matrix $R where *.matrix.size1 == $R.matrix.size2, Math::Libgsl::Vector $x where *.vector.size == $R.matrix.size2 --> Int)

This function solves the triangular system Rx = b in-place. On input $x should contain the right-hand side b, which is replaced by the solution on output.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 QRPT-decomp(Math::Libgsl::Matrix $A --> List)

This function factorizes the M-by-N matrix $A into the QRT(P) decomposition A = QRT(P).
On output the diagonal and upper triangular part of the input matrix contain the matrix R.
The function's output is a List of three objects: the Math::Libgsl::Vector $tau, the Math::Libgsl::Permutation $p and the sign of the permutation Int $signum.
In case of error a failure object is returned.

=head3 QRPT-decomp2(Math::Libgsl::Matrix $A --> List)

This function factorizes the matrix $A into the decomposition A = QRT(P) without modifying $A itself.
The function returns a List: the Math::Libgsl::Matrix $Q, the Math::Libgsl::Matrix $R, the Math::Libgsl::Permutation $p, and the sign of the permutation Int $signum.
In case of error a failure object is returned.

=head3 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)

This function solves the square system Ax = b using the QRT(P) decomposition of A held in ($QR, $tau, $p) which must have been computed previously by QRPT-decomp.
The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 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)

This function solves the square system Ax = b in-place using the QRT(P) decomposition of A held in ($QR, $tau, $p). On input $x should contain the right-hand side b, which is replaced by the solution on output.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function finds the least squares solution to the overdetermined system Ax = b where the matrix A has more rows than columns and is assumed to have full rank. The least squares solution minimizes the Euclidean norm of the residual, ||b − Ax||. The routine requires as input the QR decomposition of A into ($QR, $tau, $p) given by QRPT-decomp.
The function returns a List of two Math::Libgsl::Vector objects: the solution x and the residual.
In case of error a failure object is returned.

=head3 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)

This function finds the least squares solution to the overdetermined system Ax = b where the matrix A has more rows than columns and has rank given by the input rank. If the user does not know the rank of A, it may be estimated by calling QRPT-rank.
The routine requires as input the QR decomposition of A into ($QR, $tau, $p) given by QRPT-decomp.
The function returns a List of two Math::Libgsl::Vector objects: the solution x and the residual.
In case of error a failure object is returned.

=head3 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)

This function solves the square system RT(P) x = T(Q) b for x. It can be used when the QR decomposition of a matrix is available in unpacked form as ($Q, $R).
The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 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)

This function performs a rank-1 update wT(v) of the QRT(P) decomposition ($Q, $R, $p). The update is given by Q' R' = Q(R + wT(v) P) where the output matrices Q' and R' are also orthogonal and right triangular. Note that $w is destroyed by the update.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 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)

This function solves the triangular system RT(P) x = b for the N-by-N matrix R contained in $QR.
The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 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)

This function solves the triangular system RT(P) x = b in-place for the N-by-N matrix R contained in $QR. On input $x should contain the right-hand side b, which is replaced by the solution on output.
This function returns GSL_SUCCESS if successful, or one of the error codes listed in Math::Libgsl::Constants::gsl-error.

=head3 QRPT-rank(Math::Libgsl::Matrix $QR, Num() $tolerance --> Int)

This function returns the rank of the triangular matrix R contained in $QR.

=head3 QRPT-rcond(Math::Libgsl::Matrix $QR where *.matrix.size1 ≥ $QR.matrix.size2 --> Num)

This function returns the reciprocal condition number (using the 1-norm) of the R factor, stored in the upper triangle of $QR.
In case of error a failure object is returned.

=head3 LQ-decomp(Math::Libgsl::Matrix $A --> Math::Libgsl::Vector)

This function factorizes the M-by-N matrix $A into the LQ decomposition A = LQ. On output the diagonal and lower trapezoidal part of the input matrix contain the matrix L.
This function returns the Math::Libgsl::Vector $tau.
The vector $tau and the elements above the diagonal of the matrix $A contain the Householder coefficients and Householder vectors which encode the orthogonal matrix Q.
In case of error a failure object is returned.

=head3 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)

This function finds the solution to the system Ax = b. The routine requires as input the LQ decomposition of A into ($LQ, $tau) given by LQ-decomp.
The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 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)

This function finds the solution to the system Ax = b. The routine requires as input the LQ decomposition of A into ($LQ, $tau) given by LQ-decomp. On input $x should contain the right-hand side b, which is replaced by the solution on output.
In case of error a failure object is returned.

=head3 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)

This function finds the minimum norm least squares solution to the underdetermined system Ax = b, where the M-by-N matrix A has M ≤ N. The routine requires as input the LQ decomposition of A into ($LQ, $tau) given by LQ-decomp.
The function returns a List of two Math::Libgsl::Vector objects: the solution x and the residual.
In case of error a failure object is returned.

=head3 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)

The function returns the Math::Libgsl::Vector $x.
In case of error a failure object is returned.

=head3 LQ-Lsvx-T(Math::Libgsl::Matrix $LQ where *.matrix.size1 == $LQ.matrix.size2 --> Math::Libgsl::Vector)

On input $x should contain the right-hand side b, which is replaced by