Real(3pm)

📦

version 0.38 · 3 versions

0.38 · ubuntu noble (latest)
0.26 · ubuntu jammy
0.35 · debian bookworm

💖 Real(3pm) — PDL::LinearAlgebra::Real 💖

User Contributed Perl Documentation — yes bestie, we're in the contributed docs section, come sit with me 🫶

📛 NAME

PDL::LinearAlgebra::Real — it's the PDL interface to the real lapack linear algebra programming library. That's the whole vibe. Linear algebra. Real numbers. Lapack. Slay. ✨

🎀 SYNOPSIS

okay so first things first, you gotta load it:

use PDL::LinearAlgebra::Real;

$a = random (100,100);
$s = zeroes(100);
$u = zeroes(100,100);
$v = zeroes(100,100);
$info = 0;
$job = 0;
gesdd($a, $job, $info, $s , $u, $v);

look at that. random to make a matrix, zeroes to make the outputs, then gesdd does the thing and drops the answers into $s, $u, $v. main character energy fr fr 💅

🌸 DESCRIPTION

This module is your lil doorway into parts of the real lapack library. Like, it's an interface. That's its whole job.

And here's the cute part: these routines are flexible, they accept either float or double ndarrays. So you don't gotta fight it, just feed it whichever precision you're feeling. 🌈

🔥 FUNCTIONS

🦋 gtsv

Signature: ([phys]DL(n); [phys]D(n); [phys]DU(n); [io,phys]B(n,nrhs); int [o,phys]info())

So what does it do? It solves this equation right here:

A * X = B

...where A is an "n" by "n" tridiagonal matrix, solved by Gaussian elimination with partial pivoting, and B is an "n" by "nrhs" matrix. Iconic, honestly.

Pro tip bestie: you can solve A**T*X = B too — just swap the order of the arguments DU and DL. That's it. That's the trick. 🤯

NB — this differs from the LINPACK function dgtsl: here DL starts from its first element, whereas the LINPACK equivalent starts from its second element. Not me catching that off-by-one for you. 🙈

📝 Arguments

DL — On entry, DL must contain the (n-1) sub-diagonal elements of A. On exit, DL is overwritten by the (n-2) elements of the second super-diagonal of the upper triangular matrix U from the LU factorization of A, sitting in DL(1), ..., DL(n-2).
D — On entry, D must contain the diagonal elements of A. On exit, D gets overwritten by the n diagonal elements of U.
DU — On entry, DU must contain the (n-1) super-diagonal elements of A. On exit, DU is overwritten by the (n-1) elements of the first super-diagonal of the U.
B — On entry, the n by nrhs matrix of right hand side matrix B. On exit, if info = 0, the n by nrhs solution matrix X. 💖
info — the receipts, always the receipts:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value 😭
- > 0: if info = i, U(i,i) is exactly zero, and the solution has not been computed. The factorization has not been completed unless i = n.

$dl = random(float, 9);
$d = random(float, 10);
$du = random(float, 9);
$b = random(10,5);
gtsv($dl, $d, $du, $b, ($info=null));
print "X is:\n$b" unless $info;

And the bad-value tea ☕: gtsv ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. iykyk.

🧚 gesvd

Signature: ([io]A(m,n); int jobu(); int jobvt(); [o]s(minmn); [o]U(p,p); [o]VT(s,s); int [o]info())

Computes the singular value decomposition (SVD) of a real M-by-N matrix A. The moment. The icon. 🎀

The SVD is written like this:

A = U * SIGMA * V'

where SIGMA is an M-by-N matrix that's zero everywhere except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A — they're real and non-negative, and they come back in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A.

⚠️ heads up: the routine returns VT = V', not V. Don't be fooled bestie, the transpose is the whole point.

🎯 Arguments

jobu — Specifies options for computing all or part of the matrix U:
- = 0: no columns of U (no left singular vectors) are computed.
- = 1: all M columns of U are returned in array U:
- = 2: the first min(m,n) columns of U (the left singular vectors) are returned in the array U;
- = 3: the first min(m,n) columns of U (the left singular vectors) are overwritten on the array A;
jobvt — Specifies options for computing all or part of the matrix V':
- = 0: no rows of V' (no right singular vectors) are computed.
- = 1: all N rows of V' are returned in the array VT;
- = 2: the first min(m,n) rows of V' (the right singular vectors) are returned in the array VT;
- = 3: the first min(m,n) rows of V' (the right singular vectors) are overwritten on the array A;

🚨 WARNING: jobvt and jobu cannot both be 3. You cannot do both. Pick one. Periodt.

A — On entry, the M-by-N matrix A. On exit:
- if jobu = 3, A is overwritten with the first min(m,n) columns of U (the left singular vectors, stored columnwise);
- if jobvt = 3, A is overwritten with the first min(m,n) rows of V' (the right singular vectors, stored rowwise);
- if jobu != 3 and jobvt != 3, the contents of A are destroyed. 💀
s — The singular values of A, sorted so that s(i) >= s(i+1).
U — If jobu = 1, U contains the M-by-M orthogonal matrix U; if jobu = 3, U contains the first min(m,n) columns of U (the left singular vectors, stored columnwise); if jobu = 0 or 3, U is not referenced. Min size = [1,1].
VT — If jobvt = 1, VT contains the N-by-N orthogonal matrix V'; if jobvt = 2, VT contains the first min(m,n) rows of V' (the right singular vectors, stored rowwise); if jobvt = 0 or 3, VT is not referenced. Min size = [1,1].
info —
- = 0: successful exit.
- < 0: if info = -i, the i-th argument had an illegal value.
- > 0: if bdsqr did not converge, info specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero.

$a = random (float, 100,100);
$s = zeroes(float, 100);
$u = zeroes(float, 100,100);
$vt = zeroes(float, 100,100);
$info = pdl(long, 0);
gesvd($a, 2, 2, $s , $u, $vt, $info);

Bad-value behaviour, same as always: gesvd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🌸

💅 gesdd

Signature: ([io]A(m,n); int jobz(); [o]s(minmn); [o]U(p,p); [o]VT(s,s); int [o]info(); int [t]iwork(iworkn))

Also computes the singular value decomposition (SVD) of a real M-by-N matrix A — but this one is built different.

This routine use the Coppen's divide and conquer algorithm. It is much faster than the simple driver for large matrices, but uses more workspace. Speed demon behaviour, but it wants its space. 🏃‍♀️💨

🎯 Arguments

jobz — Specifies options for computing all or part of matrix:
- = 0: no columns of U or rows of V' are computed;
- = 1: all M columns of U and all N rows of V' are returned in the arrays U and VT;
- = 2: the first min(M,N) columns of U and the first min(M,N) rows of V' are returned in the arrays U and VT;
- = 3: If M >= N, the first N columns of U are overwritten on the array A and all rows of V' are returned in the array VT; otherwise, all columns of U are returned in the array U and the first M rows of V' are overwritten on the array A.
A — On entry, the M-by-N matrix A. On exit:
- if jobz = 3, A is overwritten with the first N columns of U (the left singular vectors, stored columnwise) if M >= N; A is overwritten with the first M rows of V' (the right singular vectors, stored rowwise) otherwise.
- if jobz != 3, the contents of A are destroyed. 💀
s — The singular values of A, sorted so that s(i) >= s(i+1).
U — If jobz = 1 or jobz = 3 and M < N, U contains the M-by-M orthogonal matrix U; if jobz = 2, U contains the first min(M,N) columns of U (the left singular vectors, stored columnwise); if jobz = 3 and M >= N, or jobz = 0, U is not referenced. Min size = [1,1].
VT — If jobz = 1 or jobz = 3 and M >= N, VT contains the N-by-N orthogonal matrix V'; if jobz = 2, VT contains the first min(M,N) rows of V' (the right singular vectors, stored rowwise); if jobz = 3 and M < N, or jobz = 0, VT is not referenced. Min size = [1,1].
info —
- = 0: successful exit.
- < 0: if info = -i, the i-th argument had an illegal value.
- > 0: bdsdc did not converge, updating process failed. 😭

$lines = 50;
$columns = 100;
$a = random (float, $lines, $columns);
$min = $lines < $columns ? $lines : $columns;
$s = zeroes(float, $min);
$u = zeroes(float, $lines, $lines);
$vt = zeroes(float, $columns, $columns);
$info = long (0);
gesdd($a, 1, $s , $u, $vt, $info);

Bad-value behaviour, one more time for the people in the back: gesdd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. ✨

🌈 ggsvd

Signature: ([io]A(m,n); int jobu(); int jobv(); int jobq(); [io]B(p,n); int [o]k(); int [o]l();[o]alpha(n);[o]beta(n); [o]U(q,q); [o]V(r,r); [o]Q(s,s); int [o]iwork(n); int [o]info())

Computes the generalized singular value decomposition (GSVD) of an M-by-N real matrix A and a P-by-N real matrix B:

U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )

where U, V and Q are orthogonal matrices, and Z' is the transpose
of Z.

Let K+L = the effective numerical rank of the matrix (A',B')', then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the following structures, respectively:

If M-K-L >= 0,

        K  L
D1 =     K ( I  0 )
        L ( 0  C )
    M-K-L ( 0  0 )

      K  L
D2 =   L ( 0  S )
    P-L ( 0  0 )

        N-K-L  K    L
( 0 R ) = K (  0   R11  R12 )
          L (  0    0   R22 )

where

C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
S = diag( BETA(K+1),  ... , BETA(K+L) ),
C**2 + S**2 = I.

R is stored in A(1:K+L,N-K-L+1:N) on exit.

If M-K-L < 0,

              K M-K K+L-M
D1 =   K ( I  0    0   )
    M-K ( 0  C    0   )

                K M-K K+L-M
D2 =   M-K ( 0  S    0  )
    K+L-M ( 0  0    I  )
      P-L ( 0  0    0  )

okay so this is the shape of the decomposition, bestie — it's giving triangular realness:

                              N-K-L  K   M-K  K+L-M
             ( 0 R ) =     K ( 0    R11  R12  R13  )
                         M-K ( 0     0   R22  R23  )
                       K+L-M ( 0     0    0   R33  )

where

  C = diag( ALPHA(K+1), ... , ALPHA(M) ),
  S = diag( BETA(K+1),  ... , BETA(M) ),
  C**2 + S**2 = I.

