PDL::L...lgebra(3pm)

📦

version 0.437 · 6 versions

💖 NAME 💖

okay so bestie, welcome to PDL::LinearAlgebra(3pm) — it lives in the User Contributed Perl Documentation section, and it is exactly what it sounds like: Linear Algebra utils for PDL 🎀✨

🦋 SYNOPSIS 🦋

 use PDL::LinearAlgebra;

 $a = random (100,100);
 ($U, $s, $V) = mdsvd($a);

like. three lines and you're already doing a singular value decomposition. iconic. 💅

🌸 DESCRIPTION 🌸

This module is your convenient little interface to PDL::LinearAlgebra::Real and PDL::LinearAlgebra::Complex. 🫶

Now here's the tea bestie: BLAS and LAPACK use a column major ordering scheme, so some of the routines in here have to transpose matrices before they call the fortran routines and then transpose back afterwards (each routine's own documentation tells you when this is happening, check it 👀).

And if you want optimized code? Go straight to PDL::LinearAlgebra::Real and PDL::LinearAlgebra::Complex directly. No cap. 💯

🔥 FUNCTIONS 🔥

✨ setlaerror ✨

Sets the action type used when an error gets encountered, and it returns the previous type. The available values are NO, WARN and BARF (these are predefined constants, bestie, they're just there for you 💖).

So like — say you're computing an inverse and the routine detects singularity. Depending on what you set, it can:

silently return values from the computation (see the manuals for what you actually get 🙈)
WARN you about the singularity
or straight up BARF 😭

And BARF is the default value. She's dramatic. Periodt.

 # h : x -> g(f(x))

 $a = sequence(5,5);
 $err = setlaerror(NO);
 ($b, $info)= f($a);
 setlaerror($err);
 $info ? barf "can't compute h" : return g($b);

see how you save the old error setting, do your thing, then put it back? main character error handling. 👑

🎀 getlaerror 🎀

Gets the action type used when an error is encountered. That's it, that's the whole vibe:

  0 => NO,
  1 => WARN,
  2 => BARF

💜 issym 💜

 PDL = issym(PDL, SCALAR|PDL(tol),SCALAR(hermitian))
 tol : tolerance value, default: 1e-8 for double else 1e-5
 hermitian : Hermitian = 1 | Symmetric = 0, default = 0;

Checks the symmetricity/Hermitianicity of a matrix. 🦋

tol — the tolerance value, default is 1e-8 for double, otherwise 1e-5
hermitian — Hermitian = 1 or Symmetric = 0, default is 0

Supports broadcasting. ✨

🌈 diag 🌈

Returns the i-th diagonal if you hand it a matrix, or a matrix with the i-th diagonal if you hand it an entry. And the i-th diagonal it returns flows data back&forth — so yes bestie, it can be used as an lvalue sub if your perl supports it. 👀

Supports broadcasting.

 PDL = diag(PDL, SCALAR(i), SCALAR(vector)))
 i  : i-th diagonal, default = 0
 vector  : create diagonal matrices by broadcasting over row vectors, default = 0

i — which diagonal, default 0
vector — creates diagonal matrices by broadcasting over row vectors, default 0

 my $a = random(5,5);
 my $diag  = diag($a,2);
 # If your perl support lvaluable subroutines.
 $a->diag(-2) .= pdl(1,2,3);
 # Construct a (5,5,5) PDL (5 matrices) with
 # diagonals from row vectors of $a
 $a->diag(0,1)

the .= assignment into a diagonal? that's the lvalue sub flexing. obsessed. 💅

🫶 positivise 🫶

Returns the entry pdl with the sign changed by row, so that the average of the positive signs is > 0. In other words: it broadcasts among dimension 1, and does row = -row whenever sum(sign(row)) < 0. 🧚

Heads up: this only makes sense for real ndarrays.

Works inplace. 🔥

 my $a = random(10,10);
 $a -= 0.5;
 $a->xchg(0,1)->inplace->positivise;

⭐ mcrossprod ⭐

Computes the cross-product of two matrices: A' x B. If you only give it one matrix, it takes B to be the same as A. Supports broadcasting. Uses crossprod or ccrossprod.

 PDL = mcrossprod(PDL(A), (PDL(B))

 my $a = random(10,10);
 my $crossproduct = mcrossprod($a);

💫 mrank 💫

Computes the rank of a matrix using a singular value decomposition, and returns it as a plain Perl scalar. Straight from LAPACK, bestie. ✨

 SCALAR = mrank(PDL, SCALAR(TOL))
 TOL:  tolerance value, default : mnorm(dims(PDL),'inf') * mnorm(PDL) * EPS

TOL is the tolerance value, and its default is mnorm(dims(PDL),'inf') * mnorm(PDL) * EPS.

 my $a = random(10,10);
 my $b = mrank($a, 1e-5);

🌟 mnorm 🌟

Computes the norm of a real or complex matrix. Supports broadcasting. 😌

 PDL(norm) = mnorm(PDL, SCALAR(ord));
 ord :
  0|'inf' : Infinity norm
  1|'one' : One norm
  2|'two'  : norm 2 (default)
  3|'fro' : frobenius norm

0 or 'inf' — Infinity norm
1 or 'one' — One norm
2 or 'two' — norm 2, and this one is the default
3 or 'fro' — frobenius norm

 my $a = random(10,10);
 my $norm = mnorm($a);

🎯 mdet 🎯

Computes the determinant of a general square matrix using LU factorization. Supports broadcasting. Uses getrf or cgetrf from LAPACK.

 PDL(determinant) = mdet(PDL);

 my $a = random(10,10);
 my $det = mdet($a);

🍓 mposdet 🍓

Computes the determinant of a symmetric or Hermitian positive definite square matrix using Cholesky factorization. Supports broadcasting. Uses potrf or cpotrf from LAPACK.

 (PDL, PDL) = mposdet(PDL, SCALAR)
 SCALAR : UPPER = 0 | LOWER = 1, default = 0

so the scalar is UPPER = 0 or LOWER = 1, and the default is 0.

 my $a = random(10,10);
 my $det = mposdet($a);

🌺 mcond 🌺

Computes the condition number (two-norm) of a general matrix.

The condition number in two-n is defined as:

 norm (a) * norm (inv (a)).

Uses a singular value decomposition. Supports broadcasting. 💅

 PDL = mcond(PDL)

 my $a = random(10,10);
 my $cond = mcond($a);

💀 mrcond 💀

Estimates the reciprocal condition number of a general square matrix using LU factorization, in either the 1-norm or the infinity-norm.

The reciprocal condition number is defined as:

 1/(norm (a) * norm (inv (a)))

Supports broadcasting. And heads up: it works on transposed array(s).

 PDL = mrcond(PDL, SCALAR(ord))
 ord :
  0 : Infinity norm (default)
  1 : One norm

0 — Infinity norm, and this is the default
1 — One norm

 my $a = random(10,10);
 my $rcond = mrcond($a,1);

🦄 morth 🦄

Returns an orthonormal basis of the range space of matrix A. Also, important: no broadcast, bestie. She does not do that.

 PDL = morth(PDL(A), SCALAR(tol))
 tol : tolerance for determining rank, default: 1e-8 for double else 1e-5

tol is the tolerance for determining rank, defaulting to 1e-8 for double, else 1e-5.

 my $a = sequence(10,10);
 my $ortho = morth($a, 1e-8);

🧚 mnull 🧚

Returns an orthonormal basis of the null space of matrix A. And again: no broadcast.

 PDL = mnull(PDL(A), SCALAR(tol))
 tol : tolerance for determining rank, default: 1e-8 for double else 1e-5

tol is the tolerance for determining rank — 1e-8 for double, else 1e-5.

 my $a = sequence(10,10);
 my $null = mnull($a, 1e-8);

👑 minv 👑

Computes the inverse of a general square matrix using LU factorization. Supports inplace and broadcasting. Uses getrf and getri or cgetrf and cgetri from LAPACK, and it returns "inverse, info" in array context. 👀

 PDL(inv)  = minv(PDL)

 my $a = random(10,10);
 my $inv = minv($a);

🌈 mtriinv 🌈

Computes the inverse of a triangular matrix. Supports inplace and broadcasting. Uses trtri or ctrtri from LAPACK. Returns "inverse, info" in array context.

 (PDL, PDL(info))) = mtriinv(PDL, SCALAR(uplo), SCALAR|PDL(diag))
 uplo : UPPER = 0 | LOWER = 1, default = 0
 diag : UNITARY DIAGONAL = 1, default = 0