(R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
( 0  R22 R23 )
in B(M-K+1:L,N+M-K-L+1:N) on exit.

This routine is out here computing C, S, R, and optionally the orthogonal transformation matrices U, V and Q.

In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A and B implicitly gives the SVD of A*inv(B):

A*inv(B) = U*(D1*inv(D2))*V'.

If ( A',B')' has orthonormal columns, then the GSVD of A and B is also equal to the CS decomposition of A and B. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem:

A'*A x = lambda* B'*B x.

In some literature, the GSVD of A and B is presented in the form

U'*A*X = ( 0 D1 ),   V'*B*X = ( 0 D2 )
where U and V are orthogonal and X is nonsingular, D1 and D2 are "diagonal".

The former GSVD form can be converted to the latter form by taking the nonsingular matrix X as

X = Q*( I   0    )
      ( 0 inv(R) ).

Arguments ✨

jobu — = 0: U is not computed. = 1: Orthogonal matrix U is computed;
jobv — = 0: V is not computed. = 1: Orthogonal matrix V is computed;
jobq — = 0: Q is not computed. = 1: Orthogonal matrix Q is computed;
k:
l — On exit, k and l specify the dimension of the subblocks described in the Purpose section. k + l = effective numerical rank of (A',B')'.
A — On entry, the M-by-N matrix A. On exit, A contains the triangular matrix R, or part of R.
B — On entry, the P-by-N matrix B. On exit, B contains the triangular matrix R if M-k-l < 0.
alpha:
beta — On exit, alpha and beta contain the generalized singular value pairs of A and B;
- alpha(1:k) = 1,
- beta(1:k) = 0,
- and if M-k-l >= 0,
  - alpha(k+1:k+l) = C,
  - beta(k+1:k+l) = S,
- or if M-k-l < 0,
  - alpha(k+1:M)=C, alpha(M+1:k+l)=0
  - beta(k+1:M) =S, beta(M+1:k+l) =1
- and
  - alpha(k+l+1:N) = 0
  - beta(k+l+1:N) = 0
U — If jobu = 1, U contains the M-by-M orthogonal matrix U. If jobu = 0, U is not referenced. Need a minimum array of (1,1) if jobu = 0;
V — If jobv = 1, V contains the P-by-P orthogonal matrix V. If jobv = 0, V is not referenced. Need a minimum array of (1,1) if jobv = 0;
Q — If jobq = 1, Q contains the N-by-N orthogonal matrix Q. If jobq = 0, Q is not referenced. Need a minimum array of (1,1) if jobq = 0;
iwork — On exit, iwork stores the sorting information. More precisely, the following loop will sort alpha

for I = k+1, min(M,k+l)
    swap alpha(I) and alpha(iwork(I))
endfor

such that alpha(1) >= alpha(2) >= ... >= alpha(N).

info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. > 0: if info = 1, the Jacobi-type procedure failed to converge. For further details, see subroutine tgsja.

        $k = null;
        $l = null;
        $A = random(5,6);
        $B = random(7,6);
        $alpha = zeroes(6);
        $beta = zeroes(6);
        $U = zeroes(5,5);
        $V = zeroes(7,7);
        $Q = zeroes(6,6);
        $iwork = zeroes(long, 6);
        $info = null;
        ggsvd($A,1,1,1,$B,$k,$l,$alpha, $beta,$U, $V, $Q, $iwork,$info);

ggsvd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

geev 🔥

(A(n,n); int jobvl(); int jobvr(); [o]wr(n); [o]wi(n); [o]vl(m,m); [o]vr(p,p); int [o]info())

Computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors.

The right eigenvector v(j) of A satisfies:

A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue.

The left eigenvector u(j) of A satisfies:

u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate
transpose of u(j).

The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real.

Arguments ✨

jobvl — = 0: left eigenvectors of A are not computed; = 1: left eigenvectors of A are computed.
jobvr — = 0: right eigenvectors of A are not computed; = 1: right eigenvectors of A are computed.
A — A is overwritten.
wr:
wi — wr and wi contain the real and imaginary parts, respectively, of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first.
vl — If jobvl = 1, the left eigenvectors u(j) are stored one after another in the columns of vl, in the same order as their eigenvalues else vl is not referenced. If the j-th eigenvalue is real, then u(j) = vl(:,j), the j-th column of vl. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and u(j+1) = vl(:,j) - i*vl(:,j+1). Min size = [1].
vr — If jobvr = 1, the right eigenvectors v(j) are stored one after another in the columns of vr, in the same order as their eigenvalues else vr is not referenced. If the j-th eigenvalue is real, then v(j) = vr(:,j), the j-th column of vr. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and v(j+1) = vr(:,j) - i*vr(:,j+1). Min size = [1].
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. > 0: if info = i, the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed; elements i+1:N of wr and wi contain eigenvalues which have converged.

        $a = random (5, 5);
        $wr = zeroes(5);
        $wi = zeroes($wr);
        $vl = zeroes($a);
        $vr = zeroes($a);
        $info = null;
        geev($a, 1, 1, $wr, $wi, $vl, $vr, $info);

geev ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

geevx ✨

([io]A(n,n);  int jobvl(); int jobvr(); int balance(); int sense(); [o]wr(n); [o]wi(n); [o]vl(m,m); [o]vr(p,p); int [o]ilo(); int [o]ihi(); [o]scale(n); [o]abnrm(); [o]rconde(q); [o]rcondv(r); int [o]info(); int [t]iwork(iworkn))

Computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors.

Optionally also, it computes a balancing transformation to improve the conditioning of the eigenvalues and eigenvectors (ilo, ihi, scale, and abnrm), reciprocal condition numbers for the eigenvalues (rconde), and reciprocal condition numbers for the right eigenvectors (rcondv).

The right eigenvector v(j) of A satisfies:

A * v(j) = lambda(j) * v(j)
where lambda(j) is its eigenvalue.

The left eigenvector u(j) of A satisfies:

u(j)**H * A = lambda(j) * u(j)**H
where u(j)**H denotes the conjugate transpose of u(j).

The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real.

Balancing a matrix means permuting the rows and columns to make it more nearly upper triangular, and applying a diagonal similarity transformation D * A * D**(-1), where D is a diagonal matrix, to make its rows and columns closer in norm and the condition numbers of its eigenvalues and eigenvectors smaller. The computed reciprocal condition numbers correspond to the balanced matrix. Permuting rows and columns will not change the condition numbers (in exact arithmetic) but diagonal scaling will. For further explanation of balancing, see section 4.10.2 of the LAPACK Users' Guide.

Arguments ✨

balance — Indicates how the input matrix should be diagonally scaled and/or permuted to improve the conditioning of its eigenvalues.
- = 0: Do not diagonally scale or permute;
- = 1: Perform permutations to make the matrix more nearly upper triangular. Do not diagonally scale;
- = 2: Diagonally scale the matrix, i.e. replace A by D*A*D**(-1), where D is a diagonal matrix chosen to make the rows and columns of A more equal in norm. Do not permute;
- = 3: Both diagonally scale and permute A. Computed reciprocal condition numbers will be for the matrix after balancing and/or permuting. Permuting does not change condition numbers (in exact arithmetic), but balancing does.
jobvl — = 0: left eigenvectors of A are not computed; = 1: left eigenvectors of A are computed. If sense = 1 or 3, jobvl must = 1.
jobvr; — = 0: right eigenvectors of A are not computed; = 1: right eigenvectors of A are computed. If sense = 1 or 3, jobvr must = 1.
sense — Determines which reciprocal condition numbers are computed.
- = 0: None are computed;
- = 1: Computed for eigenvalues only;
- = 2: Computed for right eigenvectors only;
- = 3: Computed for eigenvalues and right eigenvectors. If sense = 1 or 3, both left and right eigenvectors must also be computed (jobvl = 1 and jobvr = 1).
A — The N-by-N matrix. It is overwritten. If jobvl = 1 or jobvr = 1, A contains the real Schur form of the balanced version of the input matrix A.
wr, wi — okay so these two hold the real and imaginary parts of the eigenvalues we computed, respectively 💖. Here's the tea: complex conjugate eigenvalue pairs show up back to back, and the one with the positive imaginary part comes first. Periodt.
vl — if jobvl = 1, the left eigenvectors u(j) get stored one after another in the columns of vl, in the exact same order as their eigenvalues. If not, vl isn't referenced at all, she's just sitting there 🙈.
- If the j-th eigenvalue is real, then u(j) = vl(:,j) — that's the j-th column of vl, bestie.
- If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and u(j+1) = vl(:,j) - i*vl(:,j+1) ✨.
- Min size = [1].
vr — if jobvr = 1, the right eigenvectors v(j) get stored one after another in the columns of vr, same order as their eigenvalues. Otherwise vr is not referenced, no cap 💅.
- If the j-th eigenvalue is real, then v(j) = vr(:,j), the j-th column of vr.
- If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and v(j+1) = vr(:,j) - i*vr(:,j+1).
- Min size = [1].
ilo, ihi — integer values that get determined when A was balanced 🎀. The balanced A(i,j) = 0 if I > J and J = 1,...,ilo-1 or I = ihi+1,...,N.
scale — all the deets on the permutations and scaling factors that got applied when balancing A 🌸. If P(j) is the index of the row and column interchanged with row and column j, and D(j) is the scaling factor applied to row and column j, then:

scale(J) = P(J),    for J = 1,...,ilo-1
         = D(J),    for J = ilo,...,ihi
         = P(J)     for J = ihi+1,...,N.

The order in which the interchanges are made is N to ihi+1, then 1 to ilo-1 — iconic little sequence, honestly 🦋.

abnrm — the one-norm of the balanced matrix, aka the maximum of the sum of absolute values of elements of any column. It's giving "biggest column energy" ⭐.
rconde — rconde(j) is the reciprocal condition number of the j-th eigenvalue 👑.
rcondv — rcondv(j) is the reciprocal condition number of the j-th right eigenvector 💜.
info — the vibe check at the end:

= 0:  successful exit
< 0:  if info = -i, the i-th argument had an illegal value.
> 0:  if info = i, the QR algorithm failed to compute all the
      eigenvalues, and no eigenvectors or condition numbers
      have been computed; elements 1:ilo-1 and i+1:N of wr
      and wi contain eigenvalues which have converged.

Okay and here's how you'd actually call it, in all her glory:

$a = random (5,5);
$wr = zeroes(5);
$wi = zeroes(5);
$vl = zeroes(5,5);
$vr = zeroes(5,5);
$ilo = null;
$ihi = null;
$scale  = zeroes(5);
$abnrm = null;
$rconde = zeroes(5);
$rcondv = zeroes(5);
$info = null;
geevx($a, 1,1,3,3,$wr, $wi, $vl, $vr, $ilo, $ihi, $scale, $abnrm,$rconde, $rcondv, $info);

geevx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. Not her problem, she's just passing the vibe along 😭.

ggev 🦋

Signature: (A(n,n); int [phys]jobvl();int [phys]jobvr();B(n,n);[o]alphar(n);[o]alphai(n);[o]beta(n);[o]VL(m,m);[o]VR(p,p);int [o]info())

So for a pair of N-by-N real nonsymmetric matrices (A,B), this girl computes the generalized eigenvalues, and optionally the left and/or right generalized eigenvectors 🔥.

A generalized eigenvalue for a pair of matrices (A,B) is a scalar lambda, or a ratio alpha/beta = lambda, such that A - lambda*B is singular. It's usually represented as the pair (alpha,beta), because there's a reasonable interpretation for beta=0, and even for both of them being zero 💀. (Yeah. Both zero. She said what she said.)

The right eigenvector v(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies:

A * v(j) = lambda(j) * B * v(j).

The left eigenvector u(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies:

u(j)**H * A  = lambda(j) * u(j)**H * B .

where u(j)**H is the conjugate-transpose of u(j) ✨.

Arguments 💫

jobvl — pick your fighter:

= 0:  do not compute the left generalized eigenvectors;
= 1:  compute the left generalized eigenvectors.

jobvr — same energy, other side:

= 0:  do not compute the right generalized eigenvectors;
= 1:  compute the right generalized eigenvectors.

A — on entry, the matrix A in the pair (A,B). On exit? A has been overwritten. She took it, bestie, it's gone 🫶.
B — on entry, the matrix B in the pair (A,B). On exit, B has been overwritten too. No survivors 😭.
alphar, alphai, beta — on exit, (alphar(j) + alphai(j)*i)/beta(j), for j=1,...,N, will be the generalized eigenvalues. If alphai(j) is zero, then the j-th eigenvalue is real; if it's positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with alphai(j+1) negative 👑.

Note: the quotients alphar(j)/beta(j) and alphai(j)/beta(j) may easily over- or underflow, and beta(j) may even be zero 💀. So please, the user should avoid naively computing the ratio alpha/beta. However — plot twist — alphar and alphai will always be less than and usually comparable with norm(A) in magnitude, and beta always less than and usually comparable with norm(B) 🌈.

VL — if jobvl = 1, the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1) 🌸. Each eigenvector will be scaled so the largest component has abs(real part)+abs(imag. part)=1. Not referenced if jobvl = 0.
VR — if jobvr = 1, the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1) 💜. Each eigenvector will be scaled so the largest component has abs(real part)+abs(imag. part)=1. Not referenced if jobvr = 0.
info — the final report card:

= 0:  successful exit
< 0:  if info = -i, the i-th argument had an illegal value.
= 1,...,N:
      The QZ iteration failed.  No eigenvectors have been
      calculated, but alphar(j), alphai(j), and beta(j)
      should be correct for j=info+1,...,N.
> N:  =N+1: other than QZ iteration failed in hgeqz.
      =N+2: error return from tgevc.

And the call, nice and clean:

$a = random(5,5);
$b = random(5,5);
$alphar = zeroes(5);
$alphai = zeroes(5);
$beta = zeroes(5);
$vl = zeroes(5,5);
$vr = zeroes(5,5);
ggev($a, 1, 1, $b, $alphar, $alphai, $beta, $vl, $vr, ($info=null));

ggev ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. Same rules as always, she's consistent 💅.

ggevx 🎀

Signature: ([io,phys]A(n,n);int balanc();int jobvl();int jobvr();int sense();[io,phys]B(n,n);[o]alphar(n);[o]alphai(n);[o]beta(n);[o]VL(m,m);[o]VR(p,p);int [o]ilo();int [o]ihi();[o]lscale(n);[o]rscale(n);[o]abnrm();[o]bbnrm();[o]rconde(r);[o]rcondv(s);int [o]info(); int [t]bwork(bworkn); int [t]iwork(iworkn))

Okay so for a pair of N-by-N real nonsymmetric matrices (A,B), this one computes the generalized eigenvalues, and optionally the left and/or right generalized eigenvectors 🦋.

Optionally, also, she computes a balancing transformation to improve the conditioning of the eigenvalues and eigenvectors (ilo, ihi, lscale, rscale, abnrm, and bbnrm), reciprocal condition numbers for the eigenvalues (rconde), and reciprocal condition numbers for the right eigenvectors (rcondv) ✨. She's doing the most, and honestly? She ate.

The right eigenvector v(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies:

A * v(j) = lambda(j) * B * v(j) .

The left eigenvector u(j) corresponding to the eigenvalue lambda(j) of (A,B) satisfies:

u(j)**H * A  = lambda(j) * u(j)**H * B.

where u(j)**H is the conjugate-transpose of u(j) ⭐.

Further Details 🌈

Balancing a matrix pair (A,B) includes, first, permuting rows and columns to isolate eigenvalues, and second, applying diagonal similarity transformation to the rows and columns to make the rows and columns as close in norm as possible 👑. The computed reciprocal condition numbers correspond to the balanced matrix. Permuting rows and columns will not change the condition numbers (in exact arithmetic) but diagonal scaling will. For further explanation of balancing, see section 4.11.1.2 of LAPACK Users' Guide 📖.

An approximate error bound on the chordal distance between the i-th computed generalized eigenvalue w and the corresponding exact eigenvalue lambda is:

chord(w, lambda) <= EPS * norm(abnrm, bbnrm) / rconde(I)

An approximate error bound for the angle between the i-th computed eigenvector vl(i) or vr(i) is given by:

EPS * norm(abnrm, bbnrm) / DIF(i).

For further explanation of the reciprocal condition numbers rconde and rcondv, see section 4.11 of LAPACK User's Guide 💖.

Arguments 💫

balanc — specifies the balance option to be performed:

= 0:  do not diagonally scale or permute;
= 1:  permute only;
= 2:  scale only;
= 3:  both permute and scale.

Computed reciprocal condition numbers will be for the matrices after permuting and/or balancing. Permuting does not change condition numbers (in exact arithmetic), but balancing does — important, don't skip that part 🎀.

jobvl — your pick:

= 0:  do not compute the left generalized eigenvectors;
= 1:  compute the left generalized eigenvectors.

jobvr — and the other side:

= 0:  do not compute the right generalized eigenvectors;
= 1:  compute the right generalized eigenvectors.

sense — okay so this decides which reciprocal condition numbers get computed, bestie. = 0: none are computed; = 1: computed for eigenvalues only; = 2: computed for eigenvectors only; = 3: computed for eigenvalues and eigenvectors. 💅
A — On entry, the matrix A in the pair (A,B). On exit, A has been overwritten. If jobvl=1 or jobvr=1 or both, then A contains the first part of the real Schur form of the "balanced" versions of the input A and B. ✨
B — On entry, the matrix B in the pair (A,B). On exit, B has been overwritten. If jobvl=1 or jobvr=1 or both, then B contains the second part of the real Schur form of the "balanced" versions of the input A and B. 🦋
alphar, alphai, beta — On exit, (alphar(j) + alphai(j)*i)/beta(j), j=1,...,N, will be the generalized eigenvalues. If alphai(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with alphai(j+1) negative. 🌸

Note: the quotients alphar(j)/beta(j) and alphai(j)/beta(j) may easily over- or underflow, and beta(j) may even be zero. Thus, the user should avoid naively computing the ratio ALPHA/beta. However, alphar and alphai will be always less than and usually comparable with norm(A) in magnitude, and beta always less than and usually comparable with norm(B). 💀

vl — If jobvl = 1, the left eigenvectors u(j) are stored one after another in the columns of vl, in the same order as their eigenvalues. If the j-th eigenvalue is real, then u(j) = vl(:,j), the j-th column of vl. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then u(j) = vl(:,j)+i*vl(:,j+1) and u(j+1) = vl(:,j)-i*vl(:,j+1). Each eigenvector will be scaled so the largest component have abs(real part) + abs(imag. part) = 1. Not referenced if jobvl = 0. 🎀
vr — If jobvr = 1, the right eigenvectors v(j) are stored one after another in the columns of vr, in the same order as their eigenvalues. If the j-th eigenvalue is real, then v(j) = vr(:,j), the j-th column of vr. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then v(j) = vr(:,j)+i*vr(:,j+1) and v(j+1) = vr(:,j)-i*vr(:,j+1). Each eigenvector will be scaled so the largest component have abs(real part) + abs(imag. part) = 1. Not referenced if jobvr = 0. 🌈
ilo, ihi — ilo and ihi are integer values such that on exit A(i,j) = 0 and B(i,j) = 0 if i > j and j = 1,...,ilo-1 or i = ihi+1,...,N. If balanc = 0 or 2, ilo = 1 and ihi = N. 👑
lscale — Details of the permutations and scaling factors applied to the left side of A and B. If PL(j) is the index of the row interchanged with row j, and DL(j) is the scaling factor applied to row j, then

lscale(j) = PL(j)  for j = 1,...,ilo-1
          = DL(j)  for j = ilo,...,ihi
          = PL(j)  for j = ihi+1,...,N.

The order in which the interchanges are made is N to ihi+1, then 1 to ilo-1. 💅

rscale — Details of the permutations and scaling factors applied to the right side of A and B. If PR(j) is the index of the column interchanged with column j, and DR(j) is the scaling factor applied to column j, then

rscale(j) = PR(j)  for j = 1,...,ilo-1
          = DR(j)  for j = ilo,...,ihi
          = PR(j)  for j = ihi+1,...,N

The order in which the interchanges are made is N to ihi+1, then 1 to ilo-1. ✨

abnrm — The one-norm of the balanced matrix A. 🍓
bbnrm — The one-norm of the balanced matrix B. 🧚
rconde — If sense = 1 or 3, the reciprocal condition numbers of the selected eigenvalues, stored in consecutive elements of the array. For a complex conjugate pair of eigenvalues two consecutive elements of rconde are set to the same value. Thus rconde(j), rcondv(j), and the j-th columns of vl and vr all correspond to the same eigenpair (but not in general the j-th eigenpair, unless all eigenpairs are selected). If sense = 2, rconde is not referenced. 💜
rcondv — If sense = 2 or 3, the estimated reciprocal condition numbers of the selected eigenvectors, stored in consecutive elements of the array. For a complex eigenvector two consecutive elements of rcondv are set to the same value. If the eigenvalues cannot be reordered to compute rcondv(j), rcondv(j) is set to 0; this can only occur when the true value would be very small anyway. If sense = 1, rcondv is not referenced. 🫶
info —

= 0:  successful exit
< 0:  if info = -i, the i-th argument had an illegal value.
= 1,...,N:
      The QZ iteration failed.  No eigenvectors have been
      calculated, but alphar(j), alphai(j), and beta(j)
      should be correct for j=info+1,...,N.
> N:  =N+1: other than QZ iteration failed in hgeqz.
      =N+2: error return from tgevc.

Okay now the example, no cap:

$a = random(5,5);
$b = random(5,5);
$alphar = zeroes(5);
$alphai = zeroes(5);
$beta = zeroes(5);
$vl = zeroes(5,5);
$vr = zeroes(5,5);
$lscale = zeroes(5);
$rscale = zeroes(5);
$ilo = null;
$ihi = null;
$abnrm = null;
$bbnrm = null;
$rconde = zeroes(5);
$rcondv = zeroes(5);
ggevx($a, 3, 1, 1, 3, $b, $alphar, $alphai, $beta, $vl, $vr,
$ilo, $ihi, $lscale, $rscale, $abnrm, $bbnrm, $rconde,$rcondv,($info=null));

ggevx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 😭

`gees`

Signature: ([io]A(n,n);  int jobvs(); int sort(); [o]wr(n); [o]wi(n); [o]vs(p,p); int [o]sdim(); int [o]info(); int [t]bwork(bworkn); SV* select_func)

Computes for an N-by-N real nonsymmetric matrix A, the eigenvalues, the real Schur form T, and, optionally, the matrix of Schur vectors Z. This gives the Schur factorization A = Z*T*Z'. 🔥

Optionally, it also orders the eigenvalues on the diagonal of the real Schur form so that selected eigenvalues are at the top left. The leading columns of Z then form an orthonormal basis for the invariant subspace corresponding to the selected eigenvalues. 👑

A matrix is in real Schur form if it is upper quasi-triangular with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the form

[  a  b  ]
[  c  a  ]
where b*c < 0.

The eigenvalues of such a block are a +- sqrt(bc). 🌸

Arguments

jobvs — = 0: Schur vectors are not computed; = 1: Schur vectors are computed. ✨
sort — Specifies whether or not to order the eigenvalues on the diagonal of the Schur form. = 0: Eigenvalues are not ordered; = 1: Eigenvalues are ordered (see select_func). 💅
select_func — If sort = 1, select_func is used to select eigenvalues to sort to the top left of the Schur form. If sort = 0, select_func is not referenced. An eigenvalue wr(j)+sqrt(-1)*wi(j) is selected if select_func(SCALAR(wr(j)), SCALAR(wi(j))) is true; i.e., if either one of a complex conjugate pair of eigenvalues is selected, then both complex eigenvalues are selected. Note that a selected complex eigenvalue may no longer satisfy select_func(wr(j),wi(j)) = 1 after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned); in this case info is set to N+2 (see info below). 🦋
A — The N-by-N matrix A. On exit, A has been overwritten by its real Schur form T. 🎀
sdim — If sort = 0, sdim = 0. If sort = 1, sdim = number of eigenvalues (after sorting) for which select_func is true. (Complex conjugate pairs for which select_func is true for either eigenvalue count as 2.) 💜
wr, wi — wr and wi contain the real and imaginary parts, respectively, of the computed eigenvalues in the same order that they appear on the diagonal of the output Schur form T. Complex conjugate pairs of eigenvalues will appear consecutively with the eigenvalue having the positive imaginary part first. 🌈
vs — If jobvs = 1, vs contains the orthogonal matrix Z of Schur vectors else vs is not referenced. ⭐
info —

= 0: successful exit
< 0: if info = -i, the i-th argument had an illegal value.
> 0: if info = i, and i is
   <= N: the QR algorithm failed to compute all the
         eigenvalues; elements 1:ILO-1 and i+1:N of wr and wi
         contain those eigenvalues which have converged; if
         jobvs = 1, vs contains the matrix which reduces A
         to its partially converged Schur form.
   = N+1: the eigenvalues could not be reordered because some
         eigenvalues were too close to separate (the problem
         is very ill-conditioned);
   = N+2: after reordering, roundoff changed values of some
         complex eigenvalues so that leading eigenvalues in
         the Schur form no longer satisfy select_func = 1  This
         could also be caused by underflow due to scaling.

Example time, bestie:

sub select_function{
       my ($a, $b ) = @_;
       # Stable "continuous time" eigenspace
       return $a < 0 ? 1 : 0;
}
$A = random (5,5);
$wr= zeroes(5);
$wi = zeroes(5);
$vs = zeroes(5,5);
$sdim  = null;
$info = null;
gees($A, 1,1, $wr, $wi, $vs, $sdim, $info,\&select_function);

gees ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💖

🌈 geesx 🌈

Signature: ([io]A(n,n);  int jobvs(); int sort(); int sense(); [o]wr(n); [o]wi(n); [o]vs(p,p); int [o]sdim(); [o]rconde();[o]rcondv(); int [o]info(); int [t]bwork(bworkn); SV* select_func)

okay bestie, geesx is the full-glam version of the Schur decomposition. For an N-by-N real nonsymmetric matrix A, it computes the eigenvalues, the real Schur form T, and — if you want it — the matrix of Schur vectors Z. That's the Schur factorization:

A = Z*T*Z'

and optionally it also orders the eigenvalues along the diagonal of the real Schur form so your selected ones get yeeted to the top left. On top of that it can compute a reciprocal condition number for the average of the selected eigenvalues (rconde) 💅 and a reciprocal condition number for the right invariant subspace that goes with those selected eigenvalues (rcondv). The leading columns of Z form an orthonormal basis for that invariant subspace. main character basis energy fr.

For the deeper lore on the reciprocal condition numbers rconde and rcondv, go look at Section 4.10 of the LAPACK Users' Guide — that's where these quantities are called s and sep respectively.

A real matrix is in real Schur form if it's upper quasi-triangular with 1-by-1 and 2-by-2 blocks. The 2-by-2 blocks get standardized into the form

        [  a  b  ]
        [  c  a  ]

where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). math but make it fashion. ✨

Arguments

jobvs — decides if you're getting Schur vectors at all:
- = 0: Schur vectors are not computed;
- = 1: Schur vectors are computed. 👑
sort — whether or not to order the eigenvalues on the diagonal of the Schur form:
- = 0: Eigenvalues are not ordered;
- = 1: Eigenvalues are ordered (see select_func).
select_func — if sort = 1, this is what picks which eigenvalues get sorted to the top left of the Schur form; otherwise select_func is not referenced. An eigenvalue wr(j)+sqrt(-1)*wi(j) is selected if select_func(wr(j),wi(j)) is true — i.e., if either one of a complex conjugate pair of eigenvalues is selected, then both are. Heads up though: a selected complex eigenvalue may no longer satisfy select_func(wr(j),wi(j)) = 1 after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned); in that case info may be set to N+3 (see info below). no cap, that's a real thing that can happen.
sense — which reciprocal condition numbers get 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 sense = 1, 2 or 3, sort must equal 1. 🫶

A — on entry, the N-by-N matrix A. On exit, A is overwritten by its real Schur form T. bye original A, you were cute.
sdim — if sort = 0, sdim = 0. If sort = 1, sdim = number of eigenvalues (after sorting) for which select_func is 1. (Complex conjugate pairs for which select_func is 1 for either eigenvalue count as 2.)
wr, wi — these two hold the real and imaginary parts, respectively, of the computed eigenvalues, in the same order that they appear on the diagonal of the output Schur form T. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. 💫
vs — if jobvs = 1, vs contains the orthogonal matrix Z of Schur vectors; else vs is not referenced.
rconde — if sense = 1 or 3, rconde contains the reciprocal condition number for the average of the selected eigenvalues. Not referenced if sense = 0 or 2.
rcondv — if sense = 2 or 3, rcondv contains the reciprocal condition number for the selected right invariant subspace. Not referenced if sense = 0 or 1.
info — the receipts 📋:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value.
- > 0: if info = i, and i is
  - <= N: the QR algorithm failed to compute all the eigenvalues; elements 1:ilo-1 and i+1:N of wr and wi contain those eigenvalues which have converged; if jobvs = 1, vs contains the transformation which reduces A to its partially converged Schur form.
  - = N+1: the eigenvalues could not be reordered because some eigenvalues were too close to separate (the problem is very ill-conditioned);
  - = N+2: after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Schur form no longer satisfy select_func=1. This could also be caused by underflow due to scaling.

ok so here's the vibe in Perl, bestie:

sub select_function{
       my ($a, $b) = @_;
       # Stable "discrete time" eigenspace
       return sqrt($a**2 + $b**2) < 1 ? 1 : 0;
}
$A = random (5,5);
$wr= zeroes(5);
$wi = zeroes(5);
$vs = zeroes(5,5);
$sdim  = null;
$rconde = null;
$rcondv = null;
$info = null;
geesx($A, 1,1, 3, $wr, $wi, $vs, $sdim, $rconde, $rcondv, $info, \&select_function);

geesx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. it's giving considerate, honestly.

🦋 gges 🦋

Signature: ([io]A(n,n); int jobvsl();int jobvsr();int sort();[io]B(n,n);[o]alphar(n);[o]alphai(n);[o]beta(n);[o]VSL(m,m);[o]VSR(p,p);int [o]sdim();int [o]info(); int [t]bwork(bworkn); SV* select_func)

For a pair of N-by-N real nonsymmetric matrices (A,B), gges computes the generalized eigenvalues, the generalized real Schur form (S,T), and optionally the left and/or right matrices of Schur vectors (VSL and VSR). Which gives you the generalized Schur factorization:

(A,B) = ( (VSL)*S*(VSR)', (VSL)*T*(VSR)' )

Optionally it also orders the eigenvalues so that a selected cluster of eigenvalues appears in the leading diagonal blocks of the upper quasi-triangular matrix S and the upper triangular matrix T. The leading columns of VSL and VSR then form an orthonormal basis for the corresponding left and right eigenspaces (deflating subspaces). iconic behavior. 🌈

(If you only need the generalized eigenvalues, use the driver ggev instead — it's faster. ngl, know your options.)

A generalized eigenvalue for a pair of matrices (A,B) is a scalar w or a ratio alpha/beta = w, such that A - w*B is singular. It's usually represented as the pair (alpha,beta), because there's a reasonable interpretation for beta=0 or both being zero.

A pair of matrices (S,T) is in generalized real Schur form if T is upper triangular with non-negative diagonal and S is block upper triangular with 1-by-1 and 2-by-2 blocks. The 1-by-1 blocks correspond to real generalized eigenvalues, while 2-by-2 blocks of S will be "standardized" by making the corresponding elements of T have the form:

        [  a  0  ]
        [  0  b  ]

and the pair of corresponding 2-by-2 blocks in S and T will have a complex conjugate pair of generalized eigenvalues. 💫

Arguments

jobvsl — = 0: do not compute the left Schur vectors; = 1: compute the left Schur vectors.
jobvsr — = 0: do not compute the right Schur vectors; = 1: compute the right Schur vectors.
sort — whether or not to order the eigenvalues on the diagonal of the generalized Schur form:
- = 0: Eigenvalues are not ordered;
- = 1: Eigenvalues are ordered (see delztg).
delztg — if sort = 0, delztg is not referenced. If sort = 1, delztg is used to select eigenvalues to sort to the top left of the Schur form. An eigenvalue (alphar(j)+alphai(j))/beta(j) is selected if delztg(alphar(j),alphai(j),beta(j)) is true; i.e. if either one of a complex conjugate pair of eigenvalues is selected, then both complex eigenvalues are selected.

Note that in the ill-conditioned case, a selected complex eigenvalue may no longer satisfy delztg(alphar(j),alphai(j), beta(j)) = 1 after ordering. info is to be set to N+2 in this case. yes, that's a whole warning, take it seriously.

A — on entry, the first of the pair of matrices. On exit, A has been overwritten by its generalized Schur form S.
B — on entry, the second of the pair of matrices. On exit, B has been overwritten by its generalized Schur form T.
sdim — if sort = 0, sdim = 0. If sort = 1, sdim = number of eigenvalues (after sorting) for which delztg is true. (Complex conjugate pairs for which delztg is true for either eigenvalue count as 2.)
alphar, alphai, beta — on exit, (alphar(j) + alphai(j)*i)/beta(j), j=1,...,N, will be the generalized eigenvalues. alphar(j) + alphai(j)*i, and beta(j), j=1,...,N are the diagonals of the complex Schur form (S,T) that would result if the 2-by-2 diagonal blocks of the real Schur form of (A,B) were further reduced to triangular form using 2-by-2 complex unitary transformations. If alphai(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with alphai(j+1) negative. 💅

Note: the quotients alphar(j)/beta(j) and alphai(j)/beta(j) may easily over- or underflow, and beta(j) may even be zero. Thus, the user should avoid naively computing the ratio. However, alphar and alphai will be always less than and usually comparable with norm(A) in magnitude, and beta always less than and usually comparable with norm(B). don't do the naive thing bestie, I'm begging.

VSL — if jobvsl = 1, VSL will contain the left Schur vectors. Not referenced if jobvsl = 0. The leading dimension must always be >=1.
VSR — if jobvsr = 1, VSR will contain the right Schur vectors. Not referenced if jobvsr = 0. The leading dimension must always be >=1.
info — the vibe check on how it went:
- = 0: slayed, successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value (one of your inputs was giving nothing but red flags 🚩)
- = 1,...,N: the QZ iteration failed. (A,B) are not in Schur form, but alphar(j), alphai(j), and beta(j) should still be correct for j=info+1,...,N. 😭
- > N:
  - =N+1: something other than QZ iteration failed in hgeqz.
  - =N+2: after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Generalized Schur form no longer satisfy delztg=1. This could also be caused due to scaling. 💀
  - =N+3: reordering failed in tgsen.

sub my_select{
       my ($zr, $zi, $d) = @_;
       # stable generalized eigenvalues for continuous time
       return ( ($zr < 0 && $d > 0 ) || ($zr > 0 && $d < 0) ) ?  1 : 0;
}
$a = random(5,5);
$b = random(5,5);
$sdim = null;
$alphar = zeroes(5);
$alphai = zeroes(5);
$beta = zeroes(5);
$vsl = zeroes(5,5);
$vsr = zeroes(5,5);
gges($a, 1, 1, 1, $b, $alphar, $alphai, $beta, $vsl, $vsr, $sdim,($info=null), \&my_select);

gges ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. bestie it just does not care about your bad values, periodt.

🦋 ggesx 🦋

Signature:

([io]A(n,n); int jobvsl();int jobvsr();int sort();int sense();[io]B(n,n);[o]alphar(n);[o]alphai(n);[o]beta(n);[o]VSL(m,m);[o]VSR(p,p);int [o]sdim();[o]rconde(q=2);[o]rcondv(q=2);int [o]info(); int [t]bwork(bworkn); int [t]iwork(iworkn); SV* select_func)

okay so this one computes, for a pair of N-by-N real nonsymmetric matrices (A,B), the generalized eigenvalues, the real Schur form (S,T), and, optionally, the left and/or right matrices of Schur vectors (VSL and VSR). That gives you the generalized Schur factorization:

(A,B) = ( (VSL) S (VSR)', (VSL) T (VSR)' )

Optionally, it also orders the eigenvalues so that a selected cluster of eigenvalues appears in the leading diagonal blocks of the upper quasi-triangular matrix S and the upper triangular matrix T; computes a reciprocal condition number for the average of the selected eigenvalues (RCONDE); and computes a reciprocal condition number for the right and left deflating subspaces corresponding to the selected eigenvalues (RCONDV). The leading columns of VSL and VSR then form an orthonormal basis for the corresponding left and right eigenspaces (deflating subspaces). iconic behavior honestly.

A generalized eigenvalue for a pair of matrices (A,B) is a scalar w or a ratio alpha/beta = w, such that A - w*B is singular. It is usually represented as the pair (alpha,beta), as there is a reasonable interpretation for beta=0 or for both being zero.

A pair of matrices (S,T) is in generalized real Schur form if T is upper triangular with non-negative diagonal and S is block upper triangular with 1-by-1 and 2-by-2 blocks. The 1-by-1 blocks correspond to real generalized eigenvalues, while 2-by-2 blocks of S will be "standardized" by making the corresponding elements of T have the form:

[  a  0  ]
[  0  b  ]

and the pair of corresponding 2-by-2 blocks in S and T will have a complex conjugate pair of generalized eigenvalues. it's giving structured, it's giving ✨organized✨.

Further details

An approximate (asymptotic) bound on the average absolute error of the selected eigenvalues is

EPS * norm((A, B)) / RCONDE( 1 ).

An approximate (asymptotic) bound on the maximum angular error in the computed deflating subspaces is

EPS * norm((A, B)) / RCONDV( 2 ).

See LAPACK User's Guide, section 4.11 for more information.

Arguments

jobvsl — = 0: do not compute the left Schur vectors; = 1: compute the left Schur vectors.
jobvsr — = 0: do not compute the right Schur vectors; = 1: compute the right Schur vectors.
sort — specifies whether or not to order the eigenvalues on the diagonal of the generalized Schur form. = 0: eigenvalues are not ordered; = 1: eigenvalues are ordered (see delztg).
delztg — if sort = 0, delztg is not referenced. If sort = 1, delztg is used to select eigenvalues to sort to the top left of the Schur form. An eigenvalue (alphar(j)+alphai(j))/beta(j) is selected if delztg(alphar(j),alphai(j),beta(j)) is true; i.e. if either one of a complex conjugate pair of eigenvalues is selected, then both complex eigenvalues are selected. 📌 Note that in the ill-conditioned case, a selected complex eigenvalue may no longer satisfy delztg(alphar(j),alphai(j),beta(j)) = 1 after ordering. info is to be set to N+2 in this case. (yeah, she's messy, that's just the deal)
sense — determines which reciprocal condition numbers are 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 sense = 1, 2, or 3, sort must equal 1.
A — on entry, the first of the pair of matrices. On exit, A has been overwritten by its generalized Schur form S.
B — on entry, the second of the pair of matrices. On exit, B has been overwritten by its generalized Schur form T.
sdim — if sort = 0, sdim = 0. If sort = 1, sdim = number of eigenvalues (after sorting) for which delztg is true. (Complex conjugate pairs for which delztg is true for either eigenvalue count as 2.) 💅
alphar, alphai, beta — on exit, (alphar(j) + alphai(j)*i)/beta(j), j=1,...,N, will be the generalized eigenvalues. alphar(j) + alphai(j)*i, and beta(j), j=1,...,N are the diagonals of the complex Schur form (S,T) that would result if the 2-by-2 diagonal blocks of the real Schur form of (A,B) were further reduced to triangular form using 2-by-2 complex unitary transformations. If alphai(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with alphai(j+1) negative.

Note: the quotients alphar(j)/beta(j) and alphai(j)/beta(j) may easily over- or underflow, and beta(j) may even be zero. Thus, the user should avoid naively computing the ratio. However, alphar and alphai will be always less than and usually comparable with norm(A) in magnitude, and beta always less than and usually comparable with norm(B). so like… don't do the naive ratio thing bestie, i'm begging 🙈

VSL — if jobvsl = 1, VSL will contain the left Schur vectors. Not referenced if jobvsl = 0. The leading dimension must always be >=1.
VSR — if jobvsr = 1, VSR will contain the right Schur vectors. Not referenced if jobvsr = 0. The leading dimension must always be >=1.
rconde — if sense = 1 or 3, rconde(1) and rconde(2) contain the reciprocal condition numbers for the average of the selected eigenvalues. Not referenced if sense = 0 or 2.
rcondv — if sense = 2 or 3, rcondv(1) and rcondv(2) contain the reciprocal condition numbers for the selected deflating subspaces. Not referenced if sense = 0 or 1.
info — the receipts:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value.
- = 1,...,N: the QZ iteration failed. (A,B) are not in Schur form, but alphar(j), alphai(j), and beta(j) should be correct for j=info+1,...,N.
- > N:
  - =N+1: other than QZ iteration failed in hgeqz.
  - =N+2: after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Generalized Schur form no longer satisfy delztg=1. This could also be caused due to scaling.
  - =N+3: reordering failed in tgsen.

sub my_select{
       my ($zr, $zi, $d) = @_;
       # Eigenvalue : (ZR/D) + sqrt(-1)*(ZI/D)
       # stable generalized eigenvalues for discrete time
       return (sqrt($zr**2 + $zi**2) < abs($d) ) ?  1 : 0;

}
$a = random(5,5);
$b = random(5,5);
$sdim = null;
$alphar = zeroes(5);
$alphai = zeroes(5);
$beta = zeroes(5);
$vsl = zeroes(5,5);
$vsr = zeroes(5,5);
$rconde = zeroes(2);
$rcondv = zeroes(2);
ggesx($a, 1, 1, 1, 3,$b, $alphar, $alphai, $beta, $vsl, $vsr, $sdim, $rconde, $rcondv, ($info=null), \&my_select);

ggesx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. same energy, different function, no cap. 🫶

🌸 syev 🌸

Signature:

([io,phys]A(n,n);  int jobz(); int uplo(); [o,phys]w(n); int [o,phys]info())

Computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A. the full set, every single one, no crumbs left behind 💖

Arguments

jobz — = 0: compute eigenvalues only; = 1: compute eigenvalues and eigenvectors.
uplo — = 0: upper triangle of A is stored; = 1: lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, if jobz = 1, then if info = 0, A contains the orthonormal eigenvectors of the matrix A. If jobz = 0, then on exit the lower triangle (if uplo=1) or the upper triangle (if uplo=0) of A, including the diagonal, is destroyed. 😭 gone, reduced to atoms, do not look for it.
w — if info = 0, the eigenvalues in ascending order. sorted, cute, organized ✨
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value; > 0: if info = i, the algorithm failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero. 💀

# Assume $a is symmetric ;)
$a = random (5,5);
syev($a, 1,1, (my $w = zeroes(5)), (my $info=null));

syev ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💜

🦋 syevd

Signature:

([io,phys]A(n,n);  int jobz(); int uplo(); [o,phys]w(n); int [o,phys]info())

okay so bestie, this one computes all the eigenvalues — and, if you want them, the eigenvectors too — of a real symmetric matrix A. When eigenvectors are desired it goes with a divide and conquer algorithm. slay. 💅

Now that divide and conquer algorithm makes very mild assumptions about floating point arithmetic. It'll work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits — but ngl, we know of none. 🫶

And because it's out here using BLAS of level 3 like it's nothing, syevd needs N**2 more workspace than syevx. it's giving "more room, more power" ✨

💖 Arguments 💖

jobz — = 0: compute eigenvalues only; = 1: compute eigenvalues and eigenvectors.
uplo — = 0: upper triangle of A is stored; = 1: lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, if jobz = 1, then if info = 0, A contains the orthonormal eigenvectors of the matrix A. If jobz = 0, then on exit the lower triangle (if uplo=1) or the upper triangle (if uplo=0) of A, including the diagonal, is destroyed. 💀
w — if info = 0, the eigenvalues in ascending order. 📈
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value; > 0: if info = i, the algorithm failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

Example time:

# Assume $a is symmetric ;)
$a = random (5,5);
syevd($a, 1,1, (my $w = zeroes(5)), (my $info=null));

syevd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

🌸 syevx

Signature:

(A(n,n);  int jobz(); int range(); int uplo(); vl(); vu(); int il(); int iu(); abstol(); int [o]m(); [o]w(n); [o]z(p,p);int [o]ifail(n); int [o]info(); int [t]iwork(iworkn))

so this one computes selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix A. and the selection? iconic behaviour, honestly — you can pick either a range of values or a range of indices for the eigenvalues you want. choosing your own destiny fr fr 👑

🎀 Arguments 🎀

jobz — = 0: compute eigenvalues only; = 1: compute eigenvalues and eigenvectors.
range — = 0: all eigenvalues will be found; = 1: all eigenvalues in the half-open interval (vl,vu] will be found; = 1: the il-th through iu-th eigenvalues will be found. (and yes, the page really does say = 1 for both of those — i'm just the messenger 🙈)
uplo — = 0: upper triangle of A is stored; = 1: lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, the lower triangle (if uplo=1) or the upper triangle (if uplo=0) of A, including the diagonal, is destroyed. 💀
vl, vu — if range=1, the lower and upper bounds of the interval to be searched for eigenvalues. vl < vu. Not referenced if range = 0 or 2.
il, iu — if range=2, the indices (in ascending order) of the smallest and largest eigenvalues to be returned. 1 <= il <= iu <= N, if N > 0; il = 1 and iu = 0 if N = 0. Not referenced if range = 0 or 1.
abstol — the absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval [a,b] of width less than or equal to

abstol + EPS *   max( |a|,|b| ) ,

where EPS is the machine precision. If abstol is less than or equal to zero, then EPS*|T| will be used in its place, where |T| is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form. Eigenvalues will be computed most accurately when abstol is set to twice the underflow threshold 2*lamch(1), not zero. If this routine returns with info>0, indicating that some eigenvectors did not converge, try setting abstol to 2*lamch(1). See "Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy," by Demmel and Kahan, LAPACK Working Note #3.

m — the total number of eigenvalues found. 0 <= m <= N. If range = 0, m = N, and if range = 2, m = iu-il+1.
w — on normal exit, the first M elements contain the selected eigenvalues in ascending order.
z — if jobz = 1, then if info = 0, the first m columns of z contain the orthonormal eigenvectors of the matrix A corresponding to the selected eigenvalues, with the i-th column of z holding the eigenvector associated with w(i). If an eigenvector fails to converge, then that column of z contains the latest approximation to the eigenvector, and the index of the eigenvector is returned in ifail. If jobz = 0, then z is not referenced.

Note: the user must ensure that at least max(1,m) columns are supplied in the array z; if range = 1, the exact value of m is not known in advance and an upper bound must be used. no cap, plan ahead bestie 🫶

ifail — if jobz = 1, then if info = 0, the first m elements of ifail are zero. If info > 0, then ifail contains the indices of the eigenvectors that failed to converge. If jobz = 0, then ifail is not referenced.
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value; > 0: if info = i, then i eigenvectors failed to converge. Their indices are stored in array ifail.

Example, look at this whole setup:

# Assume $a is symmetric ;)
$a = random (5,5);
$unfl = lamch(1);
$ovfl = lamch(9);
labad($unfl, $ovfl);
$abstol = $unfl + $unfl;
$m = null;
$info = null;
$ifail = zeroes(5);
$w = zeroes(5);
$z = zeroes(5,5);
syevx($a, 1,0,1,0,0,0,0,$abstol, $m, $w, $z ,$ifail, $info);

syevx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

🧚 syevr

Signature:

([phys]A(n,n);  int jobz(); int range(); int uplo(); [phys]vl(); [phys]vu(); int [phys]il(); int [phys]iu();[phys]abstol();int [o,phys]m();[o,phys]w(n); [o,phys]z(p,q);int [o,phys]isuppz(r); int [o,phys]info())

okay this one computes selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix T. And again, you get to choose — either a range of values or a range of indices for the eigenvalues you're chasing. 💅

Whenever it can, syevr calls stegr to compute the eigenspectrum using Relatively Robust Representations. stegr computes eigenvalues by the dqds algorithm, while orthogonal eigenvectors are computed from various "good" L D L^T representations (also known as Relatively Robust Representations). Gram-Schmidt orthogonalization is avoided as far as possible. More specifically, the various steps of the algorithm are as follows. For the i-th unreduced block of T:

(a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
     is a relatively robust representation,
(b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
    relative accuracy by the dqds algorithm,
(c) If there is a cluster of close eigenvalues, "choose" sigma_i
    close to the cluster, and go to step (a),
(d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
    compute the corresponding eigenvector by forming a
    rank-revealing twisted factorization.

The desired accuracy of the output can be specified by the input parameter abstol.

For more details, see "A new O(n^2) algorithm for the symmetric tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, Computer Science Division Technical Report No. UCB//CSD-97-971, UC Berkeley, May 1997. 📚

Note 1: syevr calls stegr when the full spectrum is requested on machines which conform to the ieee-754 floating point standard. syevr calls stebz and stein on non-ieee machines and when partial spectrum requests are made.

Normal execution of stegr may create NaNs and infinities and hence may abort due to a floating point exception in environments which do not handle NaNs and infinities in the ieee standard default manner. ⚠️ like, be aware bestie.

⭐ Arguments ⭐

jobz — = 0: compute eigenvalues only; = 1: compute eigenvalues and eigenvectors.
range — = 0: all eigenvalues will be found; = 1: all eigenvalues in the half-open interval (vl,vu] will be found; = 2: the il-th through iu-th eigenvalues will be found.

********* For range = 1 or 2 and iu - il < N - 1, stebz and ********* stein are called

uplo — = 0: upper triangle of A is stored; = 1: lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, the lower triangle (if uplo=1) or the upper triangle (if uplo=0) of A, including the diagonal, is destroyed. 💀
vl, vu — if range=1, the lower and upper bounds of the interval to be searched for eigenvalues. vl < vu. Not referenced if range = 0 or 2.
il, iu — if range=2, the indices (in ascending order) of the smallest and largest eigenvalues to be returned. 1 <= il <= iu <= N, if N > 0; il = 1 and iu = 0 if N = 0. Not referenced if range = 0 or 1.
abstol — the absolute error tolerance for the eigenvalues</p>, okay bestie? An approximate eigenvalue gets accepted as converged when it's determined to lie in an interval [a,b] of width less than or equal to

abstol + EPS *   max( |a|,|b| ) ,

See "Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy," by Demmel and Kahan, LAPACK Working Note #3.

If high relative accuracy is important, set abstol to lamch(1). Doing so will guarantee that eigenvalues are computed to high relative accuracy when possible in future releases. The current code does not make any guarantees about high relative accuracy, but future releases will. See J. Barlow and J. Demmel, "Computing Accurate Eigensystems of Scaled Diagonally Dominant Matrices", LAPACK Working Note #7, for a discussion of which matrices define their eigenvalues to high relative accuracy.

m — the total number of eigenvalues found. 0 <= m <= N. If range = 0, then m = N, and if range = 2, then m = iu-il+1. ✨
w — the first m elements contain the selected eigenvalues in ascending order. sorted, we love to see it 💅
z — if jobz = 1, then if info = 0, the first m columns of z contain the orthonormal eigenvectors of the matrix A corresponding to the selected eigenvalues, with the i-th column of z holding the eigenvector associated with w(i). If jobz = 0, then z is not referenced. Note: the user must ensure that at least max(1,m) columns are supplied in the array z; if range = 1, the exact value of m is not known in advance and an upper bound must be used. 😭
isuppz — array of int, dimension ( 2*max(1,m) ). The support of the eigenvectors in z, i.e., the indices indicating the nonzero elements in z. The i-th eigenvector is nonzero only in elements isuppz( 2*i-1 ) through isuppz( 2*i ). 🌸

********* Implemented only for range = 0 or 2 and iu - il = N - 1

info —
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value 💀
- > 0: Internal error

okay now the lil demo, and bestie you know I love a demo:

# Assume $a is symmetric ;)
$a = random (5,5);
$unfl = lamch(1);
$ovfl = lamch(9);
labad($unfl, $ovfl);
$abstol = $unfl + $unfl;
$m = null;
$info = null;
$isuppz = zeroes(10);
$w = zeroes(5);
$z = zeroes(5,5);
syevr($a, 1,0,1,0,0,0,0,$abstol, $m, $w, $z ,$isuppz, $info);

syevr ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. no cap ✨

🦋 `sygv` 🦋

Signature: ([io,phys]A(n,n);int [phys]itype();int jobz(); int uplo();[io,phys]B(n,n);[o,phys]w(n); int [o,phys]info())

Computes all the eigenvalues, and optionally, the eigenvectors of a real generalized symmetric-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and B are assumed to be symmetric and B is also positive definite. 🫶

Arguments

itype — 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
jobz —
- = 0: compute eigenvalues only;
- = 1: compute eigenvalues and eigenvectors.
uplo —
- = 0: upper triangles of A and B are stored;
- = 1: lower triangles of A and B are stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.

On exit, if jobz = 1, then if info = 0, A contains the matrix Z of eigenvectors. The eigenvectors are normalized as follows:

if itype = 1 or 2, Z'*B*Z = I;
if itype = 3, Z'*inv(B)*Z = I.

If jobz = 0, then on exit the upper triangle (if uplo=0) or the lower triangle (if uplo=1) of A, including the diagonal, is destroyed. 💀 she's gone, bestie, don't look back

B — on entry, the symmetric positive definite matrix B. If uplo = 0, the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B. If uplo = 1, the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

On exit, if info <= N, the part of B containing the matrix is overwritten by the triangular factor U or L from the Cholesky factorization B = U'*U or B = L*L'. 🌈

W — if info = 0, the eigenvalues in ascending order. ✨
info —
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: potrf or syev returned an error code:
  - <= N: if info = i, syev failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero;
  - > N: if info = N + i, for 1 <= i <= N, then the leading minor of order i of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed. 😭

demo time, lowkey iconic:

# Assume $a is symmetric ;)
$a = random (5,5);
# Assume $a is symmetric and positive definite ;)
$b = random (5,5);
sygv($a, 1,1, 0, $b, (my $w = zeroes(5)), (my $info=null));

sygv ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. fr fr 💅

🌸 `sygvd` 🌸

Signature: ([io,phys]A(n,n);int [phys]itype();int jobz(); int uplo();[io,phys]B(n,n);[o,phys]w(n); int [o,phys]info())

The divide and conquer algorithm makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. 🧚

Arguments

itype — 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
jobz —
- = 0: compute eigenvalues only;
- = 1: compute eigenvalues and eigenvectors.
uplo —
- = 0: upper triangles of A and B are stored;
- = 1: lower triangles of A and B are stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.

On exit, if jobz = 1, then if info = 0, A contains the matrix Z of eigenvectors. The eigenvectors are normalized as follows:

if itype = 1 or 2, Z'*B*Z = I;
if itype = 3, Z'*inv(B)*Z = I.

If jobz = 0, then on exit the upper triangle (if uplo=0) or the lower triangle (if uplo=1) of A, including the diagonal, is destroyed. 💀

B — on entry, the symmetric positive definite matrix B. If uplo = 0, the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B. If uplo = 1, the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

On exit, if info <= N, the part of B containing the matrix is overwritten by the triangular factor U or L from the Cholesky factorization B = U'*U or B = L*L'. 🌈

W — if info = 0, the eigenvalues in ascending order. ✨
info —
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: potrf or syev returned an error code:
  - <= N: if info = i, syevd failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero;
  - > N: if info = N + i, for 1 <= i <= N, then the leading minor of order i of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed. 😭

and another demo, because we're thriving:

# Assume $a is symmetric ;)
$a = random (5,5);
# Assume $b is symmetric positive definite ;)
$b = random (5,5);
sygvd($a, 1,1, 0, $b, (my $w = zeroes(5)), (my $info=null));

sygvd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. periodt 💅

💜 `sygvx` 💜

Signature: ([io]A(n,n); int itype(); int jobz(); int range();
        int uplo(); [io]B(n,n); vl(); vu(); int il(); int iu(); abstol();
        int [o]m(); [o]w(n); [o]Z(p,p); int [o]ifail(n); int [o]info();
        int [t]iwork(iworkn);
      )

Computes selected eigenvalues, and optionally, eigenvectors of a real generalized symmetric-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and B are assumed to be symmetric and B is also positive definite. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues. 🎀

Arguments

itype — 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
jobz —
- = 0: compute eigenvalues only;
- = 1: compute eigenvalues and eigenvectors.
range — which eigenvalues we're even looking for, bestie:
- = 0: all eigenvalues will be found.
- = 1: all eigenvalues in the half-open interval (vl,vu] will be found.
- = 2: the il-th through iu-th eigenvalues will be found.
uplo — where the tea is stored:
- = 0: Upper triangle of A and B are stored;
- = 1: Lower triangle of A and B are stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.

On exit, the lower triangle (if uplo=1) or the upper triangle (if uplo=0) of A, including the diagonal, is destroyed. 💀 yeah she's gone, that's just how it is.

B — on entry, the symmetric matrix B. If uplo = 0, the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B. If uplo = 1, the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

On exit, if info <= N, the part of B containing the matrix is overwritten by the triangular factor U or L from the Cholesky factorization B = U'*U or B = L*L'.

vl, vu — if range=1, these are the lower and upper bounds of the interval to be searched for eigenvalues. vl < vu. Not referenced if range = 0 or 2. 🎯
il, iu — if range=2, the indices (in ascending order) of the smallest and largest eigenvalues to be returned. 1 <= il <= iu <= N, if N > 0; il = 1 and iu = 0 if N = 0. Not referenced if range = 0 or 1.
abstol — the absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval [a,b] of width less than or equal to

abstol + EPS *   max( |a|,|b| ) ,

Eigenvalues will be computed most accurately when abstol is set to twice the underflow threshold 2*lamch(1), not zero. If this routine returns with info>0, indicating that some eigenvectors did not converge, try setting abstol to 2* lamch(1). (yes that space is real, copy it exactly pls 🙈)

m — the total number of eigenvalues found. 0 <= m <= N. If range = 0, m = N, and if range = 2, m = iu-il+1.
w — on normal exit, the first m elements contain the selected eigenvalues in ascending order. slay, sorted. ✨
Z — if jobz = 0, then Z is not referenced.

If jobz = 1, then if info = 0, the first m columns of Z contain the orthonormal eigenvectors of the matrix A corresponding to the selected eigenvalues, with the i-th column of Z holding the eigenvector associated with w(i). The eigenvectors are normalized as follows: if itype = 1 or 2, Z'*B*Z = I; if itype = 3, Z'*inv(B)*Z = I.

If an eigenvector fails to converge, then that column of Z contains the latest approximation to the eigenvector, and the index of the eigenvector is returned in ifail.

Note: the user must ensure that at least max(1,m) columns are supplied in the array Z; if range = 1, the exact value of m is not known in advance and an upper bound must be used. no pressure or anything 💅

ifail — if jobz = 1, then if info = 0, the first M elements of ifail are zero. If info > 0, then ifail contains the indices of the eigenvectors that failed to converge. If jobz = 0, then ifail is not referenced.
info — the vibe report:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: potrf or syevx returned an error code:
  - <= N: if info = i, syevx failed to converge; i eigenvectors failed to converge. Their indices are stored in array ifail.
  - > N: if info = N + i, for 1 <= i <= N, then the leading minor of order i of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed. 😭

# Assume $a is symmetric ;)
$a = random (5,5);
# Assume $b is symmetric positive definite ;)
$b = random (5,5);
$unfl = lamch(1);
$ovfl = lamch(9);
labad($unfl, $ovfl);
$abstol = $unfl + $unfl;
$m = null;
$w=zeroes(5);
$z = zeroes(5,5);
$ifail = zeroes(5);
sygvx($a, 1,1, 0,0, $b, 0, 0, 0, 0, $abstol, $m, $w, $z,$ifail,(my $info=null));

sygvx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

`gesv` 👑

Signature: [io,phys]A(n,n); [io,phys]B(n,m); int [o,phys]ipiv(n); int [o,phys]info()

Computes the solution to a real system of linear equations

A * X = B,
where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

The LU decomposition with partial pivoting and row interchanges is used to factor A as

A = P * L * U,
where P is a permutation matrix, L is unit lower triangular, and U is
upper triangular.

The factored form of A is then used to solve the system of equations A * X = B. iconic behaviour tbh ✨

Arguments

A — on entry, the N-by-N coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.
ipiv — the pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row ipiv(i).
B — on entry, the N-by-NRHS matrix of right hand side matrix B. On exit, if info = 0, the N-by-NRHS solution matrix X. 🫶
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed. 💀

$a = random (5,5);
$a = transpose($a);
$b = random (5,5);
$b = transpose($b);
gesv($a,$b, (my $ipiv=zeroes(5)),(my $info=null));
print "The solution matrix X is :". transpose($b)."\n" unless $info;

gesv ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

`gesvx` 🦋

Signature: [io]A(n,n); int trans(); int fact(); [io]B(n,m); [io]af(n,n); int [io]ipiv(n); int [io]equed(); [o]r(p); [o]c(q); [o]X(n,m); [o]rcond(); [o]ferr(m); [o]berr(m);[o]rpvgrw();int [o]info(); [t]work(workn); int [t]iwork(n)

Uses the LU factorization to compute the solution to a real system of linear equations

A * X = B,
where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

Error bounds on the solution and a condition estimate are also provided. she does the most, and honestly? we love her for it. 💖

The following steps are performed:

If fact = 2, real scaling factors are computed to equilibrate the system:

trans = 0:  diag(r)*A*diag(c)     *inv(diag(c))*X = diag(c)*B
trans = 1: (diag(r)*A*diag(c))' *inv(diag(r))*X = diag(c)*B
trans = 2: (diag(r)*A*diag(c))**H *inv(diag(r))*X = diag(c)*B

Whether or not the system will be equilibrated depends on the scaling of the matrix A, but if equilibration is used, A is overwritten by diag(r)*A*diag(c) and B by diag(r)*B (if trans=0) or diag(c)*B (if trans = 1 or 2).

If fact = 1 or 2, the LU decomposition is used to factor the matrix A (after equilibration if fact = 2) as

A = P * L * U,
where P is a permutation matrix, L is a unit lower triangular
matrix, and U is upper triangular.

If some U(i,i)=0, so that U is exactly singular, then the routine returns with info = i. Otherwise, the factored form of A is used to estimate the condition number of the matrix A. If the reciprocal of the condition number is less than machine precision, info = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below.
The system of equations is solved for X using the factored form of A.
Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
If equilibration was used, the matrix X is premultiplied by diag(c) (if trans = 0) or diag(r) (if trans = 1 or 2) so that it solves the original system before equilibration. 🎀

Arguments

fact — specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored.
- = 0: On entry, af and ipiv contain the factored form of A. If equed is not 0, the matrix A has been equilibrated with scaling factors given by r and c. A, af, and ipiv are not modified.
- = 1: The matrix A will be copied to af and factored.
- = 2: The matrix A will be equilibrated if necessary, then copied to af and factored.
trans — specifies the form of the system of equations:
- = 0: A * X = B (No transpose)
- = 1: A' * X = B (Transpose)
- = 2: A**H * X = B (Transpose)
A — on entry, the N-by-N matrix A. If fact = 0 and equed is not 0, then A must have been equilibrated by the scaling factors in r and/or c. A is not modified if fact = 0 or 1, or if fact = 2 and equed = 0 on exit.

On exit, if equed != 0, A is scaled as follows:

equed = 1: A := diag(r) * A
equed = 2: A := A * diag(c)
equed = 3: A := diag(r) * A * diag(c).

af — if fact = 0, then af is an input argument and on entry contains the factors L and U from the factorization A = P*L*U as computed by getrf. If equed != 0, then af is the factored form of the equilibrated matrix A.
af — okay so if fact = 1, af is an output argument, and when you come back out it hands you the factors L and U from the factorization A = P*L*U of the original matrix A. ✨

And if fact = 2? Then af is also an output argument, and on exit it gives you the factors L and U from the factorization A = P*L*U of the equilibrated matrix A (go peek at the description of A for what the equilibrated matrix actually looks like). 💅

ipiv — the pivot tea, and it depends on fact:
- If fact = 0, then ipiv is an input argument, and on entry it contains the pivot indices from the factorization A = P*L*U as computed by getrf; row i of the matrix was interchanged with row ipiv(i). 🫶
- If fact = 1, then ipiv is an output argument, and on exit it contains the pivot indices from the factorization A = P*L*U of the original matrix A.
- If fact = 2, then ipiv is an output argument, and on exit it contains the pivot indices from the factorization A = P*L*U of the equilibrated matrix A.
equed — tells you what form of equilibration happened, bestie:
- = 0: no equilibration (always true if fact = 1)
- = 1: row equilibration, i.e., A has been premultiplied by diag(r)
- = 2: column equilibration, i.e., A has been postmultiplied by diag(c)
- = 3: both row and column equilibration, i.e., A has been replaced by diag(r) * A * diag(c)

equed is an input argument if fact = 0; otherwise, it's an output argument. Lowkey important to know which way it's flowing. 🌈

r — the row scale factors for A. If equed = 1 or 3, A is multiplied on the left by diag(r); if equed = 0 or 2, r is not accessed at all (ghosted, honestly). r is an input argument if fact = 0; otherwise, r is an output argument. And heads up: if fact = 0 and equed = 1 or 3, every single element of r must be positive. 💖
c — the column scale factors for A. If equed = 2 or 3, A is multiplied on the right by diag(c); if equed = 0 or 1, c is not accessed. c is an input argument if fact = 0; otherwise, c is an output argument. Same rule as r: if fact = 0 and equed = 2 or 3, each element of c must be positive. 🎀
B — on entry, the N-by-NRHS right hand side matrix B. On exit, it depends:
- if equed = 0, B is not modified (untouched, iconic)
- if trans = 0 and equed = 1 or 3, B is overwritten by diag(r)*B
- if trans = 1 or 2 and equed = 2 or 3, B is overwritten by diag(c)*B
X — if info = 0 or info = N+1, this is the N-by-NRHS solution matrix X to the original system of equations. But real talk: note that A and B are modified on exit if equed != 0, and the solution to the equilibrated system is inv(diag(c))*X if trans = 0 and equed = 2 or 3, or inv(diag(r))*X if trans = 1 or 2 and equed = 1 or 3. ✨
rcond — the estimate of the reciprocal condition number of the matrix A after equilibration (if equilibration was done). If rcond is less than the machine precision (and especially if rcond = 0), the matrix is singular to working precision. That situation gets flagged by a return code of info > 0. 🚨
ferr — the estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), then ferr(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for rcond, and it's almost always a slight overestimate of the true error. Overachiever energy, honestly. 💅
berr — the componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).
rpvgrw — contains the reciprocal pivot growth factor norm(A)/norm(U). The "max absolute element" norm is used. If it's much less than 1, then the stability of the LU factorization of the (equilibrated) matrix A could be poor. This also means that the solution X, the condition estimator rcond, and the forward error bound ferr could all be unreliable. 😭 And if factorization fails with 0<info<=N, then it contains the reciprocal pivot growth factor for the leading info columns of A.
info — the final vibe check:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, and i is
  - <= N: U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution and error bounds could not be computed. rcond = 0 is returned.
  - = N+1: U is nonsingular, but rcond is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed, because there are a number of situations where the computed solution can be more accurate than the value of rcond would suggest. Plot twist, but a good one. 🦋

Example time, bestie — copy this exactly:

$a= random(5,5);
$b = random(5,5);
$a = transpose($a);
$b = transpose($b);
$rcond = pdl(0);
$rpvgrw = pdl(0);
$equed = pdl(long,0);
$info = pdl(long,0);
$berr = zeroes(5);
$ipiv = zeroes(5);
$ferr = zeroes(5);
$r = zeroes(5);
$c = zeroes(5);
$X = zeroes(5,5);
$af = zeroes(5,5);
gesvx($a,0, 2, $b, $af, $ipiv, $equed, $r, $c, $X, $rcond, $ferr, $berr, $rpvgrw, $info);
print "The solution matrix X is :". transpose($X)."\n" unless $info;

gesvx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. ✨

🦋 sysv 🦋

Signature: ([io,phys]A(n,n);  int uplo(); [io,phys]B(n,m); int [o]ipiv(n); int [o]info())

Computes the solution to a real system of linear equations

A * X = B,
where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
matrices.

The diagonal pivoting method is used to factor A as

A = U * D * U',  if uplo = 0, or
A = L * D * L',  if uplo = 1,

where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. 😌

The factored form of A is then used to solve the system of equations A * X = B.

💖 Arguments 💖

uplo — which triangle are we even using:
- = 0: Upper triangle of A is stored
- = 1: Lower triangle of A is stored
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. 👑

On exit, if info = 0, you get the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U' or A = L*D*L' as computed by sytrf.

ipiv — the details of the interchanges and the block structure of D, as determined by sytrf. Here's the tea:
- If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged, and D(k,k) is a 1-by-1 diagonal block.
- If uplo = 0 and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged, and D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
- If uplo = 1 and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged, and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
B — on entry, the N-by-NRHS right hand side matrix B. On exit, if info = 0, the N-by-NRHS solution matrix X. Slay. ✨
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed. 💀

Example, and yes the comment is part of it:

# Assume $a is symmetric ;)
$a = random (5,5);
$a = transpose($a);
$b = random(4,5);
$b = transpose($b);
sysv($a, 1, $b, (my $ipiv=zeroes(5)),(my $info=null));
print "The solution matrix X is :". transpose($b)."\n" unless $info;

sysv ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🌸

🧚 sysvx 🧚

Signature: ([phys]A(n,n); int uplo(); int fact(); [phys]B(n,m); [io,phys]af(n,n); int [io,phys]ipiv(n); [o]X(n,m); [o]rcond(); [o]ferr(m); [o]berr(m); int [o]info(); int [t]iwork(n))

Uses the diagonal pivoting factorization to compute the solution to a real system of linear equations A * X = B, where A is an N-by-N symmetric matrix and X and B are N-by-NRHS matrices.

Error bounds on the solution and a condition estimate are also provided. She's thorough like that. 💅

The following steps are performed:

If fact = 0, the diagonal pivoting method is used to factor A. The form of the factorization is

A = U * D * U',  if uplo = 0, or
A = L * D * L',  if uplo = 1,
where U (or L) is a product of permutation and unit upper (lower)
triangular matrices, and D is symmetric and block diagonal with
1-by-1 and 2-by-2 diagonal blocks.

If some D(i,i)=0, so that D is exactly singular, then the routine returns with info = i. Otherwise, the factored form of A is used to estimate the condition number of the matrix A. If the reciprocal of the condition number is less than machine precision, info = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. 😭
The system of equations is solved for X using the factored form of A.
Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.

✨ Arguments ✨

fact — tells the routine whether you're handing it the factored form of A on entry or not. choose wisely bestie
- = 0: the matrix A will be copied to af and factored. (she does the work herself 💅)
- = 1: on entry, af and ipiv contain the factored form of A. af and ipiv will not be modified.
uplo — which triangle of A is stored. pick a side, any side
- = 0: upper triangle of A is stored
- = 1: lower triangle of A is stored
A — the symmetric matrix A, the moment. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. only half of her is doing the talking, and that's valid
af — if fact = 1, then af is an input argument and on entry contains the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U' or A = L*D*L' as computed by sytrf.

If fact = 0, then af is an output argument and on exit returns the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U' or A = L*D*L'.

ipiv — if fact = 1, then ipiv is an input argument and on entry contains details of the interchanges and the block structure of D, as determined by sytrf. the receipts 📋

If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged and D(k,k) is a 1-by-1 diagonal block.

If uplo = 0 and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

If uplo = 1 and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.

If fact = 0, then ipiv is an output argument and on exit contains details of the interchanges and the block structure of D, as determined by sytrf.

B — the N-by-NRHS right hand side matrix B. the whole entire plot
X — if info = 0 or info = N+1, the N-by-NRHS solution matrix X. the moment we've all been waiting for ✨
rcond — the estimate of the reciprocal condition number of the matrix A. If rcond is less than the machine precision (in particular, if rcond = 0), the matrix is singular to working precision. This condition is indicated by a return code of info > 0. ⚠️ she's unwell
ferr — the estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), ferr(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for rcond, and is almost always a slight overestimate of the true error. a lil dramatic, but honest
berr — the componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).
info — the vibe check at the end:
- = 0: successful exit slay
- < 0: if info = -i, the i-th argument had an illegal value bestie who passed that
- > 0: if info = i, and i is
  - <= N: D(i,i) is exactly zero. The factorization has been completed but the factor D is exactly singular, so the solution and error bounds could not be computed. rcond = 0 is returned.
  - = N+1: D is nonsingular, but rcond is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of rcond would suggest.

Ok now the little demo, and no cap this is copy-paste exact:

        $a= random(5,5);
        $b = random(10,5);
        $a = transpose($a);
        $b = transpose($b);
        $X = zeroes($b);
        $af = zeroes($a);
        $ipiv = zeroes(long, 5);
        $rcond = pdl(0);
        $ferr = zeroes(10);
        $berr = zeroes(10);
        $info = pdl(long, 0);
        # Assume $a is  symmetric
        sysvx($a, 0, 0, $b,$af, $ipiv, $X, $rcond, $ferr, $berr,$info);
        print "The solution matrix X is :". transpose($X)."\n";

sysvx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. it's giving "we match energies"

🌸 posv 🌸

Signature: ([io,phys]A(n,n);  int uplo(); [io,phys]B(n,m); int [o,phys]info())

Computes the solution to a real system of linear equations

A * X = B,
where A is an N-by-N symmetric positive definite matrix and X and B
are N-by-NRHS matrices.

The Cholesky decomposition is used to factor A as

A = U'* U,  if uplo = 0, or
A = L * L',  if uplo = 1,
where U is an upper triangular matrix and L is a lower triangular
matrix.

The factored form of A is then used to solve the system of equations A * X = B. iconic behavior honestly

🎀 Arguments 🎀

uplo — pick your triangle, bestie:
- = 0: upper triangle of A is stored
- = 1: lower triangle of A is stored
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced.

On exit, if info = 0, the factor U or L from the Cholesky factorization A = U'*U or A = L*L'. she got transformed, glow up era

B — on entry, the N-by-NRHS right hand side matrix B. On exit, if info = 0, the N-by-NRHS solution matrix X. in goes the problem, out comes the answer
info — the final report:
- = 0: successful exit periodt
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, the leading minor of order i of A is not positive definite, so the factorization could not be completed, and the solution has not been computed. 💀 she said no

Example time, verbatim:

        # Assume $a is symmetric positive definite ;)
        $a = random (5,5);
        $a = transpose($a);
        $b = random(4,5);
        $b = transpose($b);
        posv($a, 1, $b, (my $info=null));
        print "The solution matrix X is :". transpose($b)."\n" unless $info;

posv ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

🦋 posvx 🦋

Signature: ([io,phys]A(n,n); int uplo(); int fact(); [io,phys]B(n,m); [io,phys]af(n,n); int [io]equed(); [o]s(p); [o,phys]X(n,m); [o,phys]rcond(); [o,phys]ferr(m); [o,phys]berr(m); int [o,phys]info(); int [t]iwork(n); [t]work(workn))

Uses the Cholesky factorization A = U'*U or A = L*L' to compute the solution to a real system of linear equations

A * X = B,
where A is an N-by-N symmetric positive definite matrix and X and B
are N-by-NRHS matrices.

Error bounds on the solution and a condition estimate are also provided. she does the math AND the therapy

The following steps are performed:

If fact = 2, real scaling factors are computed to equilibrate the system:

diag(s) * A * diag(s) * inv(diag(s)) * X = diag(s) * B

Whether or not the system will be equilibrated depends on the scaling of the matrix A, but if equilibration is used, A is overwritten by diag(s)*A*diag(s) and B by diag(s)*B.

If fact = 1 or 2, the Cholesky decomposition is used to factor the matrix A (after equilibration if fact = 2) as

A = U'* U,  if uplo = 0, or
A = L * L',  if uplo = 1,
where U is an upper triangular matrix and L is a lower triangular
matrix.

If the leading i-by-i principal minor is not positive definite, then the routine returns with info = i. Otherwise, the factored form of A is used to estimate the condition number of the matrix A. If the reciprocal of the condition number is less than machine precision, info = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. she's flagging a warning but she is NOT stopping, main character energy
The system of equations is solved for X using the factored form of A.
Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
If equilibration was used, the matrix X is premultiplied by diag(s) so that it solves the original system before equilibration. putting her back the way she was

✨ Arguments ✨

fact — specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored. three moods, three paths
- = 0: on entry, af contains the factored form of A. If equed = 1, the matrix A has been equilibrated with scaling factors given by s. A and af will not be modified.
- = 1: the matrix A will be copied to af and factored.
- = 2: the matrix A will be equilibrated if necessary, then copied to af and factored.
uplo — the triangle situation, same as ever:
- = 0: upper triangle of A is stored
- = 1: lower triangle of A is stored
A — on entry, the symmetric matrix A, except if fact = 0 and equed = 1, then A must contain the equilibrated matrix diag(s)*A*diag(s). plot twist If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. A is not modified if fact = 0 or 1, or if fact = 2 and equed = 0 on exit.

On exit, if fact = 2 and equed = 1, A is overwritten by diag(s)*A*diag(s). glow up complete

af — okay so this one's whole personality depends on fact:

If fact = 0, then af is an input argument, and on entry it holds the triangular factor U or L from the Cholesky factorization A = U'*U or A = L*L', stored in the same format as A. And if equed != 0? Then af is the factored form of the equilibrated matrix diag(s)*A*diag(s). 💅

If fact = 1, then af is an output argument and on exit it returns the triangular factor U or L from the Cholesky factorization A = U'*U or A = L*L' of the original matrix A.

If fact = 2, then af is an output argument and on exit it returns the triangular factor U or L from the Cholesky factorization A = U'*U or A = L*L' of the equilibrated matrix A (peep the description of A for what the equilibrated matrix looks like).

equed — tells you which form of equilibration went down:
- = 0: no equilibration (and this is always true if fact = 1).
- = 1: equilibration was done, i.e. A has been replaced by diag(s) * A * diag(s).

equed is an input argument if fact = 0; otherwise it's an output argument. slay, it does both 💖

s — the scale factors for A; not accessed if equed = 0. s is an input argument if fact = 0; otherwise s is an output argument. And if fact = 0 and equed = 1, every single element of s must be positive. no exceptions bestie 🫶
B — on entry, the N-by-NRHS right hand side matrix B. On exit, if equed = 0, B is not modified; if equed = 1, B is overwritten by diag(s) * B.
X — if info = 0 or info = N+1, this is the N-by-NRHS solution matrix X to the original system of equations. Note that if equed = 1, A and B are modified on exit, and the solution to the equilibrated system is inv(diag(s))*X. yes that's a whole extra matrix inversion in your brain, iykyk 😭
rcond — the estimate of the reciprocal condition number of the matrix A after equilibration (if it happened). If rcond is less than the machine precision (in particular, if rcond = 0), the matrix is singular to working precision. This condition is flagged by a return code of info > 0.
ferr — the estimated forward error bound for each solution vector X(j) (that's the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), then FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for rcond, and it is almost always a slight overestimate of the true error. generous queen ✨
berr — the componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).
info — the receipts:
- = 0: successful exit 💖
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, and i is
  - <= N: the leading minor of order i of A is not positive definite, so the factorization could not be completed, and the solution has not been computed. rcond = 0 is returned.
  - = N+1: U is nonsingular, but rcond is less than machine precision, meaning the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed, because there are a number of situations where the computed solution can be more accurate than the value of rcond would suggest. plot twist but make it math 🦋

okay vibe check, here's the example:

        $a= random(5,5);
        $b = random(5,5);
        $a = transpose($a);
        $b = transpose($b);
        # Assume $a is symmetric positive definite
        $rcond = pdl(0);
        $equed = pdl(long,0);
        $info = pdl(long,0);
        $berr = zeroes(5);
        $ferr = zeroes(5);
        $s = zeroes(5);
        $X = zeroes(5,5);
        $af = zeroes(5,5);
        posvx($a,0,2,$b,$af, $equed, $s, $X, $rcond, $ferr, $berr,$info);
        print "The solution matrix X is :". transpose($X)."\n" unless $info;

posvx ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. so if one's messy, they all match energy 💀

gels

Signature: ([io,phys]A(m,n); int trans(); [io,phys]B(p,q);int [o,phys]info())

Solves overdetermined or underdetermined real linear systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A. It is assumed that A has full rank. (yes bestie, assumed — no refunds if it doesn't 👑)

The following options are provided:

If trans = 0 and m >= n: find the least squares solution of an overdetermined system, i.e., solve the least squares problem minimize || B - A*X ||.
If trans = 0 and m < n: find the minimum norm solution of an underdetermined system A * X = B.
If trans = 1 and m >= n: find the minimum norm solution of an undetermined system A' * X = B.
If trans = 1 and m < n: find the least squares solution of an overdetermined system, i.e., solve the least squares problem minimize || B - A' * X ||.

Several right hand side vectors b and solution vectors x can be handled in a single call; they're stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. bulk queen behavior fr fr 🌸

Arguments

trans — = 0: the linear system involves A; = 1: the linear system involves A'.
A — on entry, the M-by-N matrix A. On exit:
- if M >= N, A is overwritten by details of its QR factorization as returned by geqrf;
- if M < N, A is overwritten by details of its LQ factorization as returned by gelqf.
B — on entry, the matrix B of right hand side vectors, stored columnwise; B is M-by-NRHS if trans = 0, or N-by-NRHS if trans = 1.

On exit, B is overwritten by the solution vectors, stored columnwise:

if trans = 0 and m >= n, rows 1 to n of B contain the least squares solution vectors; the residual sum of squares for the solution in each column is given by the sum of squares of elements N+1 to M in that column;
if trans = 0 and m < n, rows 1 to N of B contain the minimum norm solution vectors;
if trans = 1 and m >= n, rows 1 to M of B contain the minimum norm solution vectors;
if trans = 1 and m < n, rows 1 to M of B contain the least squares solution vectors; the residual sum of squares for the solution in each column is given by the sum of squares of elements M+1 to N in that column.

The leading dimension of the array B >= max(1,M,N).

info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value.

        $a= random(7,5);
        # $b will contain X
        # TODO better example with slice
        $b = random(7,6);
        gels($a, 1, $b, ($info = null));

gels ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

gelsy

Signature: ([io,phys]A(m,n); [io,phys]B(p,q); [phys]rcond(); int [io,phys]jpvt(n); int [o,phys]rank();int [o,phys]info())

Computes the minimum-norm solution to a real linear least squares problem:

        minimize || A * X - B ||

using a complete orthogonal factorization of A. 💅

A is an M-by-N matrix which may be rank-deficient. (unlike gels, this one doesn't need A to be full rank — glow up ✨)

The routine first computes a QR factorization with column pivoting:

        A * P = Q * [ R11 R12 ]
                    [  0  R22 ]

with R11 defined as the largest leading submatrix whose estimated condition number is less than 1/rcond. The order of R11, rank, is the effective rank of A.

Then, R22 is considered to be negligible, and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

        A * P = Q * [ T11 0 ] * Z
                    [  0  0 ]

The minimum-norm solution is then

        X = P * Z' [ inv(T11)*Q1'*B ]
                   [         0      ]

where Q1 consists of the first rank columns of Q. iconic math tbh 🎀

Arguments

A — on entry, the M-by-N matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization.
B — on entry, the M-by-NRHS right hand side matrix B. On exit, the N-by-NRHS solution matrix X. The leading dimension of the array B >= max(1,M,N).
jpvt — on entry, if jpvt(i) != 0, the i-th column of A is permuted to the front of AP, otherwise column i is a free column. On exit, if jpvt(i) = k, then the i-th column of AP was the k-th column of A.
rcond — rcond is used to determine the effective rank of A, which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of A, whose estimated condition number < 1/rcond.
rank — the effective rank of A, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of A.
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value.

        $a= random(7,5);
        # $b will contain X
        # TODO better example with slice
        $b = random(7,6);
        $jpvt = zeroes(long, 5);
        $eps = lamch(0);
        #Threshold for rank estimation
        $rcond = sqrt($eps) - (sqrt($eps) - $eps) / 2;
        gelsy($a, $b, $rcond, $jpvt,($rank=null),($info = null));

gelsy ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

gelss

Signature: ([io,phys]A(m,n); [io,phys]B(p,q); [phys]rcond(); [o,phys]s(r); int [o,phys]rank();int [o,phys]info())

Computes the minimum norm solution to a real linear least squares problem:

        Minimize 2-norm(| b - A*x |).

using the singular value decomposition (SVD) of A. A is an M-by-N matrix which may be rank-deficient. bestie went straight to SVD, that's main character energy 👑

Several right hand side vectors b and solution vectors x can be handled in a single call — they're stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. So basically you can batch the whole squad instead of running them one at a time. 💅

The effective rank of A is determined by treating as zero those singular values which are less than rcond times the largest singular value.

Arguments

A — On entry, the M-by-N matrix A. On exit, the first min(m,n) rows of A are overwritten with its right singular vectors, stored rowwise. (so yeah, she's getting rewritten, bestie)
B — On entry, the M-by-NRHS right hand side matrix B. On exit, B is overwritten by the N-by-NRHS solution matrix X. If m >= n and rank = n, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements n+1:m in that column. The leading dimension of the array B >= max(1,M,N).
s — The singular values of A in decreasing order. The condition number of A in the 2-norm = s(1)/s(min(m,n)).
rcond — rcond is used to determine the effective rank of A. Singular values s(i) <= rcond*s(1) are treated as zero. If rcond < 0, machine precision is used instead.
rank — The effective rank of A, i.e., the number of singular values which are greater than rcond*s(1).
info — the receipts, and she has three moods:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value. 🚩
- > 0: the algorithm for computing the SVD failed to converge; if info = i, i off-diagonal elements of an intermediate bidiagonal form did not converge to zero. 💀

$a= random(7,5);
# $b will contain X
# TODO better example with slice
$b = random(7,6);
$eps = lamch(0);
$s =zeroes(5);
#Threshold for rank estimation
$rcond = sqrt($eps) - (sqrt($eps) - $eps) / 2;
gelss($a, $b, $rcond, $s, ($rank=null),($info = null));

gelss ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. no cap.

gelsd

Signature: ([io]A(m,n); [io]B(p,q); rcond(); [o]s(minmn); int [o]rank();int [o]info(); int [t]iwork(iworkn))

Computes the minimum-norm solution to a real linear least squares problem:

minimize 2-norm(| b - A*x |)

using the singular value decomposition (SVD) of A. A is an M-by-N matrix which may be rank-deficient. (rank-deficient is valid and we support her 🫶)

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X.

The problem is solved in three steps:

Reduce the coefficient matrix A to bidiagonal form with Householder transformations, reducing the original problem into a "bidiagonal least squares problem" (BLS)
Solve the BLS using a divide and conquer approach.
Apply back all the Householder transformations to solve the original least squares problem.

The effective rank of A is determined by treating as zero those singular values which are less than rcond times the largest singular value.

Arguments

A — On entry, the M-by-N matrix A. On exit, A has been destroyed. 😭 (that's just how it is)
B — On entry, the M-by-NRHS right hand side matrix B. On exit, B is overwritten by the N-by-NRHS solution matrix X. If m >= n and rank = n, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements n+1:m in that column. The leading dimension of the array B >= max(1,M,N).
s — The singular values of A in decreasing order. The condition number of A in the 2-norm = s(1)/s(min(m,n)).
rcond — rcond is used to determine the effective rank of A. Singular values s(i) <= rcond*s(1) are treated as zero. If rcond < 0, machine precision is used instead.
rank — The effective rank of A, i.e., the number of singular values which are greater than rcond*s(1).
info — bestie check this before you trust the answer:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value. 🚩
- > 0: the algorithm for computing the SVD failed to converge; if info = i, i off-diagonal elements of an intermediate bidiagonal form did not converge to zero. 💀

$a= random(7,5);
# $b will contain X
# TODO better example with slice
$b = random(7,6);
$eps = lamch(0);
$s =zeroes(5);
#Threshold for rank estimation
$rcond = sqrt($eps) - (sqrt($eps) - $eps) / 2;
gelsd($a, $b, $rcond, $s, ($rank=null),($info = null));

gelsd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

gglse

Signature: ([phys]A(m,n); [phys]B(p,n);[io,phys]c(m);[phys]d(p);[o,phys]x(n);int [o,phys]info())

Solves the linear equality-constrained least squares (LSE) problem:

minimize || c - A*x ||_2   subject to   B*x = d

where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
M-vector, and d is a given P-vector. It is assumed that
P <= N <= M+P, and

     rank(B) = P and  rank( ( A ) ) = N.
                          ( ( B ) )

These conditions ensure that the LSE problem has a unique solution, which is obtained using a GRQ factorization of the matrices B and A. one answer only, she's monogamous like that 💅

Arguments

A — On entry, the M-by-N matrix A. On exit, A is destroyed.
B — On entry, the P-by-N matrix B. On exit, B is destroyed.
c — On entry, c contains the right hand side vector for the least squares part of the LSE problem. On exit, the residual sum of squares for the solution is given by the sum of squares of elements N-P+1 to M of vector c.
d — On entry, d contains the right hand side vector for the constrained equation. On exit, d is destroyed.
x — On exit, x is the solution of the LSE problem. 👑
info — the vibe check:
- = 0: successful exit. ✨
- < 0: if info = -i, the i-th argument had an illegal value. 🚩

$a = random(7,5);
$b = random(4,5);
$c = random(7);
$d = random(4);
$x = zeroes(5);
gglse($a, $b, $c, $d, $x, ($info=null));

gglse ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

ggglm

Signature: ([phys]A(n,m); [phys]B(n,p);[phys]d(n);[o,phys]x(m);[o,phys]y(p);int [o,phys]info())

Solves a general Gauss-Markov linear model (GLM) problem:

minimize || y ||_2   subject to   d = A*x + B*y
   x

where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
given N-vector. It is assumed that M <= N <= M+P, and

       rank(A) = M    and    rank( A B ) = N.

Under these assumptions, the constrained equation is always consistent, and there is a unique solution x and a minimal 2-norm solution y, which is obtained using a generalized QR factorization of A and B.

In particular, if matrix B is square nonsingular, then the problem GLM is equivalent to the following weighted linear least squares problem

minimize || inv(B)*(d-A*x) ||_2
   x

where inv(B) denotes the inverse of B.

Arguments

A — On entry, the N-by-M matrix A. On exit, A is destroyed.
B — On entry, the N-by-P matrix B. On exit, B is destroyed.
d — On entry, d is the left hand side of the GLM equation. On exit, d is destroyed.
x — (paired with y below 👇)
y — On exit, x and y are the solutions of the GLM problem. 🫶
info — checking in:
- = 0: successful exit. ✨
- < 0: if info = -i, the i-th argument had an illegal value. 🚩

$a = random(7,5);
$b = random(7,4);
$d = random(7);
$x = zeroes(5);
$y = zeroes(4);
ggglm($a, $b, $d, $x, $y,($info=null));

ggglm ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

getrf

Signature: ([io]A(m,n); int [o]ipiv(p); int [o]info())

Computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges.

The factorization has the form

A = P * L * U

where P is a permutation matrix, L is lower triangular with unit
diagonal elements (lower trapezoidal if m > n), and U is upper
triangular (upper trapezoidal if m < n).

This is the right-looking Level 3 BLAS version of the algorithm. highkey the fancy one 🔥

Arguments

A — On entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.
ipiv — The pivot indices; for 1 <= i <= min(M,N), row i of the matrix was interchanged with row ipiv(i).
info — the tea:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value 🚩
- > 0: if info = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations. 💀

$a = random (float, 100,50);
$ipiv = zeroes(long, 50);
$info = null;
getrf($a, $ipiv, $info);

getrf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

getf2

Signature: ([io]A(m,n); int [o]ipiv(p); int [o]info())

Computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges.

The factorization has the form

A = P * L * U

where P is a permutation matrix, L is lower triangular with unit
diagonal elements (lower trapezoidal if m > n), and U is upper
triangular (upper trapezoidal if m < n).

This is the right-looking Level 2 BLAS version of the algorithm.

📋 Arguments 📋

A — on entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored (they're literally all 1, we're not wasting memory on that 💅).
ipiv — the pivot indices; for 1 <= i <= min(M,N), row i of the matrix was interchanged with row ipiv(i).
info — the vibe report:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations. 💀

$a = random (float, 100,50);
$ipiv = zeroes(long, 50);
$info = null;
getf2($a, $ipiv, $info);

getf2 ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

🦋 sytrf 🦋

Signature: ([io,phys]A(n,n); int uplo(); int [o,phys]ipiv(n); int [o,phys]info())

Computes the factorization of a real symmetric matrix A using the Bunch-Kaufman diagonal pivoting method. The form of the factorization is

A = U*D*U'  or  A = L*D*L'
where U (or L) is a product of permutation and unit upper (lower)
triangular matrices, and D is symmetric and block diagonal with
1-by-1 and 2-by-2 diagonal blocks.

This is the blocked version of the algorithm, calling Level 3 BLAS. no cap.

📋 Arguments 📋

uplo —
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced.

On exit, the block diagonal matrix D and the multipliers used to obtain the factor U or L (see below for further details 👀).

ipiv — details of the interchanges and the block structure of D.
- If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged and D(k,k) is a 1-by-1 diagonal block.
- If uplo = 0 and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
- If uplo = 1 and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
info — the receipts:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, and division by zero will occur if it is used to solve a system of equations. 🚨

🔍 Further Details 🔍

If uplo = 0, then A = U*D*U', where

   U = P(n)*U(n)* ... *P(k)U(k)* ...,

i.e., U is a product of terms P(k)*U(k), where k decreases from n to 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as defined by ipiv(k), and U(k) is a unit upper triangular matrix, such that if the diagonal block D(k) is of order s (s = 1 or 2), then

           (   I    v    0   )   k-s
   U(k) =  (   0    I    0   )   s
           (   0    0    I   )   n-k
              k-s   s   n-k

If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), and A(k,k), and v overwrites A(1:k-2,k-1:k).

If uplo = 1, then A = L*D*L', where

   L = P(1)*L(1)* ... *P(k)*L(k)* ...,

i.e., L is a product of terms P(k)*L(k), where k increases from 1 to n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as defined by ipiv(k), and L(k) is a unit lower triangular matrix, such that if the diagonal block D(k) is of order s (s = 1 or 2), then

           (   I    0     0   )  k-1
   L(k) =  (   0    I     0   )  s
           (   0    v     I   )  n-k-s+1
              k-1   s  n-k-s+1

If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).

$a = random(100,100);
$ipiv = zeroes(100);
$info = null;
# Assume $a is symmetric
sytrf($a, 0, $ipiv, $info);

sytrf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

🌸 sytf2 🌸

Signature: ([io,phys]A(n,n); int uplo(); int [o,phys]ipiv(n); int [o,phys]info())

Computes the factorization of a real symmetric matrix A using the Bunch-Kaufman diagonal pivoting method. The form of the factorization is

A = U*D*U'  or  A = L*D*L'
where U (or L) is a product of permutation and unit upper (lower)
triangular matrices, and D is symmetric and block diagonal with
1-by-1 and 2-by-2 diagonal blocks.

This is the unblocked version of the algorithm, calling Level 2 BLAS. it's giving smaller-scale but still iconic 🎀

📋 Arguments 📋

uplo —
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced.

On exit, the block diagonal matrix D and the multipliers used to obtain the factor U or L (see below for further details).

ipiv — details of the interchanges and the block structure of D.
- If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged and D(k,k) is a 1-by-1 diagonal block.
- If uplo = 0 and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
- If uplo = 1 and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
info — the receipts:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, and division by zero will occur if it is used to solve a system of equations. 💀

For further details see sytrf.

$a = random(100,100);
$ipiv = zeroes(100);
$info = null;
# Assume $a is symmetric
sytf2($a, 0, $ipiv, $info);

sytf2 ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

👑 potrf 👑

Signature: ([io,phys]A(n,n); int uplo(); int [o,phys]info())

Computes the Cholesky factorization of a real symmetric positive definite matrix A.

The factorization has the form

A = U' * U,  if uplo = 0, or
A = L  * L',  if uplo = 1,
where U is an upper triangular matrix and L is lower triangular.

This is the block version of the algorithm, calling Level 3 BLAS. main character energy fr fr 🔥

📋 Arguments 📋

uplo —
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced.

On exit, if info = 0, the factor U or L from the Cholesky factorization A = U'*U or A = L*L'.

info — the receipts:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, the leading minor of order i is not positive definite, and the factorization could not be completed. 😭

$a = random(100,100);
# Assume $a is symmetric positive definite
potrf($a, 0, ($info = null));

potrf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

🧚 potf2 🧚

Signature: ([io,phys]A(n,n); int uplo(); int [o,phys]info())

Computes the Cholesky factorization of a real symmetric positive definite matrix A.

The factorization has the form

A = U' * U,  if uplo = 0, or
A = L  * L',  if uplo = 1,
where U is an upper triangular matrix and L is lower triangular.

This is the unblocked version of the algorithm, calling Level 2 BLAS.

📋 Arguments 📋

uplo —
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — on entry, the symmetric matrix A. If uplo = 0, the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced.

On exit, if info = 0, the factor U or L from the Cholesky factorization A = U'*U or A = L*L'.

info — the vibe check on how it went, bestie:
- = 0: successful exit, we ate 💅
- < 0: if info = -i, then the i-th argument had an illegal value (like, no)
- > 0: if info = i, the leading minor of order i is not positive definite, and the factorization could not be completed 😭

$a = random(100,100);
# Assume $a is symmetric positive definite
potf2($a, 0, ($info = null));

potf2 ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💖

getri

Signature: ([io,phys]A(n,n); int [phys]ipiv(n); int [o,phys]info())

Okay so this one computes the inverse of a matrix using the LU factorization computed by getrf. 🦋

It inverts U and then computes inv(A) by solving the system

inv(A)*L = inv(U) for inv(A).

Arguments

A — On entry: the factors L and U from the factorization A = P*L*U as computed by getrf. On exit: if info = 0, the inverse of the original matrix A. ✨
ipiv — The pivot indices from getrf; for 1<=i<=N, row i of the matrix was interchanged with row ipiv(i).
info — the group chat update:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, U(i,i) is exactly zero; the matrix is singular and its inverse could not be computed 💀

$a = random (float, 100, 100);
$ipiv = zeroes(long, 100);
$info = null;
getrf($a, $ipiv, $info);
if ($info == 0){
       getri($a, $ipiv, $info);
}
print "Inverse of \$a is :\n $a" unless $info;

getri ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🌸

sytri

Signature: ([io]A(n,n); int uplo(); int ipiv(n); int [o]info(); [t]work(n))

Computes the inverse of a real symmetric indefinite matrix A using the factorization A = U*D*U' or A = L*D*L' computed by sytrf. 💜

Arguments

uplo — Specifies whether the details of the factorization are stored as an upper or lower triangular matrix.
- = 0: Upper triangular, form is A = U*D*U';
- = 1: Lower triangular, form is A = L*D*L'.
A — On entry: the block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by sytrf. On exit: if info = 0, the (symmetric) inverse of the original matrix. If uplo = 0, the upper triangular part of the inverse is formed and the part of A below the diagonal is not referenced; if uplo = 1 the lower triangular part of the inverse is formed and the part of A above the diagonal is not referenced. 👑
ipiv — Details of the interchanges and the block structure of D as determined by sytrf.
info —
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, D(i,i) = 0; the matrix is singular and its inverse could not be computed 🙈

$a = random (float, 100, 100);
# assume $a is symmetric
$ipiv = zeroes(long, 100);
sytrf($a, 0, $ipiv, ($info=null));
if ($info == 0){
       sytri($a, 0, $ipiv, $info);
}
print "Inverse of \$a is :\n $a" unless $info;

sytri ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🎀

potri

Signature: ([io,phys]A(n,n); int uplo(); int [o,phys]info())

Computes the inverse of a real symmetric positive definite matrix A using the Cholesky factorization A = U'*U or A = L*L' computed by potrf. 🌈

Arguments

uplo —
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — On entry: the triangular factor U or L from the Cholesky factorization A = U'*U or A = L*L', as computed by potrf. On exit: the upper or lower triangle of the (symmetric) inverse of A, overwriting the input factor U or L. 💫
info —
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, the (i,i) element of the factor U or L is zero, and the inverse could not be computed 😭

$a = random (float, 100, 100);
# Assume $a is symmetric positive definite
potrf($a, 0, ($info = null));
if ($info == 0){ # Hum... is it positive definite????
       potri($a, 0,$info);
}
print "Inverse of \$a is :\n $a" unless $info;

potri ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. ⭐

trtri

Signature: ([io,phys]A(n,n); int uplo(); int diag(); int [o,phys]info())

Computes the inverse of a real upper or lower triangular matrix A. 🔥

This is the Level 3 BLAS version of the algorithm.

Arguments

uplo —
- = 0: A is upper triangular;
- = 1: A is lower triangular.
diag —
- = 0: A is non-unit triangular;
- = 1: A is unit triangular.
A — On entry: the triangular matrix A. If uplo = 0, the leading N-by-N upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If diag = 1, the diagonal elements of A are also not referenced and are assumed to be 1. On exit: the (triangular) inverse of the original matrix, in the same storage format. ✨
info —
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, A(i,i) is exactly zero. The triangular matrix is singular and its inverse can not be computed 💀

$a = random (float, 100, 100);
# assume $a is upper triangular
trtri($a, 1, ($info=null));
print "Inverse of \$a is :\n transpose($a)" unless $info;

trtri ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🧚

trti2

Signature: ([io,phys]A(n,n); int uplo(); int diag(); int [o,phys]info())

Computes the inverse of a real upper or lower triangular matrix A. 🦋

This is the Level 2 BLAS version of the algorithm.

Arguments

uplo —
- = 0: A is upper triangular;
- = 1: A is lower triangular.
diag —
- = 0: A is non-unit triangular;
- = 1: A is unit triangular.
A — On entry: the triangular matrix A. If uplo = 0, the leading N-by-N upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If diag = 1, the diagonal elements of A are also not referenced and are assumed to be 1. On exit: the (triangular) inverse of the original matrix, in the same storage format. 💜
info —
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

$a = random (float, 100, 100);
# assume $a is upper triangular
trtri2($a, 1, ($info=null));
print "Inverse of \$a is :\n transpose($a)" unless $info;

trti2 ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🌸

getrs

Signature: ([phys]A(n,n); int trans(); [io,phys]B(n,m); int [phys]ipiv(n); int [o,phys]info())

Solves a system of linear equations

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

with a general N-by-N matrix A using the LU factorization computed by getrf. 🔥

Arguments

trans — Specifies the form of the system of equations:
- = 0: A * X = B (No transpose)
- = 1: A'* X = B (Transpose)
A — The factors L and U from the factorization A = P*L*U as computed by getrf.
ipiv — The pivot indices from getrf; for 1<=i<=N, row i of the matrix was interchanged with row ipiv(i).
B — On entry: the right hand side matrix B. On exit: the solution matrix X. ✨
info —
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

$a = random (float, 100, 100);
$ipiv = zeroes(long, 100);
$b = random(100,50);
getrf($a, $ipiv, ($info=null));
if ($info == 0){
       getrs($a, 0, $b, $ipiv, $info);
}
print "X is :\n $b" unless $info;

getrs ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💅

sytrs

Signature: ([phys]A(n,n); int uplo();[io,phys]B(n,m); int [phys]ipiv(n); int [o,phys]info())

Solves a system of linear equations A*X = B with a real symmetric matrix A using the factorization A = U*D*U' or A = L*D*L' computed by sytrf. 🎀

Arguments

uplo — Specifies whether the details of the factorization are stored as an upper or lower triangular matrix.
- = 0: Upper triangular, form is A = U*D*U';
- = 1: Lower triangular, form is A = L*D*L'.
A — The block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by sytrf.
ipiv — Details of the interchanges and the block structure of D as determined by sytrf.
B — On entry: the right hand side matrix B. On exit: the solution matrix X. 👑
info —
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

$a = random (float, 100, 100);
$b = random(50,100);
$a = transpose($a);
$b = transpose($b);
# Assume $a is symmetric
sytrf($a, 0, ($ipiv=zeroes(100)), ($info=null));
if ($info == 0){
       sytrs($a, 0, $b, $ipiv, $info);
}
print("X is :\n".transpose($b))unless $info;

okay so real quick bestie — sytrs literally ignores the bad-value flag of the input ndarrays 💅. But it will set the bad-value flag on all the output ndarrays if the flag is set for any of the input ndarrays. no cap, that's the whole deal.

🌸 potrs 🌸

Signature: ([phys]A(n,n); int uplo(); [io,phys]B(n,m); int [o,phys]info())

Solves a system of linear equations A*X = B where A is a symmetric positive definite matrix, using the Cholesky factorization A = U'*U or A = L*L' that potrf already computed for you. it's giving "reuse the homework you already did" ✨

Arguments

uplo — which triangle of A we're actually storing:
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — the triangular factor U or L from the Cholesky factorization A = U'*U or A = L*L', as computed by potrf. yes, the same one. don't go making a new one 🙈
B — on entry, the right hand side matrix B. on exit, the solution matrix X. glow-up arc fr fr 🦋
info — the vibe check:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value 💀

$a = random (float, 100, 100);
$b = random(50,100);
$a = transpose($a);
$b = transpose($b);
# Assume $a is symmetric positive definite
potrf($a, 0, ($info=null));
if ($info == 0){
       potrs($a, 0, $b, $info);
}
print("X is :\n".transpose($b))unless $info;

potrs ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. same energy as always, you know the drill 💜

🎀 trtrs 🎀

Signature: ([phys]A(n,n); int uplo(); int trans(); int diag();[io,phys]B(n,m); int [o,phys]info())

Solves a triangular system of the form

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

        where A is a triangular matrix of order N, and B is an N-by-NRHS
        matrix.

and bestie, it literally checks that A is nonsingular. it does the work so you don't have to 🙏

Arguments

uplo — which triangle A is:
- = 0: A is upper triangular;
- = 1: A is lower triangular.
trans — what form of the system we're solving:
- = 0: A * X = B (No transpose)
- = 1: A**T * X = B (Transpose)
diag — is the diagonal doing anything or not:
- = 0: A is non-unit triangular;
- = 1: A is unit triangular.
A — the triangular matrix A. if uplo = 0, the leading N-by-N upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. if uplo = 1, the leading N-by-N lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. and if diag = 1, the diagonal elements of A are also not referenced and are assumed to be 1. iconic, honestly — half the array is just vibing in the background 🧚
B — on entry, the right hand side matrix B. on exit, if info = 0, the solution matrix X.
info — the receipts:
- = 0: successful exit ✨
- < 0: if info = -i, the i-th argument had an illegal value 💀
- > 0: if info = i, the i-th diagonal element of A is zero, indicating that the matrix is singular and the solutions X have not been computed. so like, no solution, no X, nothing. tragic. 😭

# Assume $a is upper triangular
$a = random (float, 100, 100);
$b = random(50,100);
$a = transpose($a);
$b = transpose($b);
$info = null;
trtrs($a, 0, 0, 0, $b, $info);
print("X is :\n".transpose($b))unless $info;

trtrs ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💜

🦋 latrs 🦋

Signature: ([phys]A(n,n); int uplo(); int trans(); int diag(); int normin();[io,phys]x(n); [o,phys]scale();[io,phys]cnorm(n);int [o,phys]info())

Solves one of the triangular systems

        A *x = s*b  or  A'*x = s*b

with scaling to prevent overflow. here A is an upper or lower triangular matrix, A' denotes the transpose of A, x and b are n-element vectors, and s is a scaling factor — usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. if the unscaled problem won't cause overflow, the Level 2 BLAS routine trsv gets called. and if the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned. main character move honestly 👑

Further Details

a rough bound on x is computed; if that's less than overflow, trsv is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation. vigilance 💅

a columnwise scheme is used for solving A*x = b. the basic algorithm if A is lower triangular is

         x[1:n] := b[1:n]
         for j = 1, ..., n
              x(j) := x(j) / A(j,j)
              x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
         end

define bounds on the components of x after j iterations of the loop:

       M(j) = bound on x[1:j]
       G(j) = bound on x[j+1:n]

initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

then for iteration j+1 we have

       M(j+1) <= G(j) / | A(j+1,j+1) |
       G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
              <= G(j) ( 1 + cnorm(j+1) / | A(j+1,j+1) | )

where cnorm(j+1) is greater than or equal to the infinity-norm of column j+1 of A, not counting the diagonal.

hence

       G(j) <= G(0) product ( 1 + cnorm(i) / | A(i,i) | )
                    1<=i<=j
    and

       |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + cnorm(i) / |A(i,i)| )
                                     1<=i< j

since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).

the bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. if the computed bound is greater than a large constant, x is scaled to prevent overflow — but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. so yeah, it doesn't just give up, it finds you a whole other answer. 🌈

similarly, a row-wise scheme is used to solve A'*x = b. the basic algorithm for A upper triangular is

         for j = 1, ..., n
              x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
         end

we simultaneously compute two bounds

         G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
         M(j) = bound on x(i), 1<=i<=j

the initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. then the bound on x(j) is

         M(j) <= M(j-1) * ( 1 + cnorm(j) ) / | A(j,j) |

              <= M(0) * product ( ( 1 + cnorm(i) ) / |A(i,i)| )
                        1<=i<=j

and we can safely call trsv if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow). ✨

Arguments

uplo — specifies whether the matrix A is upper or lower triangular.
- = 0: Upper triangular
- = 1: Lower triangular
trans — specifies the operation applied to A.
- = 0: Solve A * x = s*b (No transpose)
- = 1: Solve A'* x = s*b (Transpose)
diag — specifies whether or not the matrix A is unit triangular.
- = 0: Non-unit triangular
- = 1: Unit triangular
normin — specifies whether cnorm has been set or not.
- = 1: cnorm contains the column norms on entry
- = 0: cnorm is not set on entry. on exit, the norms will be computed and stored in cnorm.
A — the triangular matrix A. if uplo = 0, the leading n by n upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. if uplo = 1, the leading n by n lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. if diag = 1, the diagonal elements of A are also not referenced and are assumed to be 1.
x — on entry, the right hand side b of the triangular system. on exit, x is overwritten by the solution vector x. bye b, hello x 💅
scale — the scaling factor s for the triangular system A * x = s*b or A'* x = s*b. if scale = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0. that's a whole situation, so watch for it.
cnorm — if normin = 0, cnorm is an output argument and cnorm(j) returns the 1-norm of the offdiagonal part of the j-th column of A. if normin = 1, cnorm is an input argument and cnorm(j) contains the norm of the off-diagonal part of the j-th column of A. if trans = 0, cnorm(j) must be greater than or equal to the infinity-norm, and if trans = 1, cnorm(j) must be greater than or equal to the 1-norm. input or output depending on normin, so like, read the room. 🙈
info — the classic:
- = 0: successful exit 🎉
- < 0: if info = -k, the k-th argument had an illegal value 💀

# Assume $a is upper triangular
$a = random (float, 100, 100);
$b = random(100);
$a = transpose($a);
$info = null;
$scale= null;
$cnorm = zeroes(100);
latrs($a, 0, 0, 0, 0,$b, $scale, $cnorm,$info);

latrs ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

💖 gecon 💖

Signature: (A(n,n); int norm(); anorm(); [o]rcond();int [o]info(); int [t]iwork(n); [t]work(workn))

Estimates the reciprocal of the condition number of a general real matrix A, in either the 1-norm or the infinity-norm, using the LU factorization computed by getrf.

an estimate is obtained for norm(inv(A)), and the reciprocal of the condition number is computed as

       rcond = 1 / ( norm(A) * norm(inv(A)) ).

it's giving "how sketchy is this matrix, on a scale" and honestly we love a matrix that knows itself 🎀

Arguments

norm — specifies whether the 1-norm condition number or the infinity-norm condition number is required:
- = 0: Infinity-norm.
- = 1: 1-norm;
A — the factors L and U from the factorization A = P*L*U as computed by getrf.
anorm — if norm = 0, the infinity-norm of the original matrix A. if norm = 1, the 1-norm of the original matrix A. note: the original matrix, not the factors, iykyk
rcond — the reciprocal of the condition number of the matrix A, computed as rcond = 1/(norm(A) * norm(inv(A))). the moment of truth ✨
info — you already know:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value 💀

$a = random (float, 100, 100);
$anorm  = $a->lange(1);
$ipiv = zeroes(long, 100);
$info = null;
getrf($a, $ipiv, $info);
($rcond, $info) = gecon($a, 1, $anorm) unless $info != 0;

and that's where this part lands bestie 🌈💖

Also, real quick, one last thing about <code>gecon</code>: it ignores the bad-value flag of the input ndarrays 🙈. But it will set the bad-value flag on all the output ndarrays if the flag is set for any of the input ndarrays. That's just how it do 💅

🌸 sycon 🌸

Signature: ([phys]A(n,n); int uplo(); int ipiv(n); [phys]anorm(); [o,phys]rcond();int [o,phys]info(); int [t]iwork(n); [t]work(workn))

So what it does is estimate the reciprocal of the condition number (in the 1-norm, bestie) of a real symmetric matrix A, using the factorization A = U*D*U' or A = L*D*L' that sytrf already computed for you ✨

It gets an estimate for norm(inv(A)), and then the reciprocal of the condition number is literally just rcond = 1 / (anorm * norm(inv(A))). Math girl summer 💖

💫 Arguments 💫

uplo — tells you whether the factorization details are stored as an upper or lower triangular matrix.
- = 0: Upper triangular, form is A = U*D*U';
- = 1: Lower triangular, form is A = L*D*L'.
A — the block diagonal matrix D and the multipliers used to get the factor U or L, exactly as computed by sytrf.
ipiv — the details of the interchanges and the block structure of D, as determined by sytrf.
anorm — the 1-norm of the original matrix A. iconic.
rcond — the reciprocal of the condition number of the matrix A, computed as rcond = 1/(anorm * aimvnm), where ainvnm is an estimate of the 1-norm of inv(A) computed in this routine.
info — the exit status:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value.

# Assume $a is symmetric
$a = random (float, 100, 100);
$anorm  = $a->lansy(1,1);
$ipiv = zeroes(long, 100);
$info = null;
sytrf($a, 1,$ipiv, $info);
($rcond, $info) = sycon($a, 1, $anorm) unless $info != 0;

sycon ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. no cap 🫶

🎀 pocon 🎀

Signature: (A(n,n); int uplo(); anorm(); [o]rcond();int [o]info(); int [t]iwork(n); [t]work(workn))

Estimates the reciprocal of the condition number (again in the 1-norm, we stay consistent) of a real symmetric positive definite matrix, using the Cholesky factorization A = U'*U or A = L*L' that potrf computed. obsessed with her honestly 💖

An estimate is obtained for norm(inv(A)), and the reciprocal of the condition number comes out as rcond = 1 / (anorm * norm(inv(A))). periodt.

💫 Arguments 💫

uplo — which triangle we're working with:
- = 0: Upper triangle of A is stored;
- = 1: Lower triangle of A is stored.
A — the triangular factor U or L from the Cholesky factorization A = U'*U or A = L*L', as computed by potrf.
anorm — the 1-norm of the matrix A.
rcond — the reciprocal of the condition number of the matrix A, computed as rcond = 1/(anorm * ainvnm), where ainvnm is an estimate of the 1-norm of inv(A) computed in this routine.
info — the exit status:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value

# Assume $a is symmetric positive definite
$a = random (float, 100, 100);
$anorm  = $a->lansy(1,1);
$info = null;
potrf($a,  0, $info);
($rcond, $info) = pocon($a, 1, $anorm) unless $info != 0;

pocon ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💅

🦋 trcon 🦋

Signature: (A(n,n); int norm();int uplo();int diag(); [o]rcond();int [o]info(); int [t]iwork(n); [t]work(workn))

Estimates the reciprocal of the condition number of a triangular matrix A — and you get to pick, either the 1-norm or the infinity-norm, main character choice fr fr ✨

The norm of A gets computed, an estimate is obtained for norm(inv(A)), and then the reciprocal of the condition number is computed as

        rcond = 1 / ( norm(A) * norm(inv(A)) ).

💫 Arguments 💫

norm — picks whether you want the 1-norm condition number or the infinity-norm condition number:
- = 0: Infinity-norm.
- = 1: 1-norm;
uplo — where the triangle lives:
- = 0: A is upper triangular;
- = 1: A is lower triangular.
diag — whether the diagonal is part of the deal:
- = 0: A is non-unit triangular;
- = 1: A is unit triangular.
A — the triangular matrix A. If uplo = 0, the leading N-by-N upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If uplo = 1, the leading N-by-N lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If diag = 1, the diagonal elements of A are also not referenced and are assumed to be 1. lowkey so much going on but we love her 🫶
rcond — the reciprocal of the condition number of the matrix A, computed as rcond = 1/(norm(A) * norm(inv(A))).
info — the exit status:
- = 0: successful exit 🎉
- < 0: if info = -i, the i-th argument had an illegal value

# Assume $a is upper triangular
$a = random (float, 100, 100);
$info = null;
($rcond, $info) = trcon($a, 1, 1, 0) unless $info != 0;

trcon ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💜

⭐ geqp3 ⭐

Signature: ([io,phys]A(m,n); int [io,phys]jpvt(n); [o,phys]tau(k); int [o,phys]info())

geqp3 computes a QR factorization using Level 3 BLAS with column pivoting of a matrix A:

                A*P = Q*R

The matrix Q is represented as a product of elementary reflectors

        Q = H(1) H(2) . . . H(k), where k = min(m,n).

Each H(i) has the form

        H(i) = I - tau * v * v'

        where tau is a real/complex scalar, and v is a real/complex vector
        with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
        A(i+1:m,i), and tau in tau(i).

💫 Arguments 💫

A — on entry, the M-by-N matrix A. On exit, the upper triangle of the array contains the min(M,N)-by-N upper trapezoidal matrix R; the elements below the diagonal, together with the array tau, represent the orthogonal matrix Q as a product of min(M,N) elementary reflectors. glow up ✨
jpvt — on entry, if jpvt(J)!=0, the J-th column of A is permuted to the front of A*P (a leading column); if jpvt(J)=0, the J-th column of A is a free column. On exit, if jpvt(J)=K, then the J-th column of A*P was the the K-th column of A. (yes, "the the" — we keep it exactly as written, iykyk 🙈)
tau — the scalar factors of the elementary reflectors.
info — the exit status:
- = 0: successful exit.
- < 0: if info = -i, the i-th argument had an illegal value.

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
$jpvt = zeroes(long, 50);
geqp3($a, $jpvt, $tau, $info);

geqp3 ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💅

🌈 geqrf 🌈

Signature: ([io,phys]A(m,n); [o,phys]tau(k); int [o,phys]info())

geqrf computes a QR factorization of a matrix A:

        A = Q * R

The matrix Q is represented as a product of elementary reflectors

        Q = H(1) H(2) . . . H(k), where k = min(m,n).

Each H(i) has the form

        H(i) = I - tau * v * v'

        where tau is a real/complex scalar, and v is a real/complex vector
        with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
        A(i+1:m,i), and tau in tau(i).

💫 Arguments 💫

A — on exit, the elements on and above the diagonal of the array contain the min(M,N)-by-N upper trapezoidal matrix R (R is upper triangular if m >= n); the elements below the diagonal, with the array TAU, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors. she ate 🔥
tau — the scalar factors of the elementary reflectors.
info — the exit status:
- = 0: successful exit.
- < 0: if info = -i, the i-th argument had an illegal value.

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
geqrf($a, $tau, $info);

geqrf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🧚 orgqr 🧚

Signature: ([io,phys]A(m,n); [phys]tau(k); int [o,phys]info())

Generates an M-by-N real matrix Q with orthonormal columns, defined as the first N columns of a product of K elementary reflectors of order M

        Q  =  H(1) H(2) . . . H(k)

        as returned by geqrf or geqp3.

💫 Arguments 💫

A — on entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by geqrf or geqp3 in the first k columns of its array argument A. On exit, the M-by-N matrix Q. the glow up is real 💖
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by geqrf or geqp3.
info — the exit status:
- = 0: successful exit
- < 0: if info = -i, the i-th argument has an illegal value

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
geqrf($a, $tau, $info);
orgqr($a, $tau, $info) unless $info != 0;

orgqr ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. ✨

👑 ormqr 👑

Signature: ([phys]A(p,k); int side(); int trans(); [phys]tau(k); [io,phys]C(m,n);int [o,phys]info())

Overwrites the general real M-by-N matrix C with

                    side = 0     side = 1
trans = 0:      Q * C          C * Q
trans = 1:      Q' * C       C * Q'

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

        Q = H(1) H(2) . . . H(k)

        as returned by geqrf or geqp3.

Q is of order M if side = 0 and of order N if side = 1. know the difference bestie, it matters 💅

💫 Arguments 💫

side — which side we're applying from:
- = 0: apply Q or Q' from the Left;
- = 1: apply Q or Q' from the Right.
trans — transpose or nah:
- = 0: No transpose, apply Q;
- = 1: Transpose, apply Q'.
A — the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by geqrf or geqp3 in the first k columns of its array argument A. A is modified by the routine but restored on exit. she's changed but she comes back, growth 🌱
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by geqrf or geqp3.
C — ok so on the way in, this is your M-by-N matrix C 🫶. On the way out? It's been fully overwritten by Q*C or Q'*C or C*Q' or C*Q. No takebacks, bestie, C is a ✨changed matrix✨ now.
info — the receipts, pls read them:
- = 0: successful exit 💖
- < 0: if info = -i, the i-th argument had an illegal value 💀

Like, here's the whole vibe in code, bestie:

$a = random (float, 50, 100);
$a = transpose($a);
$info = null;
$tau = zeroes(float, 50);
geqrf($a, $tau, $info);
$c = random(70,50);
# $c will contain the result
$c->reshape(70,100);
$c = transpose($c);
ormqr($a, $tau, $c, $info);

And one thing you gotta know: ormqr ignores the bad-value flag of the input ndarrays. But it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. So it's giving "I don't check, but I do tell" energy, lowkey iconic 🙈

🌸 `gelqf` 🌸

Signature: ([io,phys]A(m,n); [o,phys]tau(k); int [o,phys]info())

Bestie, this one computes an LQ factorization of a real M-by-N matrix A:

        A = L * Q.

The matrix Q is represented as a product of elementary reflectors

       Q = H(k) . . . H(2) H(1), where k = min(m,n).

Each H(i) has the form

        H(i) = I - tau * v * v'

where tau is a real scalar, and v is a real vector with v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), and tau in tau(i). Periodt 💅

📋 Arguments

A — on entry, the M-by-N matrix A. On exit, the elements on and below the diagonal of the array contain the m-by-min(m,n) lower trapezoidal matrix L (L is lower triangular if m <= n); the elements above the diagonal, with the array tau, represent the orthogonal matrix Q as a product of elementary reflectors. So basically A gets a ✨glow-up✨.
tau — the scalar factors of the elementary reflectors.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

Here's how you'd actually call it, no cap:

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
gelqf($a, $tau, $info);

gelqf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. Same energy as before, fr fr 🦋

🌈 `orglq` 🌈

Signature: ([io,phys]A(m,n); [phys]tau(k); int [o,phys]info())

This one generates an M-by-N real matrix Q with orthonormal rows 👑, defined as the first M rows of a product of K elementary reflectors of order N

        Q  =  H(k) . . . H(2) H(1)

        as returned by gelqf.

📋 Arguments

A — on entry, the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by gelqf in the first k rows of its array argument A. On exit, the M-by-N matrix Q. Main character transformation, honestly ✨
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by gelqf.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument has an illegal value

Watch how these two team up, it's giving collab:

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
gelqf($a, $tau, $info);
orglq($a, $tau, $info) unless $info != 0;

orglq ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💜

🎀 `ormlq` 🎀

Signature: ([phys]A(k,p); int side(); int trans(); [phys]tau(k); [io,phys]C(m,n);int [o,phys]info())

Overwrites the general real M-by-N matrix C with

                    side = 0     side = 1
    trans = 0:      Q * C          C * Q
    trans = 1:      Q' * C       C * Q'

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

          Q = H(k) . . . H(2) H(1)

as returned by gelqf.

Q is of order M if "side" = 0 and of order N if "side" = 1. Slay 🌟

📋 Arguments

side — which way we're applying it:
- = 0: apply Q or Q' from the Left;
- = 1: apply Q or Q' from the Right.
trans — transpose or nah:
- = 0: No transpose, apply Q;
- = 1: Transpose, apply Q'.
A — the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by gelqf in the first k rows of its array argument A. Heads up tho: A is modified by the routine but restored on exit 🙈 so don't panic, it comes back.
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by gelqf.
C — on entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

Example time, bestie:

$a = random (float, 50, 100);
$a = transpose($a);
$info = null;
$tau = zeroes(float, 50);
gelqf($a, $tau, $info);
$c = random(70,50);
# $c will contain the result
$c->reshape(70,100);
$c = transpose($c);
ormlq($a, $tau, $c, $info);

ormlq ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🍓

⭐ `geqlf` ⭐

Signature: ([io,phys]A(m,n); [o,phys]tau(k); int [o,phys]info())

Computes a QL factorization of a real M-by-N matrix A:

        A = Q * L

The matrix Q is represented as a product of elementary reflectors

        Q = H(k) . . . H(2) H(1), where k = min(m,n).

Each H(i) has the form

        H(i) = I - tau * v * v'

where tau is a real scalar, and v is a real vector with v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in A(1:m-k+i-1,n-k+i), and tau in TAU(i). Yes, the capitalization switched, that's literally how it's written, don't come for me 💅

📋 Arguments

A — on entry, the M-by-N matrix A. On exit:
- if m >= n, the lower triangle of the subarray A(m-n+1:m,1:n) contains the N-by-N lower triangular matrix L;
- if m <= n, the elements on and below the (n-m)-th superdiagonal contain the M-by-N lower trapezoidal matrix L;
- the remaining elements, with the array tau, represent the orthogonal matrix Q as a product of elementary reflectors.
tau — the scalar factors of the elementary reflectors.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

Run it like this:

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
geqlf($a, $tau, $info);

geqlf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🌸

💖 `orgql` 💖

Signature: ([io,phys]A(m,n); [phys]tau(k); int [o,phys]info())

Generates an M-by-N real matrix Q with orthonormal columns, defined as the last N columns of a product of K elementary reflectors of order M

        Q  =  H(k) . . . H(2) H(1)

        as returned by geqlf.

📋 Arguments

A — on entry, the (n-k+i)-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by geqlf in the last k columns of its array argument A. On exit, the M-by-N matrix Q. Glow up, babes ✨
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by geqlf.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument has an illegal value

Bestie duo behaviour:

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
geqlf($a, $tau, $info);
orgql($a, $tau, $info) unless $info != 0;

orgql ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🧚

🔥 `ormql` 🔥

Signature: ([phys]A(p,k); int side(); int trans(); [phys]tau(k); [io,phys]C(m,n);int [o,phys]info())

Overwrites the general real M-by-N matrix C with

                    side = 0     side = 1
    trans = 0:      Q * C          C * Q
    trans = 1:      Q' * C       C * Q'

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

          Q = H(k) . . . H(2) H(1)

as returned by geqlf.

Q is of order M if "side" = 0 and of order N if "side" = 1. Same deal as ormlq, just a different parent routine, iykyk 🫶

📋 Arguments

side — pick a side:
- = 0: apply Q or Q' from the Left;
- = 1: apply Q or Q' from the Right.
trans — transpose check:
- = 0: No transpose, apply Q;
- = 1: Transpose, apply Q'.
A — the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by geqlf in the last k rows of its array argument A. And again: A is modified by the routine but restored on exit, so it's fine, it's fine 😭
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by geqlf.
C — on entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

Example incoming:

$a = random (float, 50, 100);
$a = transpose($a);
$info = null;
$tau = zeroes(float, 50);
geqlf($a, $tau, $info);
$c = random(70,50);
# $c will contain the result
$c->reshape(70,100);
$c = transpose($c);
ormql($a, $tau, $c, $info);

ormql ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💜

🦋 `gerqf` 🦋

Signature: ([io,phys]A(m,n); [o,phys]tau(k); int [o,phys]info())

Computes an RQ factorization of a real M-by-N matrix A:

        A = R * Q.

The matrix Q is represented as a product of elementary reflectors

        Q = H(1) H(2) . . . H(k), where k = min(m,n).

Notice the order flipped? Yeah. It's H(1) H(2) . . . H(k) this time, not the other way round. The way this routine said "we do things differently here" 😌

Each H(i) has the form

        H(i) = I - tau * v * v'

where tau is a real scalar, and v is a real vector with v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in A(m-k+i,1:n-k+i-1), and tau in TAU(i).

📋 Arguments

A — on entry, the M-by-N matrix A. On exit:
- if m <= n, the upper triangle of the subarray A(1:m,n-m+1:n) contains the M-by-M upper triangular matrix R;
- if m >= n, the elements on and above the (m-n)-th subdiagonal contain the M-by-N upper trapezoidal matrix R;
- the remaining elements, with the array tau, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors (see Further Details).
tau — the scalar factors of the elementary reflectors.
info — the receipts:
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
gerqf($a, $tau, $info);

and yeah — gerqf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💅

🦋 orgrq 🦋

Signature: ([io,phys]A(m,n); [phys]tau(k); int [o,phys]info())

okay so this one generates an M-by-N real matrix Q with orthonormal rows, and it's defined as the last M rows of a product of K elementary reflectors of order N, bestie:

        Q  =  H(1) H(2) . . . H(k)

as returned by gerqf. ✨

📋 Arguments 📋

A — On entry, the (m-k+i)-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by gerqf in the last k rows of its array argument A. On exit, the M-by-N matrix Q. 💖
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by gerqf.
info — = 0: successful exit. < 0: if info = -i, the i-th argument has an illegal value. 😭

$a = random (float, 100, 50);
$info = null;
$tau = zeroes(float, 50);
gerqf($a, $tau, $info);
orgrq($a, $tau, $info) unless $info != 0;

orgrq ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🌸 ormrq 🌸

Signature: ([phys]A(k,p); int side(); int trans(); [phys]tau(k); [io,phys]C(m,n);int [o,phys]info())

this one overwrites the general real M-by-N matrix C with — vibe check the little grid:

                           side = 0     side = 1
           trans = 0:      Q * C          C * Q
           trans = 1:      Q' * C       C * Q'

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

                 Q = H(1) H(2) . . . H(k)

as returned by gerqf. ✨

Q is of order M if "side" = 0 and of order N if "side" = 1. periodt.

📋 Arguments 📋

side — = 0: apply Q or Q' from the Left; = 1: apply Q or Q' from the Right. 👑
trans — = 0: No transpose, apply Q; = 1: Transpose, apply Q'. 💅
A — The i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by gerqf in the last k rows of its array argument A. A is modified by the routine but restored on exit. 🎀
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by gerqf.
C — On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. 🦋
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. 😭

$a = random (float, 50, 100);
$a = transpose($a);
$info = null;
$tau = zeroes(float, 50);
gerqf($a, $tau, $info);
$c = random(70,50);
# $c will contain the result
$c->reshape(70,100);
$c = transpose($c);
ormrq($a, $tau, $c, $info);

ormrq ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💖

⭐ tzrzf ⭐

Signature: ([io,phys]A(m,n); [o,phys]tau(k); int [o,phys]info())

so basically this one reduces the M-by-N ( M <= N ) real upper trapezoidal matrix A to upper triangular form by means of orthogonal transformations. 😌

The upper trapezoidal matrix A is factored as

        A = ( R  0 ) * Z,

where Z is an N-by-N orthogonal matrix and R is an M-by-M upper triangular matrix. slay.

The factorization is obtained by Householder's method. The kth transformation matrix, Z( k ), which is used to introduce zeros into the ( m - k + 1 )th row of A, is given in the form

       Z( k ) = ( I     0   ),
                ( 0  T( k ) )

where

       T( k ) = I - tau*u( k )*u( k )',   u( k ) = (   1    ),
                                                    (   0    )
                                                    ( z( k ) )

tau is a scalar and z( k ) is an ( n - m ) element vector. tau and z( k ) are chosen to annihilate the elements of the kth row of X. 👑

The scalar tau is returned in the kth element of "tau" and the vector u( k ) in the kth row of A, such that the elements of z( k ) are in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in the upper triangular part of A. 🎀

Z is given by

       Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).

📋 Arguments 📋

A — On entry, the leading M-by-N upper trapezoidal part of the array A must contain the matrix to be factorized. On exit, the leading M-by-M upper triangular part of A contains the upper triangular matrix R, and elements M+1 to N of the first M rows of A, with the array tau, represent the orthogonal matrix Z as a product of M elementary reflectors. 💜
tau — The scalar factors of the elementary reflectors.
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. 😭

$a = random (float, 50, 100);
$info = null;
$tau = zeroes(float, 50);
tzrzf($a, $tau, $info);

tzrzf ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🌈 ormrz 🌈

Signature: ([phys]A(k,p); int side(); int trans(); [phys]tau(k); [io,phys]C(m,n);int [o,phys]info())

okay bestie, this one overwrites the general real M-by-N matrix C with

                           side = 0     side = 1
           trans = 0:      Q * C          C * Q
           trans = 1:      Q' * C       C * Q'

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

                 Q = H(1) H(2) . . . H(k)

as returned by tzrzf. ✨

Q is of order M if "side" = 0 and of order N if "side" = 1. no cap.

📋 Arguments 📋

side — = 0: apply Q or Q' from the Left; = 1: apply Q or Q' from the Right. 👑
trans — = 0: No transpose, apply Q; = 1: Transpose, apply Q'. 💅
A — The i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by tzrzf in the last k rows of its array argument A. A is modified by the routine but restored on exit. 🎀
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by tzrzf.
C — On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. 🦋
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. 😭

$a = random (float, 50, 100);
$a = transpose($a);
$info = null;
$tau = zeroes(float, 50);
tzrzf($a, $tau, $info);
$c = random(70,50);
# $c will contain the result
$c->reshape(70,100);
$c = transpose($c);
ormrz($a, $tau, $c, $info);

ormrz ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 💖

👑 gehrd 👑

Signature: ([io,phys]A(n,n); int [phys]ilo();int [phys]ihi();[o,phys]tau(k); int [o,phys]info())

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

Q' * A * Q = H .

iconic behaviour honestly. ✨

🔍 Further Details 🔍

The matrix Q is represented as a product of (ihi-ilo) elementary reflectors

        Q = H(ilo) H(ilo+1) . . . H(ihi-1).

Each H(i) has the form

        H(i) = I - tau * v * v'
        where tau is a real scalar, and v is a real vector with
        v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
        exit in A(i+2:ihi,i), and tau in tau(i).

The contents of A are illustrated by the following example, with n = 7, ilo = 2 and ihi = 6:

        on entry,                        on exit,

        ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
        (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
        (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
        (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
        (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
        (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
        (                         a )    (                          a )

        where a denotes an element of the original matrix A, h denotes a
        modified element of the upper Hessenberg matrix H, and vi denotes an
        element of the vector defining H(i).

📋 Arguments 📋

ilo, ihi — It is assumed that A is already upper triangular in rows and columns 1:ilo-1 and ihi+1:N. ilo and ihi are normally set by a previous call to gebal; otherwise they should be set to 1 and N respectively. See Further Details. 1 <= ilo <= ihi <= N, if N > 0; ilo=1 and ihi=0, if N=0. 🌸
A — On entry, the N-by-N general matrix to be reduced. On exit, the upper triangle and the first subdiagonal of A are overwritten with the upper Hessenberg matrix H, and the elements below the first subdiagonal, with the array tau, represent the orthogonal matrix Q as a product of elementary reflectors. See Further Details. 💜
tau — The scalar factors of the elementary reflectors (see Further Details). Elements 1:ilo-1 and ihi:N-1 of tau are set to zero. (dimension (N-1)) 🎀
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. 😭

$a = random (50, 50);
$info = null;
$tau = zeroes(50);
gehrd($a, 1, 50, $tau, $info);

gehrd ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🧚 orghr 🧚

Signature: ([io,phys]A(n,n); int [phys]ilo();int [phys]ihi();[phys]tau(k); int [o,phys]info())

generates a real orthogonal matrix Q which is defined as the product of ihi-ilo elementary reflectors of order N, as returned by "gehrd":

        Q = H(ilo) H(ilo+1) . . . H(ihi-1).

📋 Arguments 📋

ilo, ihi — ilo and ihi must have the same values as in the previous call of gehrd. Q is equal to the unit matrix except in the submatrix Q(ilo+1:ihi,ilo+1:ihi). 1 <= ilo <= ihi <= N, if N > 0; ilo=1 and ihi=0, if N=0. 🌸
A — On entry, the vectors which define the elementary reflectors, as returned by gehrd. On exit, the N-by-N orthogonal matrix Q. It’s giving orthogonal realness. ✨
tau — tau(i) must contain the scalar factor of the elementary reflector H(i), as returned by gehrd. (dimension (N-1)) 💅
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value. so if it’s negative, something went wrong, bestie 😭

$a = random (50, 50);
$info = null;
$tau = zeroes(50);
gehrd($a, 1, 50, $tau, $info);
orghr($a, 1, 50, $tau, $info);

orghr literally ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. iykyk 🫶

`hseqr` 👑

okay so hseqr computes the eigenvalues of a real upper Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur form), and Z is the orthogonal matrix of Schur vectors. main character energy 💖

Optionally Z may be postmultiplied into an input orthogonal matrix Q, so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. no cap.

Arguments 💅

job — = 0: compute eigenvalues only; = 1: compute eigenvalues and the Schur form T. simple as that ✨
compz — = 0: no Schur vectors are computed; = 1: Z is initialized to the unit matrix and the matrix Z of Schur vectors of H is returned; = 2: Z must contain an orthogonal matrix Q on entry, and the product Q*Z is returned. choose your fighter 💅
ilo, ihi — It is assumed that H is already upper triangular in rows and columns 1:ilo-1 and ihi+1:N. ilo and ihi are normally set by a previous call to gebal, and then passed to gehrd when the matrix output by gebal is reduced to Hessenberg form. Otherwise ilo and ihi should be set to 1 and N respectively. 1 <= ilo <= ihi <= N, if N > 0; ilo=1 and ihi=0, if N=0. the boundaries, bestie 🦋
H — On entry, the upper Hessenberg matrix H. On exit, if job = 1, H contains the upper quasi-triangular matrix T from the Schur decomposition (the Schur form); 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If job = 0, the contents of H are unspecified on exit. unspecified means don’t trust it, periodt 💀
wr, wi — The real and imaginary parts, respectively, of the computed eigenvalues. If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of wr and wi, say the i-th and (i+1)th, with wi(i) > 0 and wi(i+1) < 0. If job = 1, the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with wr(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, wi(i) = sqrt(H(i+1,i)*H(i,i+1)) and wi(i+1) = -wi(i). math girlies rise 🧚
Z — If compz = 0: Z is not referenced. If compz = 1: on entry, Z need not be set, and on exit, Z contains the orthogonal matrix Z of the Schur vectors of H. If compz = 2: on entry Z must contain an N-by-N matrix Q, which is assumed to be equal to the unit matrix except for the submatrix Z(ilo:ihi,ilo:ihi); on exit Z contains Q*Z. Normally Q is the orthogonal matrix generated by orghr after the call to gehrd which formed the Hessenberg matrix H. the lore 📖
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value; > 0: if info = i, hseqr failed to compute all of the eigenvalues in a total of 30*(ihi-ilo+1) iterations; elements 1:ilo-1 and i+1:n of wr and wi contain those eigenvalues which have been successfully computed. read the fine print 👑

$a = random (50, 50);
$info = null;
$tau = zeroes(50);
$z= zeroes(1,1);
gehrd($a, 1, 50, $tau, $info);
hseqr($a,0,0,1,50,($wr=null),($wi=null),$z,$info);

hseqr ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. fr fr 🫶

`trevc` 🦋

okay so trevc computes some or all of the right and/or left eigenvectors of a real upper quasi-triangular matrix T. eigenvector era 💅

so the right eigenvector x and the left eigenvector y of T for an eigenvalue w are defined like this:

T*x = w*x,     y'*T = w*y'
where y' denotes the conjugate transpose of the vector y.

If all eigenvectors are requested, the routine may either return the matrices X and/or Y of right or left eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an input orthogonal matrix. If T was obtained from the real-Schur factorization of an original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of right or left eigenvectors of A. it’s giving decomposition realness ✨

T has to be in Schur canonical form (like what hseqr returns), meaning block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has its diagonal elements equal and its off-diagonal elements of opposite sign. Corresponding to each 2-by-2 diagonal block is a complex conjugate pair of eigenvalues and eigenvectors; only one eigenvector of the pair is computed, namely the one corresponding to the eigenvalue with positive imaginary part. only one, bestie 💅

Further Details 🌈

The algorithm used in this program is basically backward (forward) substitution, with scaling to make the code robust against possible overflow. scaling for safety 🛡️

Each eigenvector is normalized so that the element of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|. normalize and serve 💅

Arguments 💖

side — = 0 : compute both right and left eigenvectors; = 1 : compute right eigenvectors only; = 2 : compute left eigenvectors only. pick a side, literally 😭
howmny — = 0: compute all right and/or left eigenvectors; = 1: compute all right and/or left eigenvectors, and backtransform them using the input matrices supplied in VR and/or VL; = 2: compute selected right and/or left eigenvectors, specified by the logical array select. so many options ✨
select — If howmny = 2, select specifies the eigenvectors to be computed. If howmny = 0 or 1, select is not referenced. To select the real eigenvector corresponding to a real eigenvalue w(j), select(j) must be set to TRUE. To select the complex eigenvector corresponding to a complex conjugate pair w(j) and w(j+1), either select(j) or select(j+1) must be set to TRUE; then on exit select(j) is TRUE and select(j+1) is FALSE. TRUE/FALSE drama 💅
T — The upper quasi-triangular matrix T in Schur canonical form. canonical queen 👑
VL — On entry, if side = 2 or 0 and howmny = 1, VL must contain an N-by-N matrix Q (usually the orthogonal matrix Q of Schur vectors returned by hseqr). On exit, if side = 2 or 0, VL contains: if howmny = 0, the matrix Y of left eigenvectors of T; VL has the same quasi-lower triangular form as T'. If T(i,i) is a real eigenvalue, then the i-th column VL(i) of VL is its corresponding eigenvector. If T(i:i+1,i:i+1) is a 2-by-2 block whose eigenvalues are complex-conjugate eigenvalues of T, then VL(i)+sqrt(-1)*VL(i+1) is the complex eigenvector corresponding to the eigenvalue with positive real part. if howmny = 1, the matrix Q*Y; if howmny = 2, the left eigenvectors of T specified by select, stored consecutively in the columns of VL, in the same order as their eigenvalues. A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part, and the second the imaginary part. If side = 1, VL is not referenced. left side represent 💅
VR — On entry, if side = 1 or 0 and howmny = 1, VR must contain an N-by-N matrix Q (usually the orthogonal matrix Q of Schur vectors returned by hseqr). On exit, if side = 1 or 0, VR contains: if howmny = 0, the matrix X of right eigenvectors of T; VR has the same quasi-upper triangular form as T. If T(i,i) is a real eigenvalue, then the i-th column VR(i) of VR is its corresponding eigenvector. If T(i:i+1,i:i+1) is a 2-by-2 block whose eigenvalues are complex-conjugate eigenvalues of T, then VR(i)+sqrt(-1)*VR(i+1) is the complex eigenvector corresponding to the eigenvalue with positive real part. if howmny = 1, the matrix Q*X; if howmny = 2, the right eigenvectors of T specified by select, stored consecutively in the columns of VR, in the same order as their eigenvalues. A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part and the second the imaginary part. If side = 2, VR is not referenced. right side represent 💅
m — The number of columns in the arrays VL and/or VR actually used to store the eigenvectors. If howmny = 0 or 1, m is set to N. Each selected real eigenvector occupies one column and each selected complex eigenvector occupies two columns. counting columns, periodt ✨
info — = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value. check it bestie 😭

$a = random (50, 50);
$info = null;
$tau = zeroes(50);
$z= zeroes(1,1);
gehrd($a, 1, 50, $tau, $info);
hseqr($a,0,0,1,50,($wr=null),($wi=null),$z,$info);

okay so trevc literally does not care about the bad-value flag on the ndarrays you feed it 🙈 — like, it just ignores it. But it will set the bad-value flag on all its output ndarrays if the flag is set on any of the input ndarrays. so if you send in one messy flag, bestie, it comes right back out the other side 💅

tgevc 💜

Signature: ([io]A(n,n); int side();int howmny();[io]B(n,n);int select(q);[o]VL(m,m); [o]VR(p,p);int [o]m(); int [o]info(); [t]work(workn))

So this one computes some or all of the right and/or left generalized eigenvectors of a pair of real upper triangular matrices (A,B). iconic behavior honestly ✨

The right generalized eigenvector x and the left generalized eigenvector y of (A,B) that correspond to a generalized eigenvalue w are defined like this:

(A - wB) * x = 0  and  y**H * (A - wB) = 0
where y**H denotes the conjugate transpose of y.

If an eigenvalue w is determined by zero diagonal elements of both A and B, then a unit vector is returned as the corresponding eigenvector. 🫶

If you ask for all the eigenvectors, the routine may either hand you the matrices X and/or Y of right or left eigenvectors of (A,B), or the products Z*X and/or Q*Y, where Z and Q are input orthogonal matrices. And if (A,B) came from the generalized real-Schur factorization of an original pair of matrices

(A0,B0) = (Q*A*Z**H,Q*B*Z**H),

then Z*X and Q*Y are the matrices of right or left eigenvectors of A. 💖

Heads up though: A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal blocks. Corresponding to each 2-by-2 diagonal block is a complex conjugate pair of eigenvalues and eigenvectors; only one eigenvector of the pair gets computed, namely the one corresponding to the eigenvalue with positive imaginary part. no cap. 👑

Arguments

side — = 0 : compute both right and left eigenvectors; = 1 : compute right eigenvectors only; = 2 : compute left eigenvectors only.
howmny — = 0 : compute all right and/or left eigenvectors; = 1 : compute all right and/or left eigenvectors, and backtransform them using the input matrices supplied in VR and/or VL; = 2 : compute selected right and/or left eigenvectors, specified by the logical array select.
select — If howmny=2, select specifies the eigenvectors to be computed. If howmny=0 or 1, select is not referenced. To select the real eigenvector corresponding to the real eigenvalue w(j), select(j) must be set to TRUE. To select the complex eigenvector corresponding to a complex conjugate pair w(j) and w(j+1), either select(j) or select(j+1) must be set to TRUE. 💅
A — The upper quasi-triangular matrix A.
B — The upper triangular matrix B. If A has a 2-by-2 diagonal block, then the corresponding 2-by-2 block of B must be diagonal with positive elements.
VL — On entry, if side = 2 or 0 and howmny = 1, VL must contain an N-by-N matrix Q (usually the orthogonal matrix Q of left Schur vectors returned by hgqez). On exit, if side = 2 or 0, VL contains:
- if howmny = 0, the matrix Y of left eigenvectors of (A,B);
- if howmny = 1, the matrix Q*Y;
- if howmny = 2, the left eigenvectors of (A,B) specified by select, stored consecutively in the columns of VL, in the same order as their eigenvalues.

If side = 1, VL is not referenced.

A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part, and the second the imaginary part. 🌸

VR — On entry, if side = 1 or 0 and howmny = 1, VR must contain an N-by-N matrix Q (usually the orthogonal matrix Z of right Schur vectors returned by hgeqz). On exit, if side = 1 or 0, VR contains:
- if howmny = 0, the matrix X of right eigenvectors of (A,B);
- if howmny = 1, the matrix Z*X;
- if howmny = 2, the right eigenvectors of (A,B) specified by select, stored consecutively in the columns of VR, in the same order as their eigenvalues.

If side = 2, VR is not referenced.

A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part and the second the imaginary part. 🦋

M — The number of columns in the arrays VL and/or VR actually used to store the eigenvectors. If howmny = 0 or 1, M is set to N. Each selected real eigenvector occupies one column and each selected complex eigenvector occupies two columns.
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. > 0: the 2-by-2 block (info:info+1) does not have a complex eigenvalue. 💀

=for example

$a = random (50, 50);
$info = null;
$tau = zeroes(50);
$z= zeroes(1,1);
gehrd($a, 1, 50, $tau, $info);
hseqr($a,0,0,1,50,($wr=null),($wi=null),$z,$info);

tgevc ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

gebal 🎀

Signature: ([io,phys]A(n,n); int job(); int [o,phys]ilo();int [o,phys]ihi();[o,phys]scale(n); int [o,phys]info())

This one balances a general real matrix A. what does that even mean bestie? okay so it's two steps. first, it permutes A by a similarity transformation to isolate eigenvalues in the first 1 to ilo-1 and last ihi+1 to N elements on the diagonal. and second, it applies a diagonal similarity transformation to rows and columns ilo to ihi to make the rows and columns as close in norm as possible. both steps are optional, so you can pick and choose. 💅

Balancing may reduce the 1-norm of the matrix, and improve the accuracy of the computed eigenvalues and/or eigenvectors. so like… it's giving self-care for your matrix. ✨

Further Details

The permutations consist of row and column interchanges which put the matrix in the form

        ( T1   X   Y  )
P A P = (  0   B   Z  )
        (  0   0   T2 )

where T1 and T2 are upper triangular matrices whose eigenvalues lie along the diagonal. The column indices ilo and ihi mark the starting and ending columns of the submatrix B. Balancing consists of applying a diagonal similarity transformation inv(D) * B * D to make the 1-norms of each row of B and its corresponding column nearly equal. 🌈

The output matrix is

( T1     X*D          Y    )
(  0  inv(D)*B*D  inv(D)*Z ).
(  0      0           T2   )

Information about the permutations P and the diagonal matrix D is returned in the vector scale. 👑

Arguments

job — Specifies the operations to be performed on A:
- = 0: none: simply set ilo = 1, ihi = N, scale(I) = 1.0 for i = 1,...,N;
- = 1: permute only;
- = 2: scale only;
- = 3: both permute and scale.
A — On entry, the input matrix A. On exit, A is overwritten by the balanced matrix. If job = 0, A is not referenced. See Further Details. 🌸
ilo, ihi — ilo and ihi are set to integers such that on exit A(i,j) = 0 if i > j and j = 1,...,ilo-1 or I = ihi+1,...,N. If job = 0 or 2, ilo = 1 and ihi = N.
scale — Details of the permutations and scaling factors applied to A. If P(j) is the index of the row and column interchanged with row and column j and D(j) is the scaling factor applied to row and column j, then
- scale(j) = P(j) for j = 1,...,ilo-1
- = D(j) for j = ilo,...,ihi
- = P(j) for j = ihi+1,...,N.

The order in which the interchanges are made is N to ihi+1, then 1 to ilo-1.

info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. 💜

$a = random (50, 50);
$scale = zeroes(50);
$info = null;
$ilo = null;
$ihi = null;
gebal($a, $ilo, $ihi, $scale, $info);

gebal ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

gebak 🦋

Signature: ([io,phys]A(n,m); int job(); int side();int [phys]ilo();int [phys]ihi();[phys]scale(n); int [o,phys]info())

gebak forms the right or left eigenvectors of a real general matrix by backward transformation on the computed eigenvectors of the balanced matrix output by gebal. so it's like the lil glow-up reversal, bestie. ✨

Arguments

A — On entry, the matrix of right or left eigenvectors to be transformed, as returned by hsein or trevc. On exit, A is overwritten by the transformed eigenvectors.
job — Specifies the type of backward transformation required:
- = 0, do nothing, return immediately;
- = 1, do backward transformation for permutation only;
- = 2, do backward transformation for scaling only;
- = 3, do backward transformations for both permutation and scaling.

job must be the same as the argument job supplied to gebal. no cap, they gotta match. 💅

side — = 0: V contains left eigenvectors. = 1: V contains right eigenvectors;
ilo, ihi — The integers ilo and ihi determined by gebal. 1 <= ilo <= ihi <= N, if N > 0; ilo=1 and ihi=0, if N=0. Here N is the the number of rows of the matrix A.
scale — Details of the permutation and scaling factors, as returned by gebal.
info — = 0: successful exit. < 0: if info = -i, the i-th argument had an illegal value. 😭

$a = random (50, 50);
$scale = zeroes(50);
$info = null;
$ilo = null;
$ihi = null;
gebal($a, $ilo, $ihi, $scale, $info);
# Compute eigenvectors ($ev)
gebak($ev, $ilo, $ihi, $scale, $info);

gebak ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

lange ⭐

Signature: ([phys]A(n,m); int norm(); [o]b(); [t]work(workn))

Computes the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real matrix A. pick your fighter, bestie. 💖

Description

returns the value

okay so the formula for lange finishes out like this, bestie — read it like a little choose-your-own-norm menu:

              lange  = ( max(abs(A(i,j))), norm = 0
                       (
                       ( norm1(A),         norm = 1
                       (
                       ( normI(A),         norm = 2
                       (
                       ( normF(A),         norm = 3

where norm1 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). and pls clock this: max(abs(A(i,j))) is not a matrix norm. it's just vibing there in slot 0. 🙈

✨ Arguments ✨

norm — tells lange which value to hand back, exactly as described up there ⬆️
A — the n by m matrix A 💅

$a = random (float, 100, 100);
$norm  = $a->lange(1);

lange ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🌸 `lansy` 🌸

Signature: (A(n,n); int uplo(); int norm(); [o]b(); [t]work(workn))

Computes the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix A. symmetric queens only. 👑

💖 Description 💖

returns the value

              lansy  = ( max(abs(A(i,j))), norm = 0
                       (
                       ( norm1(A),         norm = 1
                       (
                       ( normI(A),         norm = 2
                       (
                       ( normF(A),         norm = 3

where norm1 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). and again, no cap: max(abs(A(i,j))) is not a matrix norm. it just lives there. 😭

norm — which value gets returned in lansy, as described above ⬆️
uplo — says whether we're referencing the upper or lower triangular part of the symmetric matrix A
- = 0: upper triangular part of A is referenced
- = 1: lower triangular part of A is referenced
A — the symmetric matrix A. if uplo = 0, the leading n by n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. if uplo = 1, the leading n by n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. it literally just doesn't look at the other half, it's not being shady, it's being efficient. ✨

# Assume $a is symmetric
$a = random (float, 100, 100);
$norm  = $a->lansy(1, 1);

lansy ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🦋 `lantr` 🦋

Signature: (A(m,n);int uplo();int norm();int diag();[o]b(); [t]work(workn))

Computes the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix A. we love a shape-shifter. 🧚

💜 Description 💜

returns the value

              lantr  = ( max(abs(A(i,j))), norm = 0
                       (
                       ( norm1(A),         norm = 1
                       (
                       ( normI(A),         norm = 2
                       (
                       ( normF(A),         norm = 3

where norm1 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). and fr fr, max(abs(A(i,j))) is not a matrix norm. it's giving "technically a value" but not a norm. 💀

norm — which value gets returned in lantr, as described above ⬆️
uplo — says whether the matrix A is upper or lower trapezoidal
- = 0: upper triangular part of A is referenced
- = 1: lower triangular part of A is referenced
- note that A is triangular instead of trapezoidal if M = N 👀
diag — says whether or not the matrix A has a unit diagonal
- = 0: non-unit diagonal
- = 1: unit diagonal
A — the trapezoidal matrix A (A is triangular if m = n).
- if uplo = 0, the leading m by n upper trapezoidal part of the array A contains the upper trapezoidal matrix, and the strictly lower triangular part of A is not referenced.
- if uplo = 1, the leading m by n lower trapezoidal part of the array A contains the lower trapezoidal matrix, and the strictly upper triangular part of A is not referenced.
- and this bit matters: when diag = 1, the diagonal elements of A are not referenced and are assumed to be one. so don't stress about what's chillin' on that diagonal, nobody's reading it. 🙈

# Assume $a is upper triangular
$a = random (float, 100, 100);
$norm  = $a->lantr(1, 1, 0);

lantr ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🔥 `gemm` 🔥

Signature: ([phys]A(m,n); int transa(); int transb(); [phys]B(p,q);[phys]alpha(); [phys]beta(); [io,phys]C(r,s))

Performs one of the matrix-matrix operations

        C := alpha*op( A )*op( B ) + beta*C,
        where  op( X ) is one of p( X ) = X   or   op( X ) = X',
        alpha and beta are scalars, and A, B and C are matrices, with op( A )
        an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.

this is THE matrix multiplication moment bestie, main character energy, the whole reason you're here. 🌈

⭐ Parameters ⭐

transa — on entry, transa specifies the form of op( A ) to be used in the matrix multiplication, like so:
- transa = 0, op( A ) = A.
- transa = 1, op( A ) = A'.
transb — on entry, transb specifies the form of op( B ) to be used in the matrix multiplication, like so:
- transb = 0, op( B ) = B.
- transb = 1, op( B ) = B'.
alpha — on entry, alpha specifies the scalar alpha. 💅
A — before entry with transa = 0, the leading m by k part of the array A must contain the matrix A, otherwise the leading k by m part of the array A must contain the matrix A.
B — before entry with transb = 0, the leading k by n part of the array B must contain the matrix B, otherwise the leading n by k part of the array B must contain the matrix B.
beta — on entry, beta specifies the scalar beta. when beta is supplied as zero then C need not be set on input. so if beta's a zero you can just leave C blank, she's not being read.
C — before entry, the leading m by n part of the array C must contain the matrix C, except when beta is zero, in which case C need not be set on entry.
- on exit, the array C is overwritten by the m by n matrix ( alpha*op( A )*op( B ) + beta*C ). so yeah she gets fully rewritten, no take-backs.

$a = random(5,4);
$b = random(5,4);
$alpha = pdl(0.5);
$beta = pdl(0);
$c = zeroes(5,5);
gemm($a, 0, 1,$b, $alpha, $beta, $c);

gemm ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🎀 `mmult` 🎀

Signature: ([phys]A(m,n); [phys]B(p,m); [o,phys]C(p,n))

Blas matrix multiplication based on gemm. ✨

mmult ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🍓 `crossprod` 🍓

Signature: ([phys]A(n,m); [phys]B(p,m); [o,phys]C(p,n))

Blas matrix cross product based on gemm. 💖

crossprod ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

👑 `syrk` 👑

Signature: ([phys]A(m,n); int uplo(); int trans(); [phys]alpha(); [phys]beta(); [io,phys]C(p,p))

Performs one of the symmetric rank k operations

        C := alpha*A*A' + beta*C,

        C := alpha*A'*A + beta*C,

        where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
        and  A  is an  n by k  matrix in the first case and a  k by n  matrix
        in the second case.

💜 Parameters 💜

uplo — on entry, uplo specifies whether the upper or lower triangular part of the array C is to be referenced, like so:
- uplo = 0 Only the upper triangular part of C is to be referenced.
- uplo = 1 Only the lower triangular part of C is to be referenced.
- unchanged on exit. ✨
trans — on entry, trans specifies the operation to be performed, like so:
- trans = 0 C := alpha*A*A' + beta*C.
- trans = 1 C := alpha*A'*A + beta*C.
alpha — on entry, alpha specifies the scalar alpha. unchanged on exit. 💅
A — before entry with trans = 0, the leading n by k part of the array A must contain the matrix A, otherwise the leading k by n part of the array A must contain the matrix A.
beta — on entry, beta specifies the scalar beta.
C — before entry with uplo = 0, the leading n by n upper triangular part of the array C must contain the upper triangular part of the symmetric matrix and the strictly lower triangular part of C is not referenced. on exit, the upper triangular part of the array C is overwritten by the upper triangular part of the updated matrix.
- before entry with uplo = 1, the leading n by n lower triangular part of the array C must contain the lower triangular part of the symmetric matrix and the strictly upper triangular part of C is not referenced. on exit, the lower triangular part of the array C is overwritten by the lower triangular part of the updated matrix.

tl;dr: you pick a half, you feed her that half, she updates that half, the other half was never invited to the party. 🙈

$a = random(5,4);
$b = zeroes(5,5);
$alpha = 1;
$beta = 0;
syrk ($a, 1,0,$alpha, $beta , $b);

syrk ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

🧚 `dot` 🧚

Signature: ([phys]a(n);[phys]b(n);[o,phys]c())

Dot product of two vectors using Blas. ⭐

$a = random(5);
$b = random(5);
$c = dot($a, $b)

dot ignores the bad-value flag of the input ndarrays. it will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. 🫶

axpy

Signature:

Signature: ([phys]a(n);[phys] alpha();[io,phys]b(m))

Okay bestie, axpy is giving linear combination energy: it does ax + b using Blas, and the result lands in b. So like, a gets scaled by alpha and added to b. Periodt. ✨

Example:

$a = random(5);
$b = random(5);
axpy($a, 12, $b)

And about bad values: axpy ignores the bad-value flag of the input ndarrays. But if the flag is set for any input ndarray, it will set the bad-value flag of all output ndarrays. So like, bad vibes in, bad vibes out? 💀

nrm2

Signature:

Signature: ([phys]a(n);[o]b())

nrm2 gives you the Euclidean norm of a vector using Blas. So like the length of the vector, bestie. 📏✨

Example:

$a = random(5);
$norm2 = norm2($a)

nrm2 ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. So it's not checking, but it will propagate. 💅

asum

Signature:

Signature: ([phys]a(n);[o]b())

asum is the sum of absolute values of a vector using Blas. So it's like adding up all the absolute values, no negative energy allowed. 💅

Example:

$a = random(5);
$absum = asum($a)

asum ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. No cap.

scal

Signature:

Signature: ([io,phys]a(n);scale())

scal scales a vector by a constant using Blas. So you multiply every element by scale. It's giving glow-up. ✨

Example:

$a = random(5);
$a->scal(0.5)

scal ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

rot

Signature:

Signature: ([io,phys]a(n);[phys]c(); [phys]s();[io,phys]b(n))

rot applies plane rotation using Blas. So it rotates vectors a and b using c and s (cosine and sine vibes). 🔄

Example:

$a = random(5);
$b = random(5);
rot($a, 0.5, 0.7, $b)

rot ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

rotg

Signature:

Signature: ([io,phys]a();[io,phys]b();[o,phys]c(); [o,phys]s())

rotg generates plane rotation using Blas. So it takes a and b and spits out c and s. The math is mathing. ✨

Example:

$a = sequence(4);
rotg($a(0), $a(1),$a(2),$a(3))

rotg ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

lasrt

Signature:

Signature: ([io,phys]d(n); int id();int [o,phys]info())

lasrt sorts the numbers in d in increasing order (if id = 0) or in decreasing order (if id = 1). So you can sort ascending or descending, bestie. 👑

It uses Quick Sort, reverting to Insertion sort on arrays of size <= 20. Dimension of stack limits N to about 232. So like, it's quick but also handles small arrays with insertion sort. And the stack dimension limits N to about 232. No cap.

✨ Arguments ✨

id: = 0: sort d in increasing order; = 1: sort d in decreasing order.
d: On entry, the array to be sorted. On exit, d has been sorted into increasing order (d(1) <= ... <= d(N)) or into decreasing order (d(1) >= ... >= d(N)), depending on id.
info: = 0: successful exit; < 0: if info = -i, the i-th argument had an illegal value.

Example:

$a = random(5);
lasrt ($a, 0, ($info = null));

lasrt ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

lacpy

Signature:

Signature: ([phys]A(m,n); int uplo(); [o,phys]B(p,n))

lacpy copies all or part of a two-dimensional matrix A to another matrix B. So like, copy-paste but with options. 📋

✨ Arguments ✨

uplo: Specifies the part of the matrix A to be copied to B.
- = 0: Upper triangular part
- = 1: Lower triangular part
- Otherwise: All of the matrix A
A: The m by n matrix A. If uplo = 0, only the upper triangle or trapezoid is accessed; if uplo = 1, only the lower triangle or trapezoid is accessed.
B: On exit, B = A in the locations specified by uplo.

Example:

$a = random(5,5);
$b = zeroes($a);
lacpy ($a, 0, $b);

lacpy ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

laswp

Signature:

Signature: ([io,phys]A(m,n);int [phys]k1();int [phys] k2(); int [phys]ipiv(p);int [phys]inc())

laswp performs a series of row interchanges on the matrix A. One row interchange is initiated for each of rows k1 through k2 of A. Doesn't use PDL indices (start = 1). So it's like swapping rows, but note the indexing starts at 1, not 0. PDL girlies, watch out. 🙈

✨ Arguments ✨

A: On entry, the matrix of column dimension N to which the row interchanges will be applied. On exit, the permuted matrix.
k1: The first element of ipiv for which a row interchange will be done.
k2: The last element of ipiv for which a row interchange will be done.
ipiv: The vector of pivot indices. Only the elements in positions k1 through k2 of ipiv are accessed. ipiv(k) = l implies rows k and l are to be interchanged.
inc: The increment between successive values of ipiv. If ipiv is negative, the pivots are applied in reverse order.

Example:

$a = random(5,5);
# reverse row (col for PDL)
$b = pdl([5,4,3,2,1]);
$a->laswp(1,2,$b,1);

laswp ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

lamch

Signature:

Signature: (cmach(); [o]precision())

lamch determines precision machine parameters. Works inplace. So it's giving you the deets on your machine's floating-point precision. 💻

✨ Arguments ✨

cmach: Specifies the value to be returned by lamch:
- = 0 LAMCH := eps
- = 1 LAMCH := sfmin
- = 2 LAMCH := base
- = 3 LAMCH := eps*base
- = 4 LAMCH := t
- = 5 LAMCH := rnd
- = 6 LAMCH := emin
- = 7 LAMCH := rmin
- = 8 LAMCH := emax
- = 9 LAMCH := rmax

where

eps = relative machine precision
sfmin = safe minimum, such that 1/sfmin does not overflow
base = base of the machine
prec = eps*base
t = number of (base) digits in the mantissa
rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
emin = minimum exponent before (gradual) underflow
rmin = underflow threshold - base**(emin-1)
emax = largest exponent before overflow
rmax = overflow threshold - (base**emax)*(1-eps)

Example:

$a = lamch (0);
print "EPS is $a for double\n";

lamch ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

labad

Signature:

Signature: ([io,phys]small(); [io,phys]large())

labad takes the values computed by lamch for underflow and overflow, and returns the square root of each of these values if the log of large is sufficiently large. This subroutine is intended to identify machines with a large exponent range, such as the Crays, and redefine the underflow and overflow limits to be the square roots of the values computed by lamch. This subroutine is needed because lamch does not compensate for poor arithmetic in the upper half of the exponent range, as is found on a Cray. So basically, it's fixing up the limits for machines with huge exponent ranges. 🌈

✨ Arguments ✨

small: On entry, the underflow threshold as computed by lamch. On exit, if LOG10(large) is sufficiently large, the square root of small, otherwise unchanged.
large: On entry, the overflow threshold as computed by lamch. On exit, if LOG10(large) is sufficiently large, the square root of large, otherwise unchanged.

Example:

$underflow = lamch(7);
$overflow = lamch(9);
labad ($underflow, $overflow);

labad ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

tricpy

Signature:

Signature: (A(m,n);[o] C(m,n); int uplo)

tricpy(PDL(A), int(uplo), PDL(C))

Example:

use PDL::LinearAlgebra;

$c = $a->tricpy($uplo); # explicit uplo
$c = $a->tricpy;                        # default upper
or
use PDL::LinearAlgebra::Real;

tricpy($a, $uplo, $c);  # modify c

tricpy copies triangular part to another matrix. If uplo == 0 copy upper triangular part. So it's like copying only the triangle you want. 🔺

tricpy does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. So if bad values are present, it's not handling them, but it will propagate the flag. 💀

cplx_eigen

Signature:

Signature: (eigreval(n);eigimval(n); eigvec(n,p);int fortran();[o]cplx_val(c=2,n);[o]cplx_vec(c=2,n,p))

cplx_eigen outputs complex eigen-values/vectors from eigen-values/vectors as computed by geev or geevx. 'fortran' means fortran storage type. So it's converting real eigen stuff into complex eigen stuff, with a Fortran storage option. 🧬

cplx_eigen does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

augment

Signature:

Signature: (x(n); y(p);[o]out(q))

augment combines two ndarrays into a single ndarray. This routine does backward and forward dataflow automatically. So it's like merging two arrays into one, and the dataflow goes both ways. ✨

augment does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

mstack

Signature:

Signature: (x(n,m);y(n,p);[o]out(n,q))

mstack combines two ndarrays into a single ndarray. This routine does backward and forward dataflow automatically. So like stacking them side by side? It's combining, bestie. 🧚

mstack does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

charpol

Signature:

Signature: (A(n,n);[o]Y(n,n);[o]out(p))

charpol computes the adjoint matrix and the characteristic polynomial. Math queen. 👑

charpol does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.

👑 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.

✨ That's the whole vibe, bestie! This page is from perl v5.38.2, dated 2024-03-31, for Real(3pm). Stay slaying! 💅💖

📦 where this page comes from

ubuntu noble: libpdl-linearalgebra-perl 0.38-1build2 · Artistic-2
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-10-14: a language model's rewrite of the original, under its terms

Real(3pm)

💖 Real(3pm) — PDL::LinearAlgebra::Real 💖

📛 NAME

🎀 SYNOPSIS

🌸 DESCRIPTION

🔥 FUNCTIONS

🦋 gtsv

📝 Arguments

🧚 gesvd

🎯 Arguments

💅 gesdd

🎯 Arguments

🌈 ggsvd

Arguments ✨

geev 🔥

Arguments ✨

geevx ✨

Arguments ✨

ggev 🦋

Arguments 💫

ggevx 🎀

Further Details 🌈

Arguments 💫

gees

Arguments

🌈 geesx 🌈

🦋 gges 🦋

🦋 ggesx 🦋

Further details

Arguments

🌸 syev 🌸

Arguments

🦋 syevd

💖 Arguments 💖

🌸 syevx

🎀 Arguments 🎀

🧚 syevr

⭐ Arguments ⭐

🦋 sygv 🦋

Arguments

🌸 sygvd 🌸

Arguments

💜 sygvx 💜

Arguments

gesv 👑

Arguments

gesvx 🦋

Arguments

🦋 sysv 🦋

💖 Arguments 💖

🧚 sysvx 🧚

✨ Arguments ✨

🌸 posv 🌸

🎀 Arguments 🎀

🦋 posvx 🦋

✨ Arguments ✨

gels

Arguments

gelsy

Arguments

gelss

Arguments

gelsd

Arguments

gglse

Arguments

ggglm

Arguments

getrf

Arguments

getf2

📋 Arguments 📋

🦋 sytrf 🦋

📋 Arguments 📋

🔍 Further Details 🔍

🌸 sytf2 🌸

📋 Arguments 📋

👑 potrf 👑

📋 Arguments 📋

🧚 potf2 🧚

`gees`

🦋 `sygv` 🦋

🌸 `sygvd` 🌸

💜 `sygvx` 💜

`gesv` 👑

`gesvx` 🦋

🌸 `gelqf` 🌸

🌈 `orglq` 🌈

🎀 `ormlq` 🎀

⭐ `geqlf` ⭐

💖 `orgql` 💖

🔥 `ormql` 🔥

🦋 `gerqf` 🦋

`hseqr` 👑

`trevc` 🦋