uplo — UPPER = 0 or LOWER = 1, default 0
diag — UNITARY DIAGONAL = 1, default 0

 # Assume $a is upper triangular
 my $a = random(10,10);
 my $inv = mtriinv($a);

💖 msyminv 💖

Computes the inverse of a symmetric square matrix using the Bunch-Kaufman diagonal pivoting method. Supports inplace and broadcasting. Uses sytrf and sytri or csytrf and csytri from LAPACK, and returns "inverse, info" in array context.

 (PDL, (PDL(info))) = msyminv(PDL, SCALAR|PDL(uplo))
 uplo : UPPER = 0 | LOWER = 1, default = 0

uplo is UPPER = 0 or LOWER = 1, default 0.

 # Assume $a is symmetric
 my $a = random(10,10);
 my $inv = msyminv($a);

🎀 mposinv 🎀

Computes the inverse of a symmetric positive definite square matrix using Cholesky factorization. Supports inplace and broadcasting. Uses potrf and potri or cpotrf and cpotri from LAPACK, and returns "inverse, info" in array context.

 (PDL, (PDL(info))) = mposinv(PDL, SCALAR|PDL(uplo))
 uplo : UPPER = 0 | LOWER = 1, default = 0

uplo is UPPER = 0 or LOWER = 1, default 0.

 # Assume $a is symmetric positive definite
 my $a = random(10,10);
 $a = $a->crossprod($a);
 my $inv = mposinv($a);

✨ mpinv ✨

Computes the pseudo-inverse (Moore-Penrose) of a general matrix. And note: it works on a transposed array.

 PDL(pseudo-inv)  = mpinv(PDL, SCALAR(tol))
 TOL:  tolerance value, default : mnorm(dims(PDL),'inf') * mnorm(PDL) * EPS

TOL is the tolerance value, defaulting to mnorm(dims(PDL),'inf') * mnorm(PDL) * EPS.

 my $a = random(5,10);
 my $inv = mpinv($a);

🧚 mlu 🧚

Computes LU factorization. Uses getrf or cgetrf from LAPACK and returns L, U, pivot and info. Works on a transposed array.

 (PDL(l), PDL(u), PDL(pivot), PDL(info)) = mlu(PDL)

 my $a = random(10,10);
 ($l, $u, $pivot, $info) = mlu($a);

🌸 mchol 🌸

Computes the Cholesky decomposition of a symmetric matrix, also known as the symmetric square root. If the inplace flag is set, it overwrites the leading upper or lower triangular part of A; otherwise it returns a triangular matrix. Returns "cholesky, info" in array context. Supports broadcasting. Uses potrf or cpotrf from LAPACK. 💅

 PDL(Cholesky) = mchol(PDL, SCALAR)
 SCALAR : UPPER = 0 | LOWER = 1, default = 0

the scalar is UPPER = 0 or LOWER = 1, default 0.

 my $a = random(10,10);
 $a = crossprod($a, $a);
 my $u  = mchol($a);

⭐ mhessen ⭐

Reduces a square matrix to Hessenberg form H and orthogonal matrix Q.

It reduces a general matrix A to upper Hessenberg form H by an orthogonal similarity transformation:

 Q' x A x Q = H

 A = Q x H x Q'

Uses gehrd and orghr or cgehrd and cunghr from LAPACK, and returns "H" in scalar context, else "H" and "Q". Works on a transposed array.

 (PDL(h), (PDL(q))) = mhessen(PDL)

 my $a = random(10,10);
 ($h, $q) = mhessen($a);

💫 mschur 💫

Computes the Schur form, and yes, it works inplace.

 A = Z x T x Z'

Supports broadcasting for unordered eigenvalues. Uses gees or cgees from LAPACK and returns schur(T) in scalar context. Works on transposed array(s).

 ( PDL(schur), (PDL(eigenvalues), (PDL(left schur vectors), PDL(right schur vectors), $sdim), $info) ) = mschur(PDL(A), SCALAR(schur vector),SCALAR(left eigenvector), SCALAR(right eigenvector),SCALAR(select_func), SCALAR(backtransform), SCALAR(norm))
 schur vector       : Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
 left eigenvector    : Left eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
 right eigenvector   : Right eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
 select_func       : Select_func is used to select eigenvalues to sort
           to the top left of the Schur form.
           An eigenvalue is selected if PerlInt select_func(w) is true;
           (the inputs are converted to complex ndarrays for you)
           Note that a selected complex eigenvalue may no longer
           satisfy select_func(w) = 1 after ordering, since
           ordering may change the value of complex eigenvalues
           (especially if the eigenvalue is ill-conditioned).
           All eigenvalues/vectors are selected if select_func is undefined.
 backtransform       : Whether or not backtransforms eigenvectors to those of A.
           Only supported if schur vectors are computed, default = 1.
 norm                : Whether or not computed eigenvectors are normalized to have Euclidean norm equal to
           1 and largest component real, default = 1

 Returned values     :
           Schur form T (SCALAR CONTEXT),
           eigenvalues,
           Schur vectors (Z) if requested,
           left eigenvectors if requested
           right eigenvectors if requested
           sdim: Number of eigenvalues selected if select_func is defined.
           info: Info output from gees/cgees.

let's break the args down bestie, one bullet each 🫶

schur vector — are Schur vectors returned? none = 0, all = 1, selected = 2, default 0
left eigenvector — left eigenvectors returned? none = 0, all = 1, selected = 2, default 0
right eigenvector — right eigenvectors returned? none = 0, all = 1, selected = 2, default 0
select_func — used to select eigenvalues to sort to the top left of the Schur form. An eigenvalue is selected if PerlInt select_func(w) is true (the inputs are converted to complex ndarrays for you 💖). Note that a selected complex eigenvalue may no longer satisfy select_func(w) = 1 after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned 😭). All eigenvalues/vectors are selected if select_func is undefined.
backtransform — whether or not to backtransform eigenvectors to those of A. Only supported if Schur vectors are computed. Default 1.
norm — whether or not computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. Default 1.

And the returned values (just to restate the tea):

Schur form T (in SCALAR CONTEXT)
eigenvalues
Schur vectors (Z) if requested
left eigenvectors if requested
right eigenvectors if requested
sdim: the number of eigenvalues selected if select_func is defined
info: the info output from gees/cgees

 my $a = random(10,10);
 my $schur  = mschur($a);
 sub select {
  my $w = shift;
  # select "discrete time" eigenspace
  return $w->Cabs < 1 ? 1 : 0;
 }
 my ($schur,$eigen,$svectors,$evectors)  = mschur($a,1,1,0,\&select);

that select sub is lowkey elegant, she's out here picking the "discrete time" eigenspace and I respect it. 🎀

🦋 mschurx 🦋

Computes the Schur form, and it works inplace. Uses geesx or cgeesx from LAPACK and returns schur(T) in scalar context. Works on a transposed array.

`mschurx` 🌸 (the rest of it, bestie)

( PDL(schur) (,PDL(eigenvalues))  (, PDL(schur vectors), HASH(result)) ) = mschurx(PDL, SCALAR(schur vector), SCALAR(left eigenvector), SCALAR(right eigenvector),SCALAR(select_func), SCALAR(sense), SCALAR(backtransform), SCALAR(norm))

Okay so here are the arguments, one by one, no skips 💖

schur vector — Schur vectors returned: none = 0 | all = 1 | selected = 2, default = 0
left eigenvector — Left eigenvectors returned: none = 0 | all = 1 | selected = 2, default = 0
right eigenvector — Right eigenvectors returned: none = 0 | all = 1 | selected = 2, default = 0
select_func — select_func is used to select eigenvalues to sort to the top left of the Schur form. An eigenvalue is selected if PerlInt select_func(w) is true; (the inputs are converted to complex ndarrays for you) Note that a selected complex eigenvalue may no longer satisfy select_func(w) = 1 after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned). All eigenvalues/vectors are selected if select_func is undefined.
sense — Determines which reciprocal condition numbers will be computed.
- 0: None are computed
- 1: Computed for average of selected eigenvalues only
- 2: Computed for selected right invariant subspace only
- 3: Computed for both If select_func is undefined, sense is not used.
backtransform — Whether or not backtransforms eigenvectors to those of A. Only supported if schur vector are computed, default = 1
norm — Whether or not computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real, default = 1

Returned values ✨

Schur form T (SCALAR CONTEXT),
eigenvalues,
Schur vectors if requested,
HASH{VL}: left eigenvectors if requested
HASH{VR}: right eigenvectors if requested
HASH{info}: info output from gees/cgees.
if select_func is defined:
- HASH{n}: number of eigenvalues selected,
- HASH{rconde}: reciprocal condition numbers for the average of the selected eigenvalues if requested,
- HASH{rcondv}: reciprocal condition numbers for the selected right invariant subspace if requested.

my $a = random(10,10);
my $schur  = mschurx($a);
sub select {
 my $m = shift;
 # select "discrete time" eigenspace
 return $m->Cabs < 1 ? 1 : 0;
}
my ($schur,$eigen, $vectors,%ret)  = mschurx($a,1,0,0,\&select);

`mgschur` 🦋

Computes generalized Schur decomposition of the pair (A,B).

 A = Q x S x Z'
 B = Q x T x Z'

Uses gges or cgges from LAPACK. Works on transposed array. (iykyk 💅)

( PDL(schur S), PDL(schur T), PDL(alpha), PDL(beta), HASH{result}) = mgschur(PDL(A), PDL(B), SCALAR(left schur vector),SCALAR(right schur vector),SCALAR(left eigenvector), SCALAR(right eigenvector), SCALAR(select_func), SCALAR(backtransform), SCALAR(scale))

left schur vector — Left Schur vectors returned: none = 0 | all = 1 | selected = 2, default = 0
right schur vector — Right Schur vectors returned: none = 0 | all = 1 | selected = 2, default = 0
left eigenvector — Left eigenvectors returned: none = 0 | all = 1 | selected = 2, default = 0
right eigenvector — Right eigenvectors returned: none = 0 | all = 1 | selected = 2, default = 0
select_func — Select_func is used to select eigenvalues to sort to the top left of the Schur form. An eigenvalue w = wr(j)+sqrt(-1)*wi(j) is selected if PerlInt select_func(alpha,beta) is true; (the inputs are converted to complex ndarrays for you) Note that a selected complex eigenvalue may no longer satisfy select_func = 1 after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned). All eigenvalues/vectors are selected if select_func is undefined.
backtransform — Whether or not backtransforms eigenvectors to those of (A,B). Only supported if right and/or left schur vector are computed,
scale — Whether or not computed eigenvectors are scaled so the largest component will have abs(real part) + abs(imag. part) = 1, default = 1

Returned values 💜

Schur form S,
Schur form T,
alpha,
beta (eigenvalues = alpha/beta),
HASH{info}: info output from gges/cgges.
HASH{SL}: left Schur vectors if requested
HASH{SR}: right Schur vectors if requested
HASH{VL}: left eigenvectors if requested
HASH{VR}: right eigenvectors if requested
HASH{n} : Number of eigenvalues selected if select_func is defined.

my $a = random(10,10);
my $b = random(10,10);
my ($S,$T) = mgschur($a,$b);
sub select {
 my ($alpha,$beta) = @_;
 return $alpha->Cabs < abs($beta) ? 1 : 0;
}
my ($S, $T, $alpha, $beta, %res)  = mgschur( $a, $b, 1, 1, 1, 1,\&select);

`mgschurx` 💅

Computes generalized Schur decomposition of the pair (A,B).

 A = Q x S x Z'
 B = Q x T x Z'

Uses ggesx or cggesx from LAPACK. Works on transposed array.

( PDL(schur S), PDL(schur T), PDL(alpha), PDL(beta), HASH{result}) = mgschurx(PDL(A), PDL(B), SCALAR(left schur vector),SCALAR(right schur vector),SCALAR(left eigenvector), SCALAR(right eigenvector), SCALAR(select_func), SCALAR(sense), SCALAR(backtransform), SCALAR(scale))

left schur vector — Left Schur vectors returned: none = 0 | all = 1 | selected = 2, default = 0
right schur vector — Right Schur vectors returned: none = 0 | all = 1 | selected = 2, default = 0
left eigenvector — Left eigenvectors returned: none = 0 | all = 1 | selected = 2, default = 0
right eigenvector — Right eigenvectors returned: none = 0 | all = 1 | selected = 2, default = 0
select_func — Select_func is used to select eigenvalues to sort to the top left of the Schur form. An eigenvalue w = wr(j)+sqrt(-1)*wi(j) is selected if PerlInt select_func(alpha,beta) is true; (the inputs are converted to complex ndarrays for you) Note that a selected complex eigenvalue may no longer satisfy select_func = 1 after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned). All eigenvalues/vectors are selected if select_func is undefined.
sense — Determines which reciprocal condition numbers will be computed.
- 0: None are computed
- 1: Computed for average of selected eigenvalues only
- 2: Computed for selected deflating subspaces only
- 3: Computed for both If select_func is undefined, sense is not used.
backtransform — Whether or not backtransforms eigenvectors to those of (A,B). Only supported if right and/or left schur vector are computed, default = 1
scale — Whether or not computed eigenvectors are scaled so the largest component will have abs(real part) + abs(imag. part) = 1, default = 1

Returned values 🌈

Schur form S,
Schur form T,
alpha,
beta (eigenvalues = alpha/beta),
HASH{info}: info output from gges/cgges.
HASH{SL}: left Schur vectors if requested
HASH{SR}: right Schur vectors if requested
HASH{VL}: left eigenvectors if requested
HASH{VR}: right eigenvectors if requested
HASH{rconde}: reciprocal condition numbers for average of selected eigenvalues if requested
HASH{rcondv}: reciprocal condition numbers for selected deflating subspaces if requested
HASH{n} : Number of eigenvalues selected if select_func is defined.

my $a = random(10,10);
my $b = random(10,10);
my ($S,$T) = mgschurx($a,$b);
sub select {
 my ($alpha,$beta) = @_;
 return $alpha->Cabs < abs($beta) ? 1 : 0;
}
my ($S, $T, $alpha, $beta, %res)  = mgschurx( $a, $b, 1, 1, 1, 1,\&select,3);

`mqr` 🎀

Computes QR decomposition. Handles complex data. Uses geqrf and orgqr or cgeqrf and cungqr from LAPACK and returns "Q" in scalar context. Works on transposed array.

(PDL(Q), PDL(R), PDL(info)) = mqr(PDL, SCALAR)
SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

my $a = random(10,10);
my ( $q, $r )  = mqr($a);
# Can compute full decomposition if nrow > ncol
$a = random(5,7);
( $q, $r )  = $a->mqr(1);

`mrq` 🍓

Computes RQ decomposition. Handles complex data. Uses gerqf and orgrq or cgerqf and cungrq from LAPACK and returns "Q" in scalar context. Works on transposed array.

(PDL(R), PDL(Q), PDL(info)) = mrq(PDL, SCALAR)
SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

my $a = random(10,10);
my ( $r, $q )  = mrq($a);
# Can compute full decomposition if nrow < ncol
$a = random(5,7);
( $r, $q )  = $a->mrq(1);

`mql` 🧚

Computes QL decomposition. Handles complex data. Uses geqlf and orgql or cgeqlf and cungql from LAPACK and returns "Q" in scalar context. Works on transposed array.

(PDL(Q), PDL(L), PDL(info)) = mql(PDL, SCALAR)
SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

my $a = random(10,10);
my ( $q, $l )  = mql($a);
# Can compute full decomposition if nrow > ncol
$a = random(5,7);
( $q, $l )  = $a->mql(1);

`mlq` ⭐

Computes LQ decomposition. Handles complex data. Uses gelqf and orglq or cgelqf and cunglq from LAPACK and returns "Q" in scalar context. Works on transposed array.

( PDL(L), PDL(Q), PDL(info) ) = mlq(PDL, SCALAR)
SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

my $a = random(10,10);
my ( $l, $q )  = mlq($a);
# Can compute full decomposition if nrow < ncol
$a = random(5,7);
( $l, $q )  = $a->mlq(1);

`msolve` 🔥

Solves linear system of equations using LU decomposition.

 A * X = B

Returns X in scalar context else X, LU, pivot vector and info. B is overwritten by X if its inplace flag is set. Supports broadcasting. Uses gesv or cgesv from LAPACK. Works on transposed arrays.

bestie heads up: B gets overwritten by X if its inplace flag is set, so don't come crying when your original B is gone 🙈

(PDL(X), (PDL(LU), PDL(pivot), PDL(info))) = msolve(PDL(A), PDL(B) )

my $a = random(5,5);
my $b = random(10,5);
my $X = msolve($a, $b);

`msolvex` 👑

Solves linear system of equations using LU decomposition.

 A * X = B

Can optionally equilibrate the matrix. Uses gesvx or cgesvx from LAPACK. Works on transposed arrays.

(PDL, (HASH(result))) = msolvex(PDL(A), PDL(B), HASH(options))

where the options are, bestie:

transpose — solves A' * X = B. 0 = false, 1 = true.
equilibrate — equilibrates A if necessary. The form of the equilibration gets handed back to you in HASH{'equilibration'}: 0 = no equilibration, 1 = row equilibration, 2 = column equilibration. Row scale factors show up in HASH{'row'}, column scale factors in HASH{'column'}. 0 = false, 1 = true.
LU — returns the lu decomposition in HASH{LU}. 0 = false, 1 = true.
A — returns scaled A if equilibration was done, in HASH{A}. 0 = false, 1 = true.
B — returns scaled B if equilibration was done, in HASH{B}. 0 = false, 1 = true.

Returned values 💖:

X (SCALAR CONTEXT)
HASH{'pivot'}: pivot indice from LU factorization
HASH{'rcondition'}: reciprocal condition of the matrix
HASH{'ferror'}: forward error bound
HASH{'berror'}: componentwise relative backward error
HASH{'rpvgrw'}: reciprocal pivot growth factor
HASH{'info'}: info, output from gesvx

my $a = random(10,10);
my $b = random(5,10);
my %options = (
   LU=>1,
   equilibrate => 1,
   );
my( $X, %result) = msolvex($a,$b,%options);

💅 mtrisolve

Solves a linear system of equations with a triangular matrix A. It's giving:

A * X = B  or A' * X = B

B gets overwritten by X if its inplace flag is set. Supports broadcasting. Uses trtrs or ctrtrs from LAPACK. Works on transposed array(s).

(PDL(X), (PDL(info)) = mtrisolve(PDL(A), SCALAR(uplo), PDL(B), SCALAR(trans), SCALAR(diag))
uplo  : UPPER  = 0 | LOWER = 1
trans  : NOTRANSPOSE  = 0 | TRANSPOSE = 1, default = 0
uplo  : UNITARY DIAGONAL = 1, default = 0

# Assume $a is upper triagonal
my $a = random(5,5);
my $b = random(5,10);
my $X = mtrisolve($a, 0, $b);

🦋 msymsolve

Solves a linear system of equations using the diagonal pivoting method with a symmetric matrix A. Literally:

A * X = B

Returns X in scalar context, else X, the block diagonal matrix D (and the multipliers), the pivot vector and info. B gets overwritten by X if its inplace flag is set. Supports broadcasting. Uses sysv or csysv from LAPACK. Works on transposed array(s). No cap.

(PDL(X), ( PDL(D), PDL(pivot), PDL(info) ) ) = msymsolve(PDL(A), SCALAR(uplo), PDL(B) )
uplo : UPPER  = 0 | LOWER = 1, default = 0

# Assume $a is symmetric
my $a = random(5,5);
my $b = random(5,10);
my $X = msymsolve($a, 0, $b);

🌸 msymsolvex

Also solves a linear system of equations using the diagonal pivoting method with a symmetric matrix A. Same energy:

A * X = B

Uses sysvx or csysvx from LAPACK. Works on transposed array.

(PDL, (HASH(result))) = msymsolvex(PDL(A), SCALAR (uplo), PDL(B), SCALAR(d))
uplo : UPPER  = 0 | LOWER = 1, default = 0
d    : whether return diagonal matrix d and pivot vector
 FALSE  = 0 | TRUE = 1, default = 0

Returned values ✨:

X (SCALAR CONTEXT)
HASH{'D'}: block diagonal matrix D (and the multipliers) (if requested)
HASH{'pivot'}: pivot indice from LU factorization (if requested)
HASH{'rcondition'}: reciprocal condition of the matrix
HASH{'ferror'}: forward error bound
HASH{'berror'}: componentwise relative backward error
HASH{'info'}: info, output from sysvx

# Assume $a is symmetric
my $a = random(10,10);
my $b = random(5,10);
my ($X, %result) = msolvex($a, 0, $b);

🎀 mpossolve

Solves a linear system of equations using Cholesky decomposition with a symmetric positive definite matrix A. Obsessed:

A * X = B

Returns X in scalar context, else X, U or L and info. B gets overwritten by X if its inplace flag is set. Supports broadcasting. Uses posv or cposv from LAPACK. Works on transposed array(s).

(PDL, (PDL, PDL, PDL)) = mpossolve(PDL(A), SCALAR(uplo), PDL(B) )
uplo : UPPER  = 0 | LOWER = 1, default = 0

# asume $a is symmetric positive definite
my $a = random(5,5);
my $b = random(5,10);
my $X = mpossolve($a, 0, $b);

🧚 mpossolvex

Solves a linear system of equations using Cholesky decomposition with a symmetric positive definite matrix A:

A * X = B

Can optionally equilibrate the matrix. Uses posvx or cposvx from LAPACK. Works on transposed array(s).

(PDL, (HASH(result))) = mpossolvex(PDL(A), SCARA(uplo), PDL(B), HASH(options))
uplo : UPPER  = 0 | LOWER = 1, default = 0

where the options are:

equilibrate — equilibrates A if necessary. The form of equilibration comes back in HASH{'equilibration'}: 0 = no equilibration, 1 = equilibration. Scale factors come back in HASH{'scale'}. 0 = false, 1 = true.
U|L — returns Cholesky factorization in HASH{U} or HASH{L}. 0 = false, 1 = true.
A — returns scaled A if equilibration was done, in HASH{A}. 0 = false, 1 = true.
B — returns scaled B if equilibration was done, in HASH{B}. 0 = false, 1 = true.

Returned values 🫶:

X (SCALAR CONTEXT)
HASH{'rcondition'}: reciprocal condition of the matrix
HASH{'ferror'}: forward error bound
HASH{'berror'}: componentwise relative backward error
HASH{'info'}: info, output from gesvx

# Assume $a is symmetric positive definite
my $a = random(10,10);
my $b = random(5,10);
my %options = (U=>1,
   equilibrate => 1,
   );
my ($X, %result) = msolvex($a, 0, $b,%opt);

💜 mlls

Solves overdetermined or underdetermined real linear systems using QR or LQ factorization.

If M > N in the M-by-N matrix A, it returns the residual sum of squares too. Uses gels or cgels from LAPACK. Works on transposed arrays.

PDL(X) = mlls(PDL(A), PDL(B), SCALAR(trans))
trans : NOTRANSPOSE  = 0 | TRANSPOSE/CONJUGATE = 1, default = 0

$a = random(4,5);
$b = random(3,5);
($x, $res) = mlls($a, $b);

⭐ mllsy

Computes the minimum-norm solution to a real linear least squares problem using a complete orthogonal factorization.

Uses gelsy or cgelsy from LAPACK. Works on transposed arrays.

( PDL(X), ( HASH(result) ) ) = mllsy(PDL(A), PDL(B))

Returned values 🌈:

X (SCALAR CONTEXT)
HASH{'A'}: complete orthogonal factorization of A
HASH{'jpvt'}: details of columns interchanges
HASH{'rank'}: effective rank of A

my $a = random(10,10);
my $b = random(10,10);
$X = mllsy($a, $b);

💀 mllss

Computes the minimum-norm solution to a real linear least squares problem using a singular value decomposition.

Uses gelss or gelsd from LAPACK. Works on transposed arrays.

( PDL(X), ( HASH(result) ) )= mllss(PDL(A), PDL(B), SCALAR(method))
method: specifies which method to use (see LAPACK for further details)
 '(c)gelss' or '(c)gelsd', default = '(c)gelsd'

Returned values 🔥:

X (SCALAR CONTEXT)
HASH{'V'}: if method = (c)gelss, the right singular vectors, stored columnwise
HASH{'s'}: singular values from SVD
HASH{'res'}: if A has full rank, the residual sum-of-squares for the solution
HASH{'rank'}: effective rank of A
HASH{'info'}: info output from method

my $a = random(10,10);
my $b = random(10,10);
$X = mllss($a, $b);

👑 mglm

Solves a general Gauss-Markov Linear Model (GLM) problem. Supports broadcasting. Uses ggglm or cggglm from LAPACK. Works on transposed arrays.

(PDL(x), PDL(y)) = mglm(PDL(a), PDL(b), PDL(d))

where d is the left hand side of the GLM equation.

my $a = random(8,10);
my $b = random(7,10);
my $d = random(10);
my ($x, $y) = mglm($a, $b, $d);

🍓 mlse

Solves a linear equality-constrained least squares (LSE) problem. Uses gglse or cgglse from LAPACK. Works on transposed arrays.

(PDL(x), PDL(res2)) = mlse(PDL(a), PDL(b), PDL(c), PDL(d))

where

c: the right hand side vector for the least squares part of the LSE problem.
d: the right hand side vector for the constrained equation.
x: the solution of the LSE problem.
res2: the residual sum of squares for the solution (returned only in array context)

my $a = random(5,4);
my $b = random(5,3);
my $c = random(4);
my $d = random(3);
my ($x, $res2) = mlse($a, $b, $c, $d);

🦋 meigen

Computes eigenvalues and, optionally, the left and/or right eigenvectors of a general square matrix (spectral decomposition). Eigenvectors are normalized (Euclidean norm = 1) and largest component real. The eigenvalues and eigenvectors returned are complex ndarrays. If only eigenvalues are requested, info is returned in array context. Supports broadcasting. Uses geev or cgeev from LAPACK. Works on transposed arrays.

(PDL(values), (PDL(LV),  (PDL(RV)), (PDL(info)))) = meigen(PDL, SCALAR(left vector), SCALAR(right vector))
left vector  : FALSE = 0 | TRUE = 1, default = 0
right vector : FALSE = 0 | TRUE = 1, default = 0

my $a = random(10,10);
my ( $eigenvalues, $left_eigenvectors, $right_eigenvectors )  = meigen($a,1,1);

💖 meigenx

Computes eigenvalues, one-norm and, optionally, the left and/or right eigenvectors of a general square matrix (spectral decomposition). Eigenvectors are normalized (Euclidean norm = 1) and largest component real. The eigenvalues and eigenvectors returned are complex ndarrays. Uses geevx or cgeevx from LAPACK. Works on transposed arrays.

(PDL(value), (PDL(lv),  (PDL(rv)), HASH(result)), HASH(result)) = meigenx(PDL, HASH(options))

where the options are, and yes we're doing all of them 💅

vector — which eigenvectors to compute:
- 'left' — computes left eigenvectors
- 'right' — computes right eigenvectors
- 'all' — computes left and right eigenvectors
- 0 — doesn't compute (default)
rcondition — reciprocal condition numbers to compute (these get returned in HASH{'rconde'} for eigenvalues and HASH{'rcondv'} for eigenvectors):
- 'value' — computes reciprocal condition numbers for eigenvalues
- 'vector' — computes reciprocal condition numbers for eigenvectors
- 'all' — computes reciprocal condition numbers for eigenvalues and eigenvectors
- 0 — doesn't compute (default)
error — specifies whether or not it computes the error bounds (returned in HASH{'eerror'} and HASH{'verror'}):

error bound = EPS * One-norm / rcond(e|v)

(reciprocal condition numbers for eigenvalues or eigenvectors must be computed).

1 — returns error bounds
0 — not computed

scale — specifies whether or not it diagonaly scales the entry matrix (scale details returned in HASH : 'scale'):
- 1 — scales
- 0 — Doesn't scale (default)
permute — specifies whether or not it permutes row and columns (permute details returned in HASH{'balance'}):
- 1 — permutes
- 0 — Doesn't permute (default)
schur — specifies whether or not it returns the Schur form (returned in HASH{'schur'}):
- 1 — returns Schur form
- 0 — not returned

Returned values ✨

eigenvalues (SCALAR CONTEXT),
left eigenvectors if requested,
right eigenvectors if requested,
HASH{'norm'}: One-norm of the matrix
HASH{'info'}: Info: if > 0, the QR algorithm failed to compute all the eigenvalues (see syevx for further details)

ok so like, here's how you'd actually use it 🫶

my $a = random(10,10);
my %options = ( rcondition => 'all',
            vector => 'all',
            error => 1,
            scale => 1,
            permute=>1,
            schur => 1
            );
my ( $eigenvalues, $left_eigenvectors, $right_eigenvectors, %result)  = meigenx($a,%options);
print "Error bounds for eigenvalues:\n $eigenvalues\n are:\n". transpose($result{'eerror'}) unless $info;

mgeigen

Computes generalized eigenvalues and, optionally, the left and/or right generalized eigenvectors for a pair of N-by-N real nonsymmetric matrices (A,B). The alpha from ratio alpha/beta is a complex ndarray. Supports broadcasting. Uses ggev or cggev from LAPACK. Works on transposed arrays. no cap 🧚

( PDL(alpha), PDL(beta), ( PDL(LV),  (PDL(RV) ), PDL(info)) = mgeigen(PDL(A),PDL(B) SCALAR(left vector), SCALAR(right vector))

left vector : FALSE = 0 | TRUE = 1, default = 0
right vector : FALSE = 0 | TRUE = 1, default = 0

like, in real life:

my $a = random(10,10);
my $b = random(10,10);
my ( $alpha, $beta, $left_eigenvectors, $right_eigenvectors )  = mgeigen($a, $b,1, 1);

mgeigenx

Computes generalized eigenvalues, one-norms and, optionally, the left and/or right generalized eigenvectors for a pair of N-by-N real nonsymmetric matrices (A,B). The alpha from ratio alpha/beta is a complex ndarray. Uses ggevx or cggevx from LAPACK. Works on transposed arrays. it's giving ✨determined eigenvalues✨

(PDL(alpha), PDL(beta), PDL(lv),  PDL(rv), HASH(result) ) = mgeigenx(PDL(a), PDL(b), HASH(options))

where options are:

vector — eigenvectors to compute:
- 'left' — computes left eigenvectors
- 'right' — computes right eigenvectors
- 'all' — computes left and right eigenvectors
- 0 — doesn't compute (default)
rcondition — reciprocal condition numbers to compute (returned in HASH{'rconde'} for eigenvalues and HASH{'rcondv'} for eigenvectors):
- 'value' — computes reciprocal condition numbers for eigenvalues
- 'vector' — computes reciprocal condition numbers for eigenvectors
- 'all' — computes reciprocal condition numbers for eigenvalues and eigenvectors
- 0 — doesn't compute (default)
error — specifies whether or not it computes the error bounds (returned in HASH{'eerror'} and HASH{'verror'}):

error bound = EPS * sqrt(one-norm(a)**2 + one-norm(b)**2) / rcond(e|v)

(reciprocal condition numbers for eigenvalues or eigenvectors must be computed).

1 — returns error bounds
0 — not computed

scale — specifies whether or not it diagonaly scales the entry matrix (scale details returned in HASH : 'lscale' and 'rscale'):
- 1 — scales
- 0 — doesn't scale (default)
permute — specifies whether or not it permutes row and columns (permute details returned in HASH{'balance'}):
- 1 — permutes
- 0 — Doesn't permute (default)
schur — specifies whether or not it returns the Schur forms (returned in HASH{'aschur'} and HASH{'bschur'}) (right or left eigenvectors must be computed.):
- 1 — returns Schur forms
- 0 — not returned

Returned values 💖

alpha,
beta,
left eigenvectors if requested,
right eigenvectors if requested,
HASH{'anorm'}, HASH{'bnorm'}: One-norm of the matrix A and B
HASH{'info'}: Info: if > 0, the QR algorithm failed to compute all the eigenvalues (see syevx for further details)

and the example, untouched, because copying is caring 🙈

$a = random(10,10);
$b = random(10,10);
%options = (rcondition => 'all',
            vector => 'all',
            error => 1,
            scale => 1,
            permute=>1,
            schur => 1
            );
($alpha, $beta, $left_eigenvectors, $right_eigenvectors, %result)  = mgeigenx($a, $b,%options);
print "Error bounds for eigenvalues:\n $eigenvalues\n are:\n". transpose($result{'eerror'}) unless $info;

msymeigen

Computes eigenvalues and, optionally eigenvectors of a real symmetric square or complex Hermitian matrix (spectral decomposition). The eigenvalues are computed from lower or upper triangular matrix. If only eigenvalues are requested, info is returned in array context. Supports broadcasting and works inplace if eigenvectors are requested. From LAPACK, uses syev or syevd for real and cheev or cheevd for complex. Works on transposed array(s). spectral decomposition?? we love a glow-up 🌈

(PDL(values), (PDL(VECTORS)), PDL(info)) = msymeigen(PDL, SCALAR(uplo), SCALAR(vector), SCALAR(method))

uplo : UPPER = 0 | LOWER = 1, default = 0
vector : FALSE = 0 | TRUE = 1, default = 0
method : 'syev' | 'syevd' | 'cheev' | 'cheevd', default = 'syevd'|'cheevd'

# Assume $a is symmetric
my $a = random(10,10);
my ( $eigenvalues, $eigenvectors )  = msymeigen($a,0,1, 'syev');

msymeigenx

Computes eigenvalues and, optionally eigenvectors of a symmetric square matrix (spectral decomposition). The eigenvalues are computed from lower or upper triangular matrix and can be selected by specifying a range. From LAPACK, uses syevx or syevr for real and cheevx or cheevr for complex. Works on transposed arrays. picky about your range? iconic, honestly 🎀

(PDL(value), (PDL(vector)), PDL(n), PDL(info), (PDL(support)) ) = msymeigenx(PDL, SCALAR(uplo), SCALAR(vector), HASH(options))

uplo : UPPER = 0 | LOWER = 1, default = 0
vector : FALSE = 0 | TRUE = 1, default = 0

where options are:

range_type — method for selecting eigenvalues:
- indice: range of indices
- interval: range of values
- 0: find all eigenvalues and optionally all vectors
range — PDL(2), lower and upper bounds interval or smallest and largest indices
- 1<=range<=N for indice
abstol — specifies error tolerance for eigenvalues
method — specifies which method to use (see LAPACK for further details):
- 'syevx' (default)
- 'syevr'
- 'cheevx' (default)
- 'cheevr'

Returned values 👑

eigenvalues (SCALAR CONTEXT),
eigenvectors if requested,
total number of eigenvalues found (n),
info
issupz or ifail (support) according to method used and returned info, for (sy|che)evx returns support only if info != 0

# Assume $a is symmetric
my $a = random(10,10);
my $overflow = lamch(9);
my $range = cat pdl(0),$overflow;
my $abstol = pdl(1.e-5);
my %options = (range_type=>'interval',
   range => $range,
   abstol => $abstol,
   method=>'syevd');
my ( $eigenvalues, $eigenvectors, $n, $isuppz )  = msymeigenx($a,0,1, %options);

msymgeigen

Computes eigenvalues and, optionally eigenvectors of a real generalized symmetric-definite or Hermitian-definite eigenproblem. The eigenvalues are computed from lower or upper triangular matrix If only eigenvalues are requested, info is returned in array context. Supports broadcasting. From LAPACK, uses sygv or sygvd for real or chegv or chegvd for complex. Works on transposed array(s). the generalized one, bestie, don't mix them up 💀

(PDL(values), (PDL(vectors)), PDL(info)) = msymgeigen(PDL(a), PDL(b),SCALAR(uplo), SCALAR(vector), SCALAR(type), SCALAR(method))

uplo : UPPER = 0 | LOWER = 1, default = 0
vector : FALSE = 0 | TRUE = 1, default = 0
type :
- 1: A * x = (lambda) * B * x
- 2: A * B * x = (lambda) * x
- 3: B * A * x = (lambda) * x
- default = 1
method : 'sygv' | 'sygvd' for real or ,'chegv' | 'chegvd' for complex, default = 'sygvd' | 'chegvd'

# Assume $a is symmetric
my $a = random(10,10);
my $b = random(10,10);
$b = $b->crossprod($b);
my ( $eigenvalues, $eigenvectors )  = msymgeigen($a, $b, 0, 1, 1, 'sygv');

msymgeigenx

Computes eigenvalues and, optionally eigenvectors of a real generalized symmetric-definite or Hermitian eigenproblem. The eigenvalues are computed from lower or upper triangular matrix and can be selected by specifying a range. Uses sygvx or cheevx from LAPACK. Works on transposed arrays. yes the "x" goes at the end, we stan a consistent naming convention 🦋

(PDL(value), (PDL(vector)), PDL(info), PDL(n), (PDL(support)) ) = msymeigenx(PDL(a), PDL(b), SCALAR(uplo), SCALAR(vector), HASH(options))

uplo : UPPER = 0 | LOWER = 1, default = 0
vector : FALSE = 0 | TRUE = 1, default = 0

where options are:

type — Specifies the problem type to be solved:
- 1: A * x = (lambda) * B * x
- 2: A * B * x = (lambda) * x
- 3: B * A * x = (lambda) * x
- default = 1
range_type — method for selecting eigenvalues:
- indice: range of indices
- interval: range of values
- 0: find all eigenvalues and optionally all vectors
range — PDL(2), lower and upper bounds interval or smallest and largest indices
- 1<=range<=N for indice
abstol — specifies error tolerance for eigenvalues

Returned values ⭐

eigenvalues (SCALAR CONTEXT),
eigenvectors if requested,
total number of eigenvalues found (n),
info
ifail according to returned info (support).

okay so the msymgeigenx example is still vibing from before — here's the rest of it, no cliffhanger bestie:

# Assume $a is symmetric
my $a = random(10,10);
my $b = random(10,10);
$b = $b->crossprod($b);
my $overflow = lamch(9);
my $range = cat pdl(0),$overflow;
my $abstol = pdl(1.e-5);
my %options = (range_type=>'interval',
   range => $range,
   abstol => $abstol,
   type => 1);
my ( $eigenvalues, $eigenvectors, $n, $isuppz )  = msymgeigenx($a, $b, 0,1, %options);

and that's the whole serve 💅 okay moving ON, we have four more decompositions to girlboss through.

💜 `mdsvd` 💜

so this one computes the SVD using Cuppen's divide and conquer algorithm 🎀 it hands back singular values in scalar context, otherwise you get left (U), singular values, right (V' — which is hermitian for complex, iykyk) singular vectors and info. Supports broadcasting ✨ and if you're only asking for the singular values, info is only returned in array context. Under the hood it's gesdd or cgesdd from LAPACK, periodt.

 (PDL(U), (PDL(s), PDL(Vt)), PDL(info)) = mdsvd(PDL, SCALAR(job))

job — 0 computes only singular values, 1 computes full SVD (square U and Vt), 2 computes SVD (singular values, right and left singular vectors). default = 1

example time:

 my $a = random(5,10);
 my ($u, $s, $vt) = mdsvd($a);
 print $u x gurney($s, $a->dims) x $vt;

🌸 `msvd` 🌸

computes SVD, iconic behaviour honestly. it can compute singular values, either U or Vt, or neither — your call bestie. returns singular values in scalar context, otherwise left (U), singular values, right (V' (hermitian for complex)) singular vector and info. Supports broadcasting 🦋 and if only the singular values are requested, info gets returned in array context. uses gesvd or cgesvd from LAPACK.

 ( (PDL(U)), PDL(s), (PDL(Vt), PDL(info)) = msvd(PDL, SCALAR(jobu), SCALAR(jobv))

jobu — 0 Doesn't compute U, 1 computes full SVD (square U), 2 computes right singular vectors. default = 1
jobv — 0 Doesn't compute Vt, 1 computes full SVD (square Vt), 2 computes left singular vectors. default = 1

another little demo for the girls:

 my $a = random(10,10);
 my ($u, $s, $vt) = msvd($a);
 print $u x gurney($s, $a->dims) x $vt;

🦋 `mgsvd` 🦋

computes the generalized (or quotient) singular value decomposition ✨ and here's the drama: if the effective rank of (A',B')' is 0, it returns only unitary V, U, Q. it handles complex data, uses ggsvd or cggsvd from LAPACK, and works on transposed arrays.

 (PDL(sa), PDL(sb), %ret) = mgsvd(PDL(a), PDL(b), %HASH(options))

where the options are:

V — whether or not to compute V (boolean, returned in HASH{'V'})
U — whether or not to compute U (boolean, returned in HASH{'U'})
Q — whether or not to compute Q (boolean, returned in HASH{'Q'})
D1 — whether or not to compute D1 (boolean, returned in HASH{'D1'})
D2 — whether or not to compute D2 (boolean, returned in HASH{'D2'})
0R — whether or not to compute 0R (boolean, returned in HASH{'0R'})
R — whether or not to compute R (boolean, returned in HASH{'R'})
X — whether or not to compute X (boolean, returned in HASH{'X'})
all — whether or not to compute all of the above. yes she said "all of it" 💅

and what you actually get back:

sa, sb — singular value pairs of A and B (generalized singular values = sa/sb)
$ret{'rank'} — effective numerical rank of (A',B')'
$ret{'info'} — info from (c)ggsvd

vibe check on the example:

 my $a = random(5,5);
 my $b = random(5,7);
 my ($c, $s, %ret) = mgsvd($a, $b, X => 1);

👑 AUTHOR 👑

This library is free software; you can redistribute it and/or modify it under the terms of the Perl Artistic License, as in the file Artistic_2 in this distribution.

and that's the whole page babes, we made it 😭💖 this was PDL::L...lgebra(3pm) in the perl v5.42.3 distribution, dated 2026-08-21 — screenshot it, send it to the group chat, and go do some linear algebra, you absolute main character 🌈👑✨

📦 where this page comes from

ubuntu stonking: libpdl-linearalgebra-perl 0.437-1build1 · Artistic-2
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-11-03: a language model's rewrite of the original, under its terms

PDL::L...lgebra(3pm)

💖 NAME 💖

🦋 SYNOPSIS 🦋

🌸 DESCRIPTION 🌸

🔥 FUNCTIONS 🔥

✨ setlaerror ✨

🎀 getlaerror 🎀

💜 issym 💜

🌈 diag 🌈

🫶 positivise 🫶

⭐ mcrossprod ⭐

💫 mrank 💫

🌟 mnorm 🌟

🎯 mdet 🎯

🍓 mposdet 🍓

🌺 mcond 🌺

💀 mrcond 💀

🦄 morth 🦄

🧚 mnull 🧚

👑 minv 👑

🌈 mtriinv 🌈

💖 msyminv 💖

🎀 mposinv 🎀

✨ mpinv ✨

🧚 mlu 🧚

🌸 mchol 🌸

⭐ mhessen ⭐

💫 mschur 💫

🦋 mschurx 🦋

mschurx 🌸 (the rest of it, bestie)

mgschur 🦋

mgschurx 💅

mqr 🎀

mrq 🍓

mql 🧚

mlq ⭐

msolve 🔥

msolvex 👑

💅 mtrisolve

🦋 msymsolve

🌸 msymsolvex

🎀 mpossolve

🧚 mpossolvex

💜 mlls

⭐ mllsy

💀 mllss

👑 mglm

🍓 mlse

🦋 meigen

💖 meigenx

mgeigen

mgeigenx

msymeigen

msymeigenx

msymgeigen

msymgeigenx

💜 mdsvd 💜

🌸 msvd 🌸

🦋 mgsvd 🦋

👑 AUTHOR 👑

`mschurx` 🌸 (the rest of it, bestie)

`mgschur` 🦋

`mgschurx` 💅

`mqr` 🎀

`mrq` 🍓

`mql` 🧚

`mlq` ⭐

`msolve` 🔥

`msolvex` 👑

💜 `mdsvd` 💜

🌸 `msvd` 🌸

🦋 `mgsvd` 🦋