complex16OTHERauxiliary(3)

📦

version 3.11.0 · 2 versions

3.11.0 · debian bookworm (latest)
3.10.0 · ubuntu jammy

, in debian bookworm

✨ complex16OTHERauxiliary(3) ✨

LAPACK · complex16 · part 1 of 26, bestie, we're going on a JOURNEY 🦋

🎀 NAME 🎀

complex16OTHERauxiliary — it's giving complex16 💅

okay so. hi. welcome to the group chat that is the LAPACK manual 💖

This page is called complex16OTHERauxiliary, and the one thing it's telling us up top is that it belongs to the complex16 family. That's it. That's the whole NAME section. Short, sweet, main character energy, no notes 🫶

But here's the thing — this page is 26 parts long, so do not scroll away, because we are literally just getting started. The NAME is the lil intro slide before the actual lore drops. Think of it as the "hi besties" before the 47-minute video essay. 👑

Stay with me. More is coming. 🌸

✨ SYNOPSIS ✨

🧚♀️ Functions 🧚♀️

okay bestie, here's the whole squad of routines. signatures are exact, no cap. descriptions are rewritten with love, but every fact stays. 💖

subroutine clag2z (M, N, SA, LDSA, A, LDA, INFO) — CLAG2Z converts a complex single precision matrix to a complex double precision matrix. glow-up from single to double, periodt ✨
double precision function dzsum1 (N, CX, INCX) — DZSUM1 forms the 1-norm of the complex vector using the true absolute value. no fake math here 💅
integer function ilazlc (M, N, A, LDA) — ILAZLC scans a matrix for its last non-zero column. detective energy 🔍
integer function ilazlr (M, N, A, LDA) — ILAZLR scans a matrix for its last non-zero row. row hunt, but make it Fortran 🛶
subroutine zdrscl (N, SA, SX, INCX) — ZDRSCL multiplies a vector by the reciprocal of a real scalar. reciprocal realness 💖
subroutine zlabrd (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY) — ZLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form. bidiagonal? more like bi-iconic ✨
subroutine zlacgv (N, X, INCX) — ZLACGV conjugates a complex vector. flip the imaginary vibe 💫
subroutine zlacn2 (N, V, X, EST, KASE, ISAVE) — ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. reverse communication? main character energy 👑
subroutine zlacon (N, V, X, EST, KASE) — ZLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. same vibe, different signature 💅
subroutine zlacp2 (UPLO, M, N, A, LDA, B, LDB) — ZLACP2 copies all or part of a real two-dimensional array to a complex array. real-to-complex pipeline 💖
subroutine zlacpy (UPLO, M, N, A, LDA, B, LDB) — ZLACPY copies all or part of one two-dimensional array to another. copy-paste but make it Fortran ✨
subroutine zlacrm (M, N, A, LDA, B, LDB, C, LDC, RWORK) — ZLACRM multiplies a complex matrix by a square real matrix. math but make it complex 💅
subroutine zlacrt (N, CX, INCX, CY, INCY, C, S) — ZLACRT performs a linear transformation of a pair of complex vectors. transformation station 🚂
complex *16 function zladiv (X, Y) — ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. safe division, bestie 🛡️
subroutine zlaein (RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO) — ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. eigenvector era 🔥
subroutine zlaev2 (A, B, C, RT1, RT2, CS1, SN1) — ZLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. tiny matrix, big energy ✨
subroutine zlag2c (M, N, A, LDA, SA, LDSA, INFO) — ZLAG2C converts a complex double precision matrix to a complex single precision matrix. precision shift, bestie 💅
subroutine zlags2 (UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ) — ZLAGS2. original description is literally just "ZLAGS2", and honestly? iconic minimalism. no cap. 💀
subroutine zlagtm (TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB) — ZLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix, B and C are rectangular matrices, and α and β are scalars, which may be 0, 1, or -1. tridiagonal? more like tri-fabulous 🌈
subroutine zlahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO) — ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. QR algorithm? more like QR queen 👑
subroutine zlahr2 (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY) — ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. zero them out, bestie 💅
subroutine zlaic1 (JOB, J, X, SEST, W, GAMMA, SESTPR, S, C) — ZLAIC1 applies one step of incremental condition estimation. one step at a time, queen ✨
double precision function zlangt (NORM, N, DL, D, DU) — ZLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. norm realness 💖
double precision function zlanhb (NORM, UPLO, N, K, AB, LDAB, WORK) — ZLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix. band matrix, but make it Hermitian 💅
double precision function zlanhp (NORM, UPLO, N, AP, WORK) — ZLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form. packed and snatched ✨
double precision function zlanhs (NORM, N, A, LDA, WORK) — ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix. Hessenberg hive 👑
double precision function zlanht (NORM, N, D, E) — ZLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix. tridiagonal but make it Hermitian 💖
double precision function zlansb (NORM, UPLO, N, K, AB, LDAB, WORK) — ZLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix. symmetric slay ✨
double precision function zlansp (NORM, UPLO, N, AP, WORK) — ZLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form. packed and perfect 💅
double precision function zlantb (NORM, UPLO, DIAG, N, K, AB, LDAB, WORK) — ZLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular band matrix. triangular but make it banded 🎀
double precision function zlantp (NORM, UPLO, DIAG, N, AP, WORK) — ZLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. packed triangular realness ✨
double precision function zlantr (NORM, UPLO, DIAG, M, N, A, LDA, WORK) — ZLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix. trapezoidal? more like trap-queen 👑
subroutine zlapll (N, X, INCX, Y, INCY, SSMIN) — ZLAPLL measures the linear dependence of two vectors. codependency check 💅
subroutine zlapmr (FORWRD, M, N, X, LDX, K) — ZLAPMR rearranges rows of a matrix as specified by a permutation vector. row shuffle realness ✨
subroutine zlapmt (FORWRD, M, N, X, LDX, K) — ZLAPMT performs a forward or backward permutation of the columns of a matrix. column chaos, controlled 💖
subroutine zlaqhb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED) — ZLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ. scale it, bestie 💅
subroutine zlaqhp (UPLO, N, AP, S, SCOND, AMAX, EQUED) — ZLAQHP scales a Hermitian matrix stored in packed form. packed and scaled ✨
subroutine zlaqp2 (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK) — ZLAQP2 computes a QR factorization with column pivoting of the matrix block. QR with a pivot? periodt 👑
subroutine zlaqps (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF) — ZLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. level 3? highkey iconic 💅
subroutine zlaqr0 (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO) — ZLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. eigenvalues and optional Schur tea ☕
subroutine zlaqr1 (N, H, LDH, S1, S2, V) — ZLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. tiny matrix, big shift energy ✨
subroutine zlaqr2 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK) — ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). aggressive? we love a go-getter 💅
subroutine zlaqr3 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK) — ZLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). twin behavior 💖
subroutine zlaqr4 (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO) — ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. same energy, different routine ✨
subroutine zlaqr5 (WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH) — ZLAQR5 performs a single small-bulge multi-shift QR sweep. sweep it up 💅
subroutine zlaqsb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED) — ZLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. scale and slay ✨
subroutine zlaqsp (UPLO, N, AP, S, SCOND, AMAX, EQUED) — ZLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. packed and scaled, no cap 💖
subroutine zlar1v (N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK) — ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. inverse submatrix? mathematically unhinged, we love it 💅
subroutine zlar2v (N, X, Y, Z, INCX, C, S, INCC) — ZLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. rotations from both sides? balanced, as all things should be ⚖️
subroutine zlarcm (M, N, A, LDA, B, LDB, C, LDC, RWORK) — ZLARCM copies all or part of a real two-dimensional array to a complex array. real-to-complex pipeline 💅
subroutine zlarf (SIDE, M, N, V, INCV, TAU, C, LDC, WORK) — ZLARF applies an elementary reflector to a general rectangular matrix. reflector realness ✨
subroutine zlarfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK) — ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix. block reflector, but make it transpose 💅
subroutine zlarfb_gett (IDENT, M, N, K, T, LDT, A, LDA, B, LDB, WORK, LDWORK) — ZLARFB_GETT. original description is literally just "ZLARFB_GETT", and honestly that's a mood. 💀
subroutine zlarfg (N, ALPHA, X, INCX, TAU) — ZLARFG generates an elementary reflector (Householder matrix). Householder? more like house-iconic 👑
subroutine zlarfgp (N, ALPHA, X, INCX, TAU) — ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. beta but make it non-negative ✨
subroutine zlarft (DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) — ZLARFT forms the triangular factor T of a block reflector H = I - vtvH. triangular factor? more like tri-factor 💅
subroutine zlarfx (SIDE, M, N, V, TAU, C, LDC, WORK) — ZLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10. loop unrolling? performance queen ✨
subroutine zlarfy (UPLO, N, V, INCV, TAU, C, LDC, WORK) — ZLARFY. original description is literally "ZLARFY", and we respect the mystery. 🕵️♀️
subroutine zlargv (N, X, INCX, Y, INCY, C, INCC) — ZLARGV generates a vector of plane rotations with real cosines and complex sines. rotation generation station 🚂
subroutine zlarnv (IDIST, ISEED, N, X) — ZLARNV returns a vector of random numbers from a uniform or normal distribution. random but make it math 💅
subroutine zlarrv (N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, — (description cut off in this part, more coming bestie 💖)

...and we're so back, bestie 💖 the synopsis was literally mid-sentence, so let's just pick it right back up like nothing happened:

RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
WORK, IWORK, INFO)

ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. that's her whole job, and she ate 🔥

okay and here's the rest of the lineup, one icon after another:

subroutine zlartv (N, X, INCX, Y, INCY, C, S, INCC) — ZLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a pair of vectors.
subroutine zlascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO) — ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlaset (UPLO, M, N, ALPHA, BETA, A, LDA) — ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine zlasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA) — ZLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine zlaswp (N, A, LDA, K1, K2, IPIV, INCX) — ZLASWP performs a series of row interchanges on a general rectangular matrix.
subroutine zlat2c (UPLO, N, A, LDA, SA, LDSA, INFO) — ZLAT2C converts a double complex triangular matrix to a complex triangular matrix.
subroutine zlatbs (UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO) — ZLATBS solves a triangular banded system of equations.
subroutine zlatdf (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV) — ZLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate.
subroutine zlatps (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO) — ZLATPS solves a triangular system of equations with the matrix held in packed storage.
subroutine zlatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW) — ZLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an unitary similarity transformation.
subroutine zlatrs (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO) — ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
subroutine zlauu2 (UPLO, N, A, LDA, INFO) — ZLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm).
subroutine zlauum (UPLO, N, A, LDA, INFO) — ZLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm).
subroutine zrot (N, CX, INCX, CY, INCY, C, S) — ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
subroutine zspmv (UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY) — ZSPMV computes a matrix-vector product for complex vectors using a complex symmetric packed matrix
subroutine zspr (UPLO, N, ALPHA, X, INCX, AP) — ZSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix.
subroutine ztprfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK) — ZTPRFB applies a complex 'triangular-pentagonal' block reflector to a complex matrix, which is composed of two blocks.

🌈 Detailed Description 🌈

This is the group of complex16 other auxiliary routines 🫶 literally the catch-all bestie group for complex16, and we love her for it.

💅 Function Documentation 💅

🫧 `subroutine clag2z (integer M, integer N, complex, dimension( ldsa, * ) SA, integer LDSA, complex16, dimension( lda, ) A, integer LDA, integer INFO)`

CLAG2Z converts a complex single precision matrix to a complex double precision matrix. it's giving glow-up but for your precision 💅

✨ Purpose

CLAG2Z takes a COMPLEX matrix, SA, and turns it into a COMPLEX*16 matrix, A.

ngl here's a lil detail: it is possible to overflow while converting from double to single, but it is not possible to overflow when converting from single to double. so that direction? safe 🫶

also this is an auxiliary routine, so there is no argument checking. she's not double-checking your work, she trusts you.

🎀 Parameters

M — M is INTEGER. The number of lines of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
SA — SA is COMPLEX array, dimension (LDSA,N). On entry, the M-by-N coefficient matrix SA.
LDSA — LDSA is INTEGER. The leading dimension of the array SA. LDSA >= max(1,M).
A — A is COMPLEX*16 array, dimension (LDA,N). On exit, the M-by-N coefficient matrix A.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
INFO — INFO is INTEGER. = 0: successful exit

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🦋 `double precision function dzsum1 (integer N, complex16, dimension( ) CX, integer INCX)`

DZSUM1 forms the 1-norm of the complex vector using the true absolute value. main character energy fr fr ✨

✨ Purpose

DZSUM1 adds up the absolute values of a complex vector and hands you back a double precision result.

it's based on DZASUM from the Level 1 BLAS — the one change is that it uses the 'genuine' absolute value instead. upgrade behavior 💅

🎀 Parameters

N — N is INTEGER. The number of elements in the vector CX.
CX — CX is COMPLEX*16 array, dimension (N). The vector whose elements will be summed.
INCX — INCX is INTEGER. The spacing between successive values of CX. INCX > 0.

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

💜 Contributors:

Nick Higham for use with ZLACON.

🧚 `integer function ilazlc (integer M, integer N, complex16, dimension( lda, ) A, integer LDA)`

ILAZLC scans a matrix for its last non-zero column. she's the finder bestie 🔍

✨ Purpose

ILAZLC scans A for its last non-zero column.

🎀 Parameters

M — M is INTEGER. The number of rows of the matrix A.
N — N is INTEGER. The number of columns of the matrix A.
A — A is COMPLEX*16 array, dimension (LDA,N). The m by n matrix A.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🌸 `integer function ilazlr (integer M, integer N, complex16, dimension( lda, ) A, integer LDA)`

ILAZLR scans a matrix for its last non-zero row. same vibe, different axis 🌈

✨ Purpose

ILAZLR scans A for its last non-zero row.

🎀 Parameters

M — M is INTEGER. The number of rows of the matrix A.
N — N is INTEGER. The number of columns of the matrix A.
A — A is COMPLEX*16 array, dimension (LDA,N). The m by n matrix A.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🍓 `subroutine zdrscl (integer N, double precision SA, complex16, dimension( ) SX, integer INCX)`

ZDRSCL multiplies a vector by the reciprocal of a real scalar. division, but make it safe 💅

✨ Purpose

ZDRSCL multiplies an n-element complex vector x by the real scalar 1/a. and the way she does it? without overflow or underflow — as long as the final result x/a doesn't overflow or underflow. periodt.

🎀 Parameters

N — N is INTEGER. The number of components of the vector x.
SA — SA is DOUBLE PRECISION. The scalar a which is used to divide each component of x. SA must be >= 0, or the subroutine will divide by zero. 🚨 bestie watch out
SX — SX is COMPLEX*16 array, dimension (1+(N-1)*abs(INCX)). The n-element vector x.
INCX — INCX is INTEGER. The increment between successive values of the vector SX. > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

💀 `subroutine zlabrd (integer M, integer N, integer NB, complex16, dimension( lda, ) A, integer LDA, double precision, dimension( * ) D, double precision, dimension( * ) E, complex16, dimension( ) TAUQ, complex16, dimension( ) TAUP, complex16, dimension( ldx, ) X, integer LDX, complex16, dimension( ldy, ) Y, integer LDY)`

ZLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form. iconic reduction behavior 🔥

✨ Purpose

ZLABRD reduces the first NB rows and columns of a complex general m by n matrix A to upper or lower real bidiagonal form by a unitary transformation Q**H * A * P, and returns the matrices X and Y which are needed to apply the transformation to the unreduced part of A.

if m >= n, A is reduced to upper bidiagonal form; if m < n, it goes to lower bidiagonal form. choose your fighter 🦋

also heads up: this is an auxiliary routine called by ZGEBRD.

🎀 Parameters

M — M is INTEGER. The number of rows in the matrix A.
N — N is INTEGER. The number of columns in the matrix A.
NB — NB is INTEGER. The number of leading rows and columns of A to be reduced.
A — A is COMPLEX*16 array, dimension (LDA,N). On entry, the m by n general matrix to be reduced. On exit, the first NB rows and columns of the matrix are overwritten; the rest of the array is unchanged. If m >= n, elements on and below the diagonal in the first NB columns, with the array TAUQ, represent the unitary matrix Q as a product of elementary reflectors; and elements above the diagonal in the first NB rows, with the array TAUP, represent the unitary matrix P as a product of elementary reflectors. If m < n, elements below the diagonal in the first NB columns, with the array TAUQ, represent the unitary matrix Q as a product of elementary reflectors, and elements on and above the diagonal in the first NB rows, with the array TAUP, represent the unitary matrix P as a product of elementary reflectors. See Further Details. 👀
LDA —
LDA — INTEGER, bestie. It's the leading dimension of the array A. Rule of thumb: LDA >= max(1,M). Don't fight it, just respect it. 💅
D — DOUBLE PRECISION array, dimension (NB). These are the diagonal elements of the first NB rows and columns of the reduced matrix. So like, D(i) = A(i,i). Iconic little diagonal moment. ✨
E — DOUBLE PRECISION array, dimension (NB). The off-diagonal elements of the first NB rows and columns of the reduced matrix. The spicy ones. 🔥
TAUQ — COMPLEX*16 array, dimension (NB). The scalar factors of the elementary reflectors that represent the unitary matrix Q. Peep Further Details below for the tea. 🫖
TAUP — COMPLEX*16 array, dimension (NB). The scalar factors of the elementary reflectors that represent the unitary matrix P. Again, Further Details has you covered. 🎀
X — COMPLEX*16 array, dimension (LDX,NB). The m-by-nb matrix X you need to update the unreduced part of A. Very much a supporting character with main character energy. 👑
LDX — INTEGER. Leading dimension of the array X. And yes: LDX >= max(1,M). ✨
Y — COMPLEX*16 array, dimension (LDY,NB). The n-by-nb matrix Y required to update the unreduced part of A. The other half of the duo, fr fr. 💖
LDY — INTEGER. Leading dimension of the array Y. Say it with me: LDY >= max(1,N). 🌸

✨ Author ✨

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

💜 Further Details: 💜

Okay so the matrices Q and P are represented as products of elementary reflectors, and here's the whole vibe:

Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)

Each H(i) and G(i) has the form:

H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H

where tauq and taup are complex scalars, and v and u are complex vectors. Scalars and vectors, living their best lives. 🦋

If m >= n, then v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). 🌈

If m < n, then v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). ⭐

The elements of the vectors v and u together form the m-by-nb matrix V and the nb-by-n matrix U**H which are needed, with X and Y, to apply the transformation to the unreduced part of the matrix, using a block update of the form: A := A - V*Y**H - X*U**H. That's the whole plot twist right there. 😭

The contents of A on exit are illustrated by the following examples with nb = 2:

m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):

  (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
  (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
  (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
  (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
  (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
  (  v1  v2  a   a   a  )

where a denotes an element of the original matrix which is unchanged, vi denotes an element of the vector defining H(i), and ui an element of the vector defining G(i). So like, a is the untouched bestie, vi and ui are the glow-ups. 💖

🌸 `subroutine zlacgv (integer N, complex16, dimension( ) X, integer INCX)` 🌸

ZLACGV conjugates a complex vector. That's it. That's the whole job. Short, sweet, slays. 💅

Purpose:

ZLACGV conjugates a complex vector of length N. ✨

Parameters

N — INTEGER. The length of the vector X. And N >= 0, no negativity allowed here. 🙈
X — COMPLEX*16 array, dimension (1+(N-1)*abs(INCX)). On entry, the vector of length N to be conjugated. On exit, X is overwritten with conjg(X). Overwritten! Gone! Transformed! 🦋
INCX — INTEGER. The spacing between successive elements of X. 🎀

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

💜 `subroutine zlacn2 (integer N, complex16, dimension( ) V, complex16, dimension( ) X, double precision EST, integer KASE, integer, dimension( 3 ) ISAVE)` 💜

ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. Reverse communication is giving "we'll call you back" energy and honestly? It works. 😭

Purpose:

ZLACN2 estimates the 1-norm of a square, complex matrix A. Reverse communication is used for evaluating matrix-vector products. 💅

Parameters

N — INTEGER. The order of the matrix. And N >= 1. ✨
V — COMPLEX*16 array, dimension (N). On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned). W stays mysterious. We don't get to see W. W is a private person. 🙈
X — COMPLEX*16 array, dimension (N). On an intermediate return, X should be overwritten by

A * X,   if KASE=1,
A**H * X,  if KASE=2,

where A**H is the conjugate transpose of A, and ZLACN2 must be re-called with all the other parameters unchanged. UNCHANGED. I'm serious. Do not touch the other params. 🔥

EST — DOUBLE PRECISION. On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be unchanged from the previous call to ZLACN2. On exit, EST is an estimate (a lower bound) for norm(A). Lower bound btw — it's an estimate, not the whole truth. 💖
KASE — INTEGER. On the initial call to ZLACN2, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether X should be overwritten by A * X or A**H * X. On the final return from ZLACN2, KASE will again be 0. Full circle moment. 🌈
ISAVE — INTEGER array, dimension (3). ISAVE is used to save variables between calls to ZLACN2. It's the little memory box. The group chat archive, if you will. 🫶

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Further Details:

Originally named CONEST, dated March 16, 1988. She's vintage, she's classic, she's been here since the 80s. 👑

Last modified: April, 1999. ✨

This is a thread safe version of ZLACON, which uses the array ISAVE in place of a SAVE statement, as follows:

ZLACON     ZLACN2
 JUMP     ISAVE(1)
 J        ISAVE(2)
 ITER     ISAVE(3)

Thread safe = she can be in multiple conversations at once and never gets confused. Multitasking queen. 💅

Contributors:

Nick Higham, University of Manchester

References:

N.J. Higham, 'FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. 📚

🎀 `subroutine zlacon (integer N, complex16, dimension( n ) V, complex16, dimension( n ) X, double precision EST, integer KASE)` 🎀

ZLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. Same energy as ZLACN2 but this one's the original recipe, iykyk. 🦋

Purpose:

ZLACON estimates the 1-norm of a square, complex matrix A. Reverse communication is used for evaluating matrix-vector products. 💜

Parameters

N — INTEGER. The order of the matrix. And N >= 1. ✨
V — COMPLEX*16 array, dimension (N). On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned). Still no W. W is a ghost. We respect it. 👻
X — COMPLEX*16 array, dimension (N). On an intermediate return, X should be overwritten by

A * X,   if KASE=1,
A**H * X,  if KASE=2,

where A**H is the conjugate transpose of A, and ZLACON must be re-called with all the other parameters unchanged. Again with the unchanged params. This is not a suggestion, bestie. 🔥

EST — DOUBLE PRECISION. On entry with KASE = 1 or 2 and JUMP = 3, EST should be unchanged from the previous call to ZLACON. On exit, EST is an estimate (a lower bound) for norm(A). 💖
KASE — INTEGER. On the initial call to ZLACON, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether X should be overwritten by A * X or A**H * X. On the final return from ZLACON, KASE will again be 0. Zero to one to two back to zero. A whole character arc. 🌈

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Further Details:

Originally named CONEST, dated March 16, 1988. Last modified: April, 1999. ✨

Contributors:

Nick Higham, University of Manchester

References:

N.J. Higham, 'FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. 📚

💖 `subroutine zlacp2 (character UPLO, integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, complex16, dimension( ldb, ) B, integer LDB)` 💖

ZLACP2 copies all or part of a real two-dimensional array to a complex array. Glow-up pipeline: real in, complex out. 🦋

Purpose:

ZLACP2 copies all or part of a real two-dimensional matrix A to a complex matrix B. ✨

Parameters

UPLO — CHARACTER*1. Specifies the part of the matrix A to be copied to B.

= 'U':      Upper triangular part
= 'L':      Lower triangular part
Otherwise:  All of the matrix A

Pick your triangle, or pick everything. No wrong answers, just different vibes. 💅

M — INTEGER. The number of rows of the matrix A. And M >= 0. ✨
N — INTEGER. The number of columns of the matrix A. And N >= 0. 🌸
A — DOUBLE PRECISION array, dimension (LDA,N). The m by n matrix A. If UPLO = 'U', only the upper trapezium is accessed; if UPLO = 'L', only the lower trapezium is accessed. Trapezium!! Not triangle!! The manpage said trapezium and we're honoring it!! 🔥
LDA — INTEGER. The leading dimension of the array A. And you already know: LDA >= max(1,M). 💅
B — COMPLEX*16 array, dimension (LDB,N). On exit, B = A in the locations specified by UPLO. On EXIT. So B is empty until she's done. Give her a second. 🫶
LDB — INTEGER. The leading dimension of the array B. And LDB >= max(1,M). ⭐

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🌈 `subroutine zlacpy (character UPLO, integer M, integer N, complex16, dimension( lda, ) A, integer LDA, complex16, dimension( ldb, ) B, integer LDB)` 🌈

ZLACPY copies all or part of one two-dimensional array to another. Yes, the name literally is "copy". She's honest about her job. We stan a transparent queen. 👑

Purpose:

Okay so like — the UPLO of it all first, because that's what decides which chunk of A gets yeeted into B:

UPLO — CHARACTER*1. Tells us which part of matrix A gets copied into B. It's giving choose-your-own-adventure:
- = 'U': only the upper triangular part 🌸
- = 'L': only the lower triangular part 💜
- Otherwise: all of the matrix A, no crumbs left behind ✨
M — INTEGER. The number of rows of matrix A. M >= 0, periodt.
N — INTEGER. The number of columns of matrix A. N >= 0, fr fr.
A — COMPLEX*16 array, dimension (LDA,N). It's the m by n matrix A. Heads up bestie: if UPLO = 'U', only the upper trapezium is accessed; if UPLO = 'L', only the lower trapezium is accessed. So like, don't stress about the other half, she's not looking at it 🫶
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). Respect the leading dimension or perish 💀
B — COMPLEX*16 array, dimension (LDB,N). On exit, B = A in the locations specified by UPLO. Slay, she copied 💅
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,M).

Author — Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver, NAG Ltd. (the absolute icons who wrote all this down 🎀)

💖 subroutine zlacrm 💖

subroutine zlacrm (integer M, integer N, complex*16, dimension( lda, * ) A,
    integer LDA, double precision, dimension( ldb, * ) B, integer LDB,
    complex*16, dimension( ldc, * ) C, integer LDC, double precision,
    dimension( * ) RWORK)

ZLACRM multiplies a complex matrix by a square real matrix. ✨

Purpose:

It's a "very simple" matrix-matrix multiplication, and by very simple I mean it says so itself. She does:

C := A * B,

where A is M by N and complex; B is N by N and real; C is M by N and complex. So like, mixed-vibes multiplication, one of them's real, one's complex, and yet they make it work 🦋

Parameters

M — INTEGER. The number of rows of matrix A and of matrix C. M >= 0.
N — INTEGER. The number of columns and rows of matrix B and the number of columns of matrix C. N >= 0. She's pulling double duty, main character energy 👑
A — COMPLEX*16 array, dimension (LDA, N). On entry, A contains the M by N matrix A.
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — DOUBLE PRECISION array, dimension (LDB, N). On entry, B contains the N by N matrix B.
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — COMPLEX*16 array, dimension (LDC, N). On exit, C contains the M by N matrix C.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,N).
RWORK — DOUBLE PRECISION array, dimension (2*M*N). Yes bestie, 2*M*N, she needs the space.

Author — Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver, NAG Ltd.

🌈 subroutine zlacrt 🌈

subroutine zlacrt (integer N, complex*16, dimension( * ) CX, integer INCX,
    complex*16, dimension( * ) CY, integer INCY, complex*16 C, complex*16 S)

ZLACRT performs a linear transformation of a pair of complex vectors. 💅

Purpose:

She performs the operation

   (  c  s )( x )  ==> ( x )
   ( -s  c )( y )      ( y )

where c and s are complex and the vectors x and y are complex. It's giving rotation but make it complex, and honestly? Obsessed.

Parameters

N — INTEGER. The number of elements in the vectors CX and CY.
CX — COMPLEX*16 array, dimension (N). On input, the vector x. On output, CX is overwritten with c*x + s*y. She said new value, who dis 🙈
INCX — INTEGER. The increment between successive values of CX. INCX <> 0. Like literally cannot be zero, do not try it.
CY — COMPLEX*16 array, dimension (N). On input, the vector y. On output, CY is overwritten with -s*x + c*y.
INCY — INTEGER. The increment between successive values of CY. INCY <> 0. Same rule, no zeros allowed ✋
C — COMPLEX*16.
S — COMPLEX*16. C and S define the matrix

   [  C   S  ].
   [ -S   C  ]

Author — Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver, NAG Ltd.

🧚 complex*16 function zladiv 🧚

complex*16 function zladiv (complex*16 X, complex*16 Y)

ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. ✨

Purpose:

ZLADIV := X / Y, where X and Y are complex.  The computation of X / Y
will not overflow on an intermediary step unless the results
overflows.

So like, she's dividing, but carefully, and the only way it overflows is if the actual result overflows. She's not out here creating chaos mid-calculation, she's too classy for that 🎀

Parameters

X — COMPLEX*16.
Y — COMPLEX*16. The complex scalars X and Y. That's it, that's the tea.

Author — Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver, NAG Ltd.

🔥 subroutine zlaein 🔥

subroutine zlaein (logical RIGHTV, logical NOINIT, integer N, complex*16,
    dimension( ldh, * ) H, integer LDH, complex*16 W, complex*16, dimension(
    * ) V, complex*16, dimension( ldb, * ) B, integer LDB, double precision,
    dimension( * ) RWORK, double precision EPS3, double precision SMLNUM,
    integer INFO)

ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. 👑

Purpose:

ZLAEIN uses inverse iteration to find a right or left eigenvector corresponding to the eigenvalue W of a complex upper Hessenberg matrix H. Inverse iteration, bestie. She's iterating backwards into the answer and honestly it works every time.

Parameters

RIGHTV — LOGICAL. Pick your fighter:
- = .TRUE.: compute right eigenvector;
- = .FALSE.: compute left eigenvector.
NOINIT — LOGICAL. Another pick-your-fighter:
- = .TRUE.: no initial vector supplied in V
- = .FALSE.: initial vector supplied in V.
N — INTEGER. The order of the matrix H. N >= 0.
H — COMPLEX*16 array, dimension (LDH,N). The upper Hessenberg matrix H.
LDH — INTEGER. The leading dimension of the array H. LDH >= max(1,N).
W — COMPLEX*16. The eigenvalue of H whose corresponding right or left eigenvector is to be computed.
V — COMPLEX*16 array, dimension (N). On entry, if NOINIT = .FALSE., V must contain a starting vector for inverse iteration; otherwise V need not be set. On exit, V contains the computed eigenvector, normalized so that the component of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|. Not the usual magnitude formula, bestie — it's the |x| + |y| one. Write that down 📝
B — COMPLEX*16 array, dimension (LDB,N).
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
RWORK — DOUBLE PRECISION array, dimension (N).
EPS3 — DOUBLE PRECISION. A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots. She's nudging the awkward ones and swapping out the zero pivots, like a lil glow-up 💅
SMLNUM — DOUBLE PRECISION. A machine-dependent value close to the underflow threshold.
INFO — INTEGER. The receipts:
- = 0: successful exit 🎉
- = 1: inverse iteration did not converge; V is set to the last iterate. Not a slay, but she tells you honestly 😭

Author — Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver, NAG Ltd.

🌸 subroutine zlaev2 🌸

subroutine zlaev2 (complex*16 A, complex*16 B, complex*16 C, double
    precision RT1, double precision RT2, double precision CS1, complex*16
    SN1)

ZLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. 🦋

Purpose:

ZLAEV2 computes the eigendecomposition of a 2-by-2 Hermitian matrix

   [  A         B  ]
   [  CONJG(B)  C  ].

On return, RT1 is the eigenvalue of larger absolute value, RT2 is the eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right eigenvector for RT1, giving the decomposition

[ CS1  CONJG(SN1) ] [    A     B ] [ CS1 -CONJG(SN1) ] = [ RT1  0  ]
[-SN1     CS1     ] [ CONJG(B) C ] [ SN1     CS1     ]   [  0  RT2 ].

Look at that factorization, she ate and left no crumbs 😌

Parameters

A — COMPLEX*16. The (1,1) element of the 2-by-2 matrix.
B — COMPLEX*16. The (1,2) element and the conjugate of the (2,1) element of the 2-by-2 matrix.
C — COMPLEX*16. The (2,2) element of the 2-by-2 matrix.
RT1 — DOUBLE PRECISION. The eigenvalue of larger absolute value.
RT2 — DOUBLE PRECISION. The eigenvalue of smaller absolute value.
CS1 — DOUBLE PRECISION.
SN1 — COMPLEX*16. The vector (CS1, SN1) is a unit right eigenvector for RT1.

Author — Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver, NAG Ltd.

Further Details:

Bestie, listen carefully, this is the accuracy lore and it matters 💖

RT1 is accurate to a few ulps barring over/underflow.
RT2 may be inaccurate if there is massive cancellation in the determinant A*C-B*B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases. Yeah so like, RT2 can get a lil wobbly, that's just how it is 😭
CS1 and SN1 are accurate to a few ulps barring over/underflow.
Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps.

🍓 subroutine zlag2c 🍓

subroutine zlag2c (integer M, integer N, complex*16, dimension( lda, * ) A,
    integer LDA, complex, dimension( ldsa, * ) SA, integer LDSA, integer
    INFO)

ZLAG2C converts a complex double precision matrix to a complex single precision matrix. ✨

Purpose:

ZLAG2C converts a COMPLEX*16 matrix, SA, to a COMPLEX matrix, A. Downcasting the vibe, keeping the plot 🎀

RMAX is the overflow for the SINGLE PRECISION arithmetic. ZLAG2C checks that all the entries of A are between -RMAX and RMAX. If not, the conversion is aborted and a flag is raised. She does a vibe check first and she will walk out if it fails 🚩

This is an auxiliary routine so there is no argument checking. Yeah, no hand-holding here, she trusts you.

Parameters

M — and okay bestie, this is literally exactly where our part runs out, so hold that thought, M's whole moment is coming up next 💅
M — INTEGER 🌸 how many lines the matrix A has. M >= 0, bestie, always.
N — INTEGER ✨ how many columns the matrix A has. N >= 0.
A — COMPLEX*16 array, dimension (LDA,N) 💜 on entry, this is the M-by-N coefficient matrix A.
LDA — INTEGER — the leading dimension of the array A. LDA >= max(1,M), no exceptions queen.
SA — COMPLEX array, dimension (LDSA,N) 🫶 on exit, if INFO=0, this is the M-by-N coefficient matrix SA; if INFO>0, the content of SA is unspecified (translation: don't trust it, don't look at it, don't text it back).
LDSA — INTEGER — the leading dimension of the array SA. LDSA >= max(1,M).
INFO — INTEGER 👑 the vibe report:
- = 0: successful exit. slay.
- = 1: an entry of the matrix A is greater than the SINGLE PRECISION overflow threshold — and in this case, the content of SA on exit is unspecified. 💀

🫶 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ `subroutine zlags2`

subroutine zlags2 (logical UPPER, double precision A1, complex*16 A2, double
    precision A3, double precision B1, complex*16 B2, double precision B3,
    double precision CSU, complex*16 SNU, double precision CSV, complex*16
    SNV, double precision CSQ, complex*16 SNQ)

ZLAGS2 💅

🌈 Purpose

okay so ZLAGS2 computes 2-by-2 unitary matrices U, V and Q, such that if (UPPER) then

          U**H *A*Q = U**H *( A1 A2 )*Q = ( x  0  )
                            ( 0  A3 )     ( x  x  )
and
          V**H*B*Q = V**H *( B1 B2 )*Q = ( x  0  )
                           ( 0  B3 )     ( x  x  )

or if (.NOT.UPPER) then

          U**H *A*Q = U**H *( A1 0  )*Q = ( x  x  )
                            ( A2 A3 )     ( 0  x  )
and
          V**H *B*Q = V**H *( B1 0  )*Q = ( x  x  )
                            ( B2 B3 )     ( 0  x  )

where

  U = (   CSU    SNU ), V = (  CSV    SNV ),
      ( -SNU**H  CSU )      ( -SNV**H CSV )

  Q = (   CSQ    SNQ )
      ( -SNQ**H  CSQ )

The rows of the transformed A and B are parallel 🦋. Moreover, if the input 2-by-2 matrix A is not zero, then the transformed (1,1) entry of A is not zero. If the input matrices A and B are both not zero, then the transformed (2,2) element of B is not zero — except when the first rows of input A and B are parallel and the second rows are zero. (yes, that specific, that iconic.)

🎀 Parameters

UPPER — LOGICAL
- = .TRUE.: the input matrices A and B are upper triangular.
- = .FALSE.: the input matrices A and B are lower triangular.
A1 — DOUBLE PRECISION
A2 — COMPLEX*16
A3 — DOUBLE PRECISION — on entry, A1, A2 and A3 are elements of the input 2-by-2 upper (lower) triangular matrix A. 🍓
B1 — DOUBLE PRECISION
B2 — COMPLEX*16
B3 — DOUBLE PRECISION — on entry, B1, B2 and B3 are elements of the input 2-by-2 upper (lower) triangular matrix B. 🌸
CSU — DOUBLE PRECISION
SNU — COMPLEX*16 — the desired unitary matrix U. 💜
CSV — DOUBLE PRECISION
SNV — COMPLEX*16 — the desired unitary matrix V. 💖
CSQ — DOUBLE PRECISION
SNQ — COMPLEX*16 — the desired unitary matrix Q. ✨

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 `subroutine zlagtm`

subroutine zlagtm (character TRANS, integer N, integer NRHS, double
    precision ALPHA, complex*16, dimension( * ) DL, complex*16, dimension( *
    ) D, complex*16, dimension( * ) DU, complex*16, dimension( ldx, * ) X,
    integer LDX, double precision BETA, complex*16, dimension( ldb, * ) B,
    integer LDB)

ZLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix, B and C are rectangular matrices, and α and β are scalars, which may be 0, 1, or -1. 🧚

🌈 Purpose

not me about to explain this so it actually clicks — ZLAGTM performs a matrix-vector product of the form

   B := alpha * A * X + beta * B

where A is a tridiagonal matrix of order N, B and X are N by NRHS matrices, and alpha and beta are real scalars, each of which may be 0., 1., or -1. periodt. 💅

🎀 Parameters

TRANS — CHARACTER*1 — specifies the operation applied to A. 🔥
- = 'N': No transpose, B := alpha * A * X + beta * B
- = 'T': Transpose, B := alpha * A**T * X + beta * B
- = 'C': Conjugate transpose, B := alpha * A**H * X + beta * B
N — INTEGER — the order of the matrix A. N >= 0.
NRHS — INTEGER — the number of right hand sides, i.e., the number of columns of the matrices X and B. 💁♀️
ALPHA — DOUBLE PRECISION — the scalar alpha. ALPHA must be 0., 1., or -1.; otherwise, it is assumed to be 0.
DL — COMPLEX*16 array, dimension (N-1) — the (n-1) sub-diagonal elements of T. 🦋
D — COMPLEX*16 array, dimension (N) — the diagonal elements of T.
DU — COMPLEX*16 array, dimension (N-1) — the (n-1) super-diagonal elements of T.
X — COMPLEX*16 array, dimension (LDX,NRHS) — the N by NRHS matrix X. 🌸
LDX — INTEGER — the leading dimension of the array X. LDX >= max(N,1).
BETA — DOUBLE PRECISION — the scalar beta. BETA must be 0., 1., or -1.; otherwise, it is assumed to be 1.
B — COMPLEX*16 array, dimension (LDB,NRHS) — on entry, the N by NRHS matrix B. On exit, B is overwritten by the matrix expression B := alpha * A * X + beta * B. 💖
LDB — INTEGER — the leading dimension of the array B. LDB >= max(N,1).

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ `subroutine zlahqr`

subroutine zlahqr (logical WANTT, logical WANTZ, integer N, integer ILO,
    integer IHI, complex*16, dimension( ldh, * ) H, integer LDH, complex*16,
    dimension( * ) W, integer ILOZ, integer IHIZ, complex*16, dimension(
    ldz, * ) Z, integer LDZ, integer INFO)

ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. 💅

🌈 Purpose

ZLAHQR is an auxiliary routine called by CHSEQR to update the eigenvalues and Schur decomposition already computed by CHSEQR, by dealing with the Hessenberg submatrix in rows and columns ILO to IHI. 🎀

🎀 Parameters

WANTT — LOGICAL
- = .TRUE. : the full Schur form T is required;
- = .FALSE.: only eigenvalues are required.
WANTZ — LOGICAL
- = .TRUE. : the matrix of Schur vectors Z is required;
- = .FALSE.: Schur vectors are not required.
N — INTEGER — the order of the matrix H. N >= 0.
ILO — INTEGER
IHI — INTEGER — it is assumed that H is already upper triangular in rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). ZLAHQR works primarily with the Hessenberg submatrix in rows and columns ILO to IHI, but applies transformations to all of H if WANTT is .TRUE.. 1 <= ILO <= max(1,IHI); IHI <= N. 👑
H — COMPLEX*16 array, dimension (LDH,N) — on entry, the upper Hessenberg matrix H. On exit, if INFO is zero and if WANTT is .TRUE., then H is upper triangular in rows and columns ILO:IHI. If INFO is zero and if WANTT is .FALSE., then the contents of H are unspecified on exit. The output state of H in case INF is positive is below under the description of INFO. 💀
LDH — INTEGER — the leading dimension of the array H. LDH >= max(1,N).
W — COMPLEX*16 array, dimension (N) — the computed eigenvalues ILO to IHI are stored in the corresponding elements of W. If WANTT is .TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i). 🌸
ILOZ — INTEGER
IHIZ — INTEGER — specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
Z — COMPLEX*16 array, dimension (LDZ,N) — if WANTZ is .TRUE., on entry Z must contain the current matrix Z of transformations accumulated by CHSEQR, and on exit Z has been updated; transformations are applied only to the submatrix Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not referenced. 🙈
LDZ — INTEGER — the leading dimension of the array Z. LDZ >= max(1,N).
INFO — INTEGER — the receipts:
- = 0: successful exit 🎉
- > 0: if INFO = i, ZLAHQR failed to compute all the eigenvalues ILO to IHI in a total of 30 iterations per eigenvalue; elements i+1:ihi of W contain those eigenvalues which have been successfully computed. 😭

If INFO > 0 and WANTT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H.

okay so if INFO > 0 and WANTT is .TRUE., then on exit:

(*)       (initial value of H)*U  = U*(final value of H)

where U is an orthogonal matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. 💅

And if INFO > 0 and WANTZ is .TRUE., then on exit:

(final value of Z)  = (initial value of Z)*U

where U is the orthogonal matrix in (*) — and this holds regardless of the value of WANTT. no cap. 🦋

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🎀 Contributors:

02-96 Based on modifications by
David Day, Sandia National Laboratory, USA

12-04 Further modifications by
Ralph Byers, University of Kansas, USA
This is a modified version of ZLAHQR from LAPACK version 3.0.
It is (1) more robust against overflow and underflow and
(2) adopts the more conservative Ahues & Tisseur stopping
criterion (LAWN 122, 1997).

so basically this one got a glow-up from the original ZLAHQR 💖 — more robust against overflow and underflow, and it adopted that conservative Ahues & Tisseur stopping criterion (LAWN 122, 1997). iconic behavior honestly.

✨ subroutine zlahr2 (integer N, integer K, integer NB, complex16, dimension( lda, ) A, integer LDA, complex16, dimension( nb ) TAU, complex16, dimension( ldt, nb ) T, integer LDT, complex*16, dimension( ldy, nb ) Y, integer LDY)

ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. 🫶

Purpose:

ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) matrix A so that elements below the k-th subdiagonal are zero. The reduction is performed by an unitary similarity transformation Q**H * A * Q. The routine returns the matrices V and T which determine Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T.

This is an auxiliary routine called by ZGEHRD. ✨

Parameters

N — N is INTEGER. The order of the matrix A.
K — K is INTEGER. The offset for the reduction. Elements below the k-th subdiagonal in the first NB columns are reduced to zero. K < N.
NB — NB is INTEGER. The number of columns to be reduced.
A — A is COMPLEX*16 array, dimension (LDA,N-K+1). On entry, the n-by-(n-k+1) general matrix A. On exit, the elements on and above the k-th subdiagonal in the first NB columns are overwritten with the corresponding elements of the reduced matrix; the elements below the k-th subdiagonal, with the array TAU, represent the matrix Q as a product of elementary reflectors. The other columns of A are unchanged. See Further Details. 👀
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,N).
TAU — TAU is COMPLEX*16 array, dimension (NB). The scalar factors of the elementary reflectors. See Further Details.
T — T is COMPLEX*16 array, dimension (LDT,NB). The upper triangular matrix T.
LDT — LDT is INTEGER. The leading dimension of the array T. LDT >= NB.
Y — Y is COMPLEX*16 array, dimension (LDY,NB). The n-by-nb matrix Y.
LDY — LDY is INTEGER. The leading dimension of the array Y. LDY >= N.

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Further Details:

The matrix Q is represented as a product of nb elementary reflectors

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

Each H(i) has the form

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

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

The elements of the vectors v together form the (n-k+1)-by-nb matrix V which is needed, with T and Y, to apply the transformation to the unreduced part of the matrix, using an update of the form: A := (I - V*T*V**H) * (A - Y*V**H).

The contents of A on exit are illustrated by the following example with n = 7, k = 3 and nb = 2:

   ( a   a   a   a   a )
   ( a   a   a   a   a )
   ( a   a   a   a   a )
   ( h   h   a   a   a )
   ( v1  h   a   a   a )
   ( v1  v2  a   a   a )
   ( v1  v2  a   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). 🌸

This subroutine is a slight modification of LAPACK-3.0's ZLAHRD incorporating improvements proposed by Quintana-Orti and Van de Gejin. Note that the entries of A(1:K,2:NB) differ from those returned by the original LAPACK-3.0's ZLAHRD routine. (This subroutine is not backward compatible with LAPACK-3.0's ZLAHRD.) ⚠️ bestie, do not assume it swaps in cleanly.

References:

Gregorio Quintana-Orti and Robert van de Geijn, 'Improving the performance of reduction to Hessenberg form,' ACM Transactions on Mathematical Software, 32(2):180-194, June 2006.

💜 subroutine zlaic1 (integer JOB, integer J, complex16, dimension( j ) X, double precision SEST, complex16, dimension( j ) W, complex16 GAMMA, double precision SESTPR, complex16 S, complex*16 C)

ZLAIC1 applies one step of incremental condition estimation. ✨

Purpose:

ZLAIC1 applies one step of incremental condition estimation in its simplest version:

Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j lower triangular matrix L, such that

twonorm(L*x) = sest

Then ZLAIC1 computes sestpr, s, c such that the vector

         [ s*x ]
xhat =   [  c  ]

is an approximate singular vector of

         [ L       0  ]
Lhat =   [ w**H gamma ]

in the sense that

twonorm(Lhat*xhat) = sestpr.

Depending on JOB, an estimate for the largest or smallest singular value is computed. 👑

Note that [s c]**H and sestpr**2 is an eigenpair of the system

diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]
                                      [ conjg(gamma) ]

where alpha = x**H * w.

Parameters

JOB — JOB is INTEGER. = 1: an estimate for the largest singular value is computed. = 2: an estimate for the smallest singular value is computed.
J — J is INTEGER. Length of X and W.
X — X is COMPLEX*16 array, dimension (J). The j-vector x.
SEST — SEST is DOUBLE PRECISION. Estimated singular value of j by j matrix L.
W — W is COMPLEX*16 array, dimension (J). The j-vector w.
GAMMA — GAMMA is COMPLEX*16. The diagonal element gamma.
SESTPR — SESTPR is DOUBLE PRECISION. Estimated singular value of (j+1) by (j+1) matrix Lhat.
S — S is COMPLEX*16. Sine needed in forming xhat.
C — C is COMPLEX*16. Cosine needed in forming xhat.

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🦋 double precision function zlangt (character NORM, integer N, complex16, dimension( ) DL, complex16, dimension( ) D, complex16, dimension( ) DU)

ZLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. 🌈

Purpose:

ZLANGT returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex tridiagonal matrix A.

Returns

ZLANGT

ZLANGT = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. 💅

Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in ZLANGT as described above.
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, ZLANGT is set to zero.
DL — DL is COMPLEX*16 array, dimension (N-1). The (n-1) sub-diagonal elements of A.
D — D is COMPLEX*16 array, dimension (N). The diagonal elements of A.
DU — DU is COMPLEX*16 array, dimension (N-1). The (n-1) super-diagonal elements of A.

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🎀 double precision function zlanhb (character NORM, character UPLO, integer N, integer K, complex16, dimension( ldab, ) AB, integer LDAB, double precision, dimension( * ) WORK)

ZLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix. 🫶

Purpose:

ZLANHB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n hermitian band matrix A, with k super-diagonals.

Returns

ZLANHB

ZLANHB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

okay so quick translation bestie 💖 norm1 is the one norm of a matrix, aka the maximum column sum. normI is the infinity norm, aka the maximum row sum. And normF is the Frobenius norm, the square root of sum of squares. and the gag is — max(abs(A(i,j))) is not a consistent matrix norm. it's giving "technically works, mathematically messy". 💀

✨ Parameters ✨

NORM — NORM is CHARACTER*1. it tells you which value gets handed back in ZLANHB, exactly like described up above. no surprises here. 🎀
UPLO — UPLO is CHARACTER*1. it specifies whether the upper or lower triangular part of the band matrix A is the thing you're supplying.
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. it's the order of the matrix A. N >= 0. and when N = 0, ZLANHB gets set to zero. slay, edge case handled. ✨
K — K is INTEGER. the number of super-diagonals or sub-diagonals of the band matrix A. K >= 0. 🌸
AB — AB is COMPLEX*16 array, dimension (LDAB,N). it holds the upper or lower triangle of the hermitian band matrix A, stored in the first K+1 rows of AB. the j-th column of A lives in the j-th column of the array AB, like this:
- if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
- if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+k).

heads up though — the imaginary parts of the diagonal elements don't need to be set, they're just assumed to be zero. so don't stress about them bestie, they're lowkey not your problem. 💅

LDAB — LDAB is INTEGER. the leading dimension of the array AB. LDAB >= K+1. periodt.
WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'. otherwise, WORK is not referenced. it's giving "you can ignore me, I'm just here for moral support". 🙈

🌸 Author 🌸

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🦋 `double precision function zlanhp` 🦋

double precision function zlanhp (character NORM, character UPLO, integer N,
    complex*16, dimension( * ) AP, double precision, dimension( * ) WORK)

ZLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form. packed!! it's like the matrix got vacuum-sealed for travel. 🧳✨

Purpose:

ZLANHP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex hermitian matrix A, supplied in packed form.

Returns:

ZLANHP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

Parameters:

NORM — NORM is CHARACTER*1. specifies the value to be returned in ZLANHP as described above. iykyk. 💜
UPLO — UPLO is CHARACTER*1. specifies whether the upper or lower triangular part of the hermitian matrix A is supplied.
- = 'U': Upper triangular part of A is supplied
- = 'L': Lower triangular part of A is supplied
N — N is INTEGER. the order of the matrix A. N >= 0. when N = 0, ZLANHP is set to zero. 💖
AP — AP is COMPLEX*16 array, dimension (N*(N+1)/2). the upper or lower triangle of the hermitian matrix A, packed columnwise in a linear array. the j-th column of A is stored in the array AP as follows:
- if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
- if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

and again: the imaginary parts of the diagonal elements need not be set and are assumed to be zero. 💅

WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced. 🎀

Author:

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

⭐ `double precision function zlanhs` ⭐

double precision function zlanhs (character NORM, integer N, complex*16,
    dimension( lda, * ) A, integer LDA, double precision, dimension( * )
    WORK)

ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix. Hessenberg!! she's almost triangular but she keeps ONE little extra diagonal. main character energy honestly. 👑

Purpose:

ZLANHS returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hessenberg matrix A.

Returns:

ZLANHS = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

Parameters:

NORM — NORM is CHARACTER*1. specifies the value to be returned in ZLANHS as described above. 💫
N — N is INTEGER. the order of the matrix A. N >= 0. when N = 0, ZLANHS is set to zero. ✨
A — A is COMPLEX*16 array, dimension (LDA,N). the n by n upper Hessenberg matrix A. the part of A below the first sub-diagonal is not referenced — so like, don't even bother filling it in, nobody's reading it. 🙈
LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= max(N,1). 🌸
WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced. note bestie: this one only needs WORK for 'I', not for '1' or 'O' like the others. she's different like that. 💅

Author:

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🎀 `double precision function zlanht` 🎀

double precision function zlanht (character NORM, integer N, double
    precision, dimension( * ) D, complex*16, dimension( * ) E)

ZLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix. tridiagonal!! only three diagonals, she's minimal, she's chic, she ate. 🍽️✨

Purpose:

ZLANHT returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix A.

Returns:

ZLANHT = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

Parameters:

NORM — NORM is CHARACTER*1. specifies the value to be returned in ZLANHT as described above. 🫶
N — N is INTEGER. the order of the matrix A. N >= 0. when N = 0, ZLANHT is set to zero. 💖
D — D is DOUBLE PRECISION array, dimension (N). the diagonal elements of A. the moment. the center of attention. ⭐
E — E is COMPLEX*16 array, dimension (N-1). the (n-1) sub-diagonal or super-diagonal elements of A. yeah she's one shorter than D, that's not a typo, that's just how tridiagonal works bestie. 🌸

Author:

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

💜 `double precision function zlansb` 💜

double precision function zlansb (character NORM, character UPLO, integer N,
    integer K, complex*16, dimension( ldab, * ) AB, integer LDAB, double
    precision, dimension( * ) WORK)

ZLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix. band matrix squad, assemble. 🎸

Purpose:

ZLANSB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n symmetric band matrix A, with k super-diagonals.

Returns:

ZLANSB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

Parameters:

NORM — NORM is CHARACTER*1. specifies the value to be returned in ZLANSB as described above. ✨
UPLO — UPLO is CHARACTER*1. specifies whether the upper or lower triangular part of the band matrix A is supplied.
- = 'U': Upper triangular part is supplied
- = 'L': Lower triangular part is supplied
N — N is INTEGER. the order of the matrix A. N >= 0. when N = 0, ZLANSB is set to zero. 🎀
K
K — INTEGER, bestie. It's how many super-diagonals or sub-diagonals the band matrix A has going on. Requirement: K >= 0. No negatives, periodt. 💅
AB — COMPLEX*16 array, dimension (LDAB,N). This is the upper or lower triangle of the symmetric band matrix A, and it lives in the first K+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB like this:

if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)   = A(i,j) for j<=i<=min(n,j+k).

LDAB — INTEGER, the leading dimension of the array AB. Gotta be LDAB >= K+1, no cap. ✨
WORK — DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'. Otherwise? WORK is literally not referenced. She said what she said. 🫶

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 `double precision function zlansp`

double precision function zlansp (character NORM, character UPLO, integer N,
    complex*16, dimension( * ) AP, double precision, dimension( * ) WORK)

ZLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form. It's giving "I only do one job and I do it flawlessly." 💖

✨ Purpose

ZLANSP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix A, supplied in packed form.

🎀 Returns

ZLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Heads up though — max(abs(A(i,j))) is not a consistent matrix norm. Not me exposing the math like that 😭

💜 Parameters

NORM — CHARACTER*1. Specifies the value to be returned in ZLANSP as described above. Pick your fighter. 🌈
UPLO — CHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrix A is supplied.
- = 'U': Upper triangular part of A is supplied
- = 'L': Lower triangular part of A is supplied
N — INTEGER. The order of the matrix A. N >= 0. When N = 0, ZLANSP is set to zero. Empty matrix? Zero. Slay. 💅
AP — COMPLEX*16 array, dimension (N*(N+1)/2). The upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows:

if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

WORK — DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced. Same energy as before, she's consistent. 🦋

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 `double precision function zlantb`

double precision function zlantb (character NORM, character UPLO, character
    DIAG, integer N, integer K, complex*16, dimension( ldab, * ) AB, integer
    LDAB, double precision, dimension( * ) WORK)

ZLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular band matrix. Band matrix but make it triangular. Iconic. ✨

✨ Purpose

ZLANTB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n triangular band matrix A, with ( k + 1 ) diagonals.

🎀 Returns

ZLANTB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And again, fr fr: max(abs(A(i,j))) is not a consistent matrix norm. 🫶

💜 Parameters

NORM — CHARACTER*1. Specifies the value to be returned in ZLANTB as described above.
UPLO — CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
DIAG — CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
N — INTEGER. The order of the matrix A. N >= 0. When N = 0, ZLANTB is set to zero. 💅
K — INTEGER. The number of super-diagonals of the matrix A if UPLO = 'U', or the number of sub-diagonals of the matrix A if UPLO = 'L'. K >= 0.
AB — COMPLEX*16 array, dimension (LDAB,N). The upper or lower triangular band matrix A, stored in the first k+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB as follows:

if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)   = A(i,j) for j<=i<=min(n,j+k).

Note that when DIAG = 'U', the elements of the array AB corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one. She's not reading them, she's just assuming. 🙈

LDAB — INTEGER. The leading dimension of the array AB. LDAB >= K+1.
WORK — DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced. Note it's only 'I' this time, bestie, watch the difference. 👀

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 `double precision function zlantp`

double precision function zlantp (character NORM, character UPLO, character
    DIAG, integer N, complex*16, dimension( * ) AP, double precision,
    dimension( * ) WORK)

ZLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. Triangular + packed = little storage, big power. 💖

✨ Purpose

ZLANTP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix A, supplied in packed form.

🎀 Returns

ZLANTP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And yes, one more time for the people in the back: max(abs(A(i,j))) is not a consistent matrix norm. 🌈

💜 Parameters

NORM — CHARACTER*1. Specifies the value to be returned in ZLANTP as described above.
UPLO — CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
DIAG — CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
N — INTEGER. The order of the matrix A. N >= 0. When N = 0, ZLANTP is set to zero. Zero. Done. ✨
AP — COMPLEX*16 array, dimension (N*(N+1)/2). The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows:

if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

Note that when DIAG = 'U', the elements of the array AP corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one. Again — not read, just vibes (and ones). 🙈

WORK — DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 `double precision function zlantr`

double precision function zlantr (character NORM, character UPLO, character
    DIAG, integer M, integer N, complex*16, dimension( lda, * ) A, integer
    LDA, double precision, dimension( * ) WORK)

ZLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix. Trapezoidal?? She contains multitudes. 👑

✨ Purpose

ZLANTR returns 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.

🎀 Returns

ZLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'
         (
         ( norm1(A),         NORM = '1', 'O' or 'o'
         (
         ( normI(A),         NORM = 'I' or 'i'
         (
         ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And (say it with me now 🎤) max(abs(A(i,j))) is not a consistent matrix norm. 💜

Parameters ✨

NORM — NORM is CHARACTER*1. Specifies the value to be returned in ZLANTR as described above.
UPLO — UPLO is CHARACTER*1. Specifies whether the matrix A is upper or lower trapezoidal.
- = 'U': Upper trapezoidal
- = 'L': Lower trapezoidal Note that A is triangular instead of trapezoidal if M = N.
DIAG — DIAG is CHARACTER*1. Specifies whether or not the matrix A has unit diagonal.
- = 'N': Non-unit diagonal
- = 'U': Unit diagonal
M — M is INTEGER. The number of rows of the matrix A. M >= 0, and if UPLO = 'U', M <= N. When M = 0, ZLANTR is set to zero.
N — N is INTEGER. The number of columns of the matrix A. N >= 0, and if UPLO = 'L', N <= M. When N = 0, ZLANTR is set to zero.
A — A is COMPLEX*16 array, dimension (LDA,N). The trapezoidal matrix A (A is triangular if M = N). If UPLO = 'U', 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 = 'L', 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. Note that when DIAG = 'U', the diagonal elements of A are not referenced and are assumed to be one.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(M,1).
WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= M when NORM = 'I'; otherwise, WORK is not referenced.

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

`zlapll` — measuring the linear dependence of two vectors 💅

subroutine zlapll (integer N, complex*16, dimension( * ) X, integer INCX,
    complex*16, dimension( * ) Y, integer INCY, double precision SSMIN)

ZLAPLL measures the linear dependence of two vectors.

Purpose ✨

Given two column vectors X and Y, let

A = ( X Y ).

The subroutine first computes the QR factorization of A = Q*R, and then computes the SVD of the 2-by-2 upper triangular matrix R. The smaller singular value of R is returned in SSMIN, which is used as the measurement of the linear dependency of the vectors X and Y.

Parameters 💖

N — N is INTEGER. The length of the vectors X and Y.
X — X is COMPLEX*16 array, dimension (1+(N-1)*INCX). On entry, X contains the N-vector X. On exit, X is overwritten.
INCX — INCX is INTEGER. The increment between successive elements of X. INCX > 0.
Y — Y is COMPLEX*16 array, dimension (1+(N-1)*INCY). On entry, Y contains the N-vector Y. On exit, Y is overwritten.
INCY — INCY is INTEGER. The increment between successive elements of Y. INCY > 0.
SSMIN — SSMIN is DOUBLE PRECISION. The smallest singular value of the N-by-2 matrix A = ( X Y ).

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

`zlapmr` — rearranging rows like a boss 👑

subroutine zlapmr (logical FORWRD, integer M, integer N, complex*16,
    dimension( ldx, * ) X, integer LDX, integer, dimension( * ) K)

ZLAPMR rearranges rows of a matrix as specified by a permutation vector.

Purpose ✨

ZLAPMR rearranges the rows of the M by N matrix X as specified by the permutation K(1),K(2),...,K(M) of the integers 1,...,M. If FORWRD = .TRUE., forward permutation:

X(K(I),*) is moved X(I,*) for I = 1,2,...,M.

If FORWRD = .FALSE., backward permutation:

X(I,*) is moved to X(K(I),*) for I = 1,2,...,M.

Parameters 💖

FORWRD — FORWRD is LOGICAL.
- = .TRUE., forward permutation
- = .FALSE., backward permutation
M — M is INTEGER. The number of rows of the matrix X. M >= 0.
N — N is INTEGER. The number of columns of the matrix X. N >= 0.
X — X is COMPLEX*16 array, dimension (LDX,N). On entry, the M by N matrix X. On exit, X contains the permuted matrix X.
LDX — LDX is INTEGER. The leading dimension of the array X, LDX >= MAX(1,M).
K — K is INTEGER array, dimension (M). On entry, K contains the permutation vector. K is used as internal workspace, but reset to its original value on output.

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

`zlapmt` — permuting columns, periodt 💅

subroutine zlapmt (logical FORWRD, integer M, integer N, complex*16,
    dimension( ldx, * ) X, integer LDX, integer, dimension( * ) K)

ZLAPMT performs a forward or backward permutation of the columns of a matrix.

Purpose ✨

ZLAPMT rearranges the columns of the M by N matrix X as specified by the permutation K(1),K(2),...,K(N) of the integers 1,...,N. If FORWRD = .TRUE., forward permutation:

X(*,K(J)) is moved X(*,J) for J = 1,2,...,N.

If FORWRD = .FALSE., backward permutation:

X(*,J) is moved to X(*,K(J)) for J = 1,2,...,N.

Parameters 💖

FORWRD — FORWRD is LOGICAL.
- = .TRUE., forward permutation
- = .FALSE., backward permutation
M — M is INTEGER. The number of rows of the matrix X. M >= 0.
N — N is INTEGER. The number of columns of the matrix X. N >= 0.
X — X is COMPLEX*16 array, dimension (LDX,N). On entry, the M by N matrix X. On exit, X contains the permuted matrix X.
LDX — LDX is INTEGER. The leading dimension of the array X, LDX >= MAX(1,M).
K — K is INTEGER array, dimension (N). On entry, K contains the permutation vector. K is used as internal workspace, but reset to its original value on output.

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

`zlaqhb` — scaling a Hermitian band matrix, no cap 💅

subroutine zlaqhb (character UPLO, integer N, integer KD, complex*16,
    dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * )
    S, double precision SCOND, double precision AMAX, character EQUED)

ZLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ.

Purpose ✨

ZLAQHB equilibrates a Hermitian band matrix A using the scaling factors in the vector S.

Parameters 💖

UPLO — UPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrix A is stored.
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. The order of the matrix A. N >= 0.
KD — KD is INTEGER. The number of super-diagonals of the matrix A if UPLO = 'U', or the number of sub-diagonals if UPLO = 'L'. KD >= 0.
AB — AB is COMPLEX*16 array, dimension (LDAB,N). On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:
- if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
- if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**H *U or A = L*L**H of the band matrix A, in the same storage format as A.
LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= KD+1.
S — S is DOUBLE PRECISION array, dimension (N). The scale factors for A.
SCOND — SCOND is DOUBLE PRECISION. Ratio of the smallest S(i) to the largest S(i).
AMAX — AMAX is DOUBLE PRECISION. Absolute value of largest matrix entry.
EQUED — EQUED is CHARACTER*1. Specifies whether or not equilibration was done.
- = 'N': No equilibration.
- = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).

Internal Parameters 🔍

THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done.
LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

`zlaqhp` — scaling a Hermitian matrix in packed form 💅

subroutine zlaqhp (character UPLO, integer N, complex*16, dimension( * ) AP,
    double precision, dimension( * ) S, double precision SCOND, double
    precision AMAX, character EQUED)

ZLAQHP scales a Hermitian matrix stored in packed form.

Purpose ✨

ZLAQHP equilibrates a Hermitian matrix A using the scaling factors in the vector S.

Parameters 💖

UPLO — UPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the Hermitian matrix A is stored.
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. The order of the matrix A. N >= 0.
AP —
AP — COMPLEX*16 array, dimension (N*(N+1)/2). On entry, bestie, this is the upper or lower triangle of the Hermitian matrix A, packed columnwise into one long linear array. The j-th column of A sits inside AP like this:

if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

And on exit? It's the equilibrated matrix: diag(S) * A * diag(S), chilling in the exact same storage format as A. Obsessed. 💖

S — DOUBLE PRECISION array, dimension (N). The scale factors for A. That's it, that's the whole job. ✨
SCOND — DOUBLE PRECISION. The ratio of the smallest S(i) to the largest S(i). Tiny number = spicy scaling situation, iykyk. 🫶
AMAX — DOUBLE PRECISION. The absolute value of the largest matrix entry. Big mood, literally. 💅
EQUED — CHARACTER*1. Specifies whether or not equilibration was done:

= 'N':  No equilibration.
= 'Y':  Equilibration was done, i.e., A has been replaced by
        diag(S) * A * diag(S).

✨ Internal Parameters ✨

THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done. 🎀

LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done. 🌸

🫶 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `zlaqp2` (`integer M`, `integer N`, `integer OFFSET`, `complex16, dimension( lda, ) A`, `integer LDA`, `integer, dimension( * ) JPVT`, `complex16, dimension( ) TAU`, `double precision, dimension( * ) VN1`, `double precision, dimension( * ) VN2`, `complex16, dimension( ) WORK`)

ZLAQP2 computes a QR factorization with column pivoting of the matrix block. 👑

Purpose:

ZLAQP2 computes a QR factorization with column pivoting of the block A(OFFSET+1:M,1:N). The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. Not me leaving half the block un-factorized bestie, that's literally the spec. 💀

Parameters 🎀

M — INTEGER. The number of rows of the matrix A. M >= 0. ✅
N — INTEGER. The number of columns of the matrix A. N >= 0. ✅
OFFSET — INTEGER. The number of rows of the matrix A that must be pivoted but no factorized. OFFSET >= 0. (Yes, "no factorized" — that's verbatim, I'm not touching it 🫶)
A — COMPLEX*16 array, dimension (LDA,N). On entry, the M-by-N matrix A. On exit, the upper triangle of block A(OFFSET+1:M,1:N) is the triangular factor obtained; the elements in block A(OFFSET+1:M,1:N) below the diagonal, together with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors. Block A(1:OFFSET,1:N) has been accordingly pivoted, but no factorized. ✨
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). 🌸
JPVT — INTEGER array, dimension (N). On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted to the front of A*P (a leading column); if JPVT(i) = 0, the i-th column of A is a free column. On exit, if JPVT(i) = k, then the i-th column of A*P was the k-th column of A. 💅
TAU — COMPLEX*16 array, dimension (min(M,N)). The scalar factors of the elementary reflectors. 🦋
VN1 — DOUBLE PRECISION array, dimension (N). The vector with the partial column norms. 👑
VN2 — DOUBLE PRECISION array, dimension (N). The vector with the exact column norms. 💖
WORK — COMPLEX*16 array, dimension (N). 🎀

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

X. Sun, Computer Science Dept., Duke University, USA

Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

References:

LAPACK Working Note 176

subroutine `zlaqps` (`integer M`, `integer N`, `integer OFFSET`, `integer NB`, `integer KB`, `complex16, dimension( lda, ) A`, `integer LDA`, `integer, dimension( * ) JPVT`, `complex16, dimension( ) TAU`, `double precision, dimension( * ) VN1`, `double precision, dimension( * ) VN2`, `complex16, dimension( ) AUXV`, `complex16, dimension( ldf, ) F`, `integer LDF`)

ZLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. ✨ (the summary says "real" but this is the complex routine — I'm just reporting what it says, don't come for me 💀)

Purpose:

ZLAQPS computes a step of QR factorization with column pivoting of a complex M-by-N matrix A by using Blas-3. It tries to factorize NB columns from A starting from the row OFFSET+1, and updates all of the matrix with Blas-3 xGEMM. 🔥

In some cases, due to catastrophic cancellations, it cannot factorize NB columns. Hence, the actual number of factorized columns is returned in KB. catastrophic cancellations is such a dramatic phrase and honestly? valid. 😭

Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.

Parameters 🎀

M — INTEGER. The number of rows of the matrix A. M >= 0. ✅
N — INTEGER. The number of columns of the matrix A. N >= 0. ✅
OFFSET — INTEGER. The number of rows of A that have been factorized in previous steps. 🫶
NB — INTEGER. The number of columns to factorize. 💅
KB — INTEGER. The number of columns actually factorized. (Which may be less than NB! The cancellation thing. You've been warned.) 💀
A — COMPLEX*16 array, dimension (LDA,N). On entry, the M-by-N matrix A. On exit, block A(OFFSET+1:M,1:KB) is the triangular factor obtained and block A(1:OFFSET,1:N) has been accordingly pivoted, but no factorized. The rest of the matrix, block A(OFFSET+1:M,KB+1:N), has been updated. ✨
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). 🌸
JPVT — INTEGER array, dimension (N). JPVT(I) = K <==> Column K of the full matrix A has been permuted into position I in AP. 🦋
TAU — COMPLEX*16 array, dimension (KB). The scalar factors of the elementary reflectors. 💖
VN1 — DOUBLE PRECISION array, dimension (N). The vector with the partial column norms. 👑
VN2 — DOUBLE PRECISION array, dimension (N). The vector with the exact column norms. ⭐
AUXV — COMPLEX*16 array, dimension (NB). Auxiliary vector. 🎀
F — COMPLEX*16 array, dimension (LDF,NB). Matrix F**H = L * Y**H * A. 🔥
LDF — INTEGER. The leading dimension of the array F. LDF >= max(1,N). 💜

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

X. Sun, Computer Science Dept., Duke University, USA

Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

References:

LAPACK Working Note 176

subroutine `zlaqr0` (`logical WANTT`, `logical WANTZ`, `integer N`, `integer ILO`, `integer IHI`, `complex16, dimension( ldh, ) H`, `integer LDH`, `complex16, dimension( ) W`, `integer ILOZ`, `integer IHIZ`, `complex16, dimension( ldz, ) Z`, `integer LDZ`, `complex16, dimension( ) WORK`, `integer LWORK`, `integer INFO`)

ZLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. 👑✨ this one is kind of the main character ngl

Purpose:

ZLAQR0 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**H, where T is an upper triangular matrix (the Schur form), and Z is the unitary matrix of Schur vectors. 💅

Optionally Z may be postmultiplied into an input unitary 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 unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H. 🦋

Parameters 🎀

WANTT — LOGICAL:

= .TRUE. : the full Schur form T is required;
= .FALSE.: only eigenvalues are required.

WANTZ — LOGICAL:

= .TRUE. : the matrix of Schur vectors Z is required;
= .FALSE.: Schur vectors are not required.

N — INTEGER. The order of the matrix H. N >= 0. ✅
ILO — INTEGER. 🫶
IHI — INTEGER. It is assumed that H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N and, if ILO > 1, H(ILO,ILO-1) is zero. ILO and IHI are normally set by a previous call to ZGEBAL, and then passed to ZGEHRD when the matrix output by ZGEBAL is reduced to Hessenberg form. Otherwise, ILO and IHI should be set to 1 and N, respectively. If N > 0, then 1 <= ILO <= IHI <= N. If N = 0, then ILO = 1 and IHI = 0. 💖
H — COMPLEX*16 array, dimension (LDH,N). On entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper triangular matrix T from the Schur decomposition (the Schur form). If INFO = 0 and WANT is .FALSE., then the contents of H are unspecified on exit. (The output value of H when INFO > 0 is given under the description of INFO below.) 😭

This subroutine may explicitly set H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. ✨

LDH —
LDH — INTEGER. The leading dimension of the array H. LDH >= max(1,N). bestie it's just how wide the array is declared, don't overthink it 💅
W — COMPLEX*16 array, dimension (N). The computed eigenvalues of H(ILO:IHI,ILO:IHI) get parked in W(ILO:IHI). And if WANTT is .TRUE., then those eigenvalues come out in the same order as they sit on the diagonal of the Schur form returned in H, with W(i) = H(i,i). cute, consistent, iconic ✨
ILOZ — INTEGER. (no description given here, it's just vibing as an integer 😌)
IHIZ — INTEGER. This one specifies the rows of Z that transformations have to be applied to when WANTZ is .TRUE. — like, the exact rows, no freelancing. Constraint check: 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
Z — COMPLEX*16 array, dimension (LDZ,IHI). If WANTZ is .FALSE., then Z is not referenced at all, she's off duty 🙈. If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the orthogonal Schur factor of H(ILO:IHI,ILO:IHI). And heads up: the output value of Z when INFO > 0 is described down in the INFO section, so scroll, don't guess 🫶
LDZ — INTEGER. The leading dimension of the array Z. If WANTZ is .TRUE. then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1. Simple as that.
WORK — COMPLEX*16 array, dimension LWORK. On exit, if LWORK = -1, then WORK(1) returns an estimate of the optimal value for LWORK. yes queen, she reports back 📋
LWORK — INTEGER. The dimension of the array WORK. LWORK >= max(1,N) is sufficient, but honestly LWORK typically as large as 6*N may be required for optimal performance. A workspace query to determine the optimal workspace size is recommended, fr fr.

If LWORK = -1, then ZLAQR0 does a workspace query. In this case, ZLAQR0 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. so she peeks without touching, respectful 💅

INFO — INTEGER.
- = 0: successful exit 🎉 periodt.
- > 0: if INFO = i, ZLAQR0 failed to compute all of the eigenvalues. Elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.) like, genuinely rare, don't spiral.

If INFO > 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H.

If INFO > 0 and WANTT is .TRUE., then on exit

(*)  (initial value of H)*U  = U*(final value of H)

where U is a unitary matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI.

If INFO > 0 and WANTZ is .TRUE., then on exit

  (final value of Z(ILO:IHI,ILOZ:IHIZ)
   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U

where U is the unitary matrix in (*) (regardless of the value of WANTT.)

If INFO > 0 and WANTZ is .FALSE., then Z is not accessed. she said not today 🙅‍♀️

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors:

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

📚 References:

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
Performance, SIAM Journal of Matrix Analysis, volume 23, pages
929--947, 2002.

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm
Part II: Aggressive Early Deflation, SIAM Journal of Matrix
Analysis, volume 23, pages 948--973, 2002.

`subroutine zlaqr1 (integer N, complex16, dimension( ldh, ) H, integer LDH, complex16 S1, complex16 S2, complex16, dimension( ) V)`

ZLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. that's the whole job, she's a specialist 💖

Purpose:

Given a 2-by-2 or 3-by-3 matrix H, ZLAQR1 sets v to a scalar multiple of the first column of the product

(*)  K = (H - s1*I)*(H - s2*I)

scaling to avoid overflows and most underflows. safety first, always 🧚

This is useful for starting double implicit shift bulges in the QR algorithm. yes bestie, bulges, it's a whole thing, it's giving main character energy 👑

Parameters

N — INTEGER. Order of the matrix H. N must be either 2 or 3. only those two, she's picky and that's her right 💅
H — COMPLEX*16 array, dimension (LDH,N). The 2-by-2 or 3-by-3 matrix H in (*).
LDH — INTEGER. The leading dimension of H as declared in the calling procedure. LDH >= N.
S1 — COMPLEX*16. (just sitting there being a shift, we love her)
S2 — COMPLEX*16. S1 and S2 are the shifts defining K in (*) above. the duo, the icons, the moment ✨
V — COMPLEX*16 array, dimension (N). A scalar multiple of the first column of the matrix K in (*). there she is 🎀

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

`subroutine zlaqr2 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT, integer NW, complex16, dimension( ldh, ) H, integer LDH, integer ILOZ, integer IHIZ, complex16, dimension( ldz, ) Z, integer LDZ, integer NS, integer ND, complex16, dimension( ) SH, complex16, dimension( ldv, ) V, integer LDV, integer NH, complex16, dimension( ldt, ) T, integer LDT, integer NV, complex16, dimension( ldwv, ) WV, integer LDWV, complex16, dimension( ) WORK, integer LWORK)`

ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). aggressive early deflation is SUCH a name, the audacity, the confidence 🔥

Purpose:

ZLAQR2 is identical to ZLAQR3 except that it avoids recursion by calling ZLAHQR instead of ZLAQR4. noted, no recursion here, she keeps it flat 📏

Aggressive early deflation:

ZLAQR2 accepts as input an upper Hessenberg matrix H and performs an unitary similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been overwritten by a new Hessenberg matrix that is a perturbation of an unitary similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. fingers crossed bestie 🤞 zero subdiagonals are the dream 🌸

Parameters

WANTT — LOGICAL. If .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues. minimal effort, maximum result, that's balance ⭐
WANTZ — LOGICAL. If .TRUE., then the unitary matrix Z is updated so so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced.
N — INTEGER. The order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary matrix Z.
KTOP — INTEGER. It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix. an isolated block, she's independent, she doesn't need anyone 💅
KBOT — INTEGER. It is assumed without a check that either KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix. > yes, "without a check" — she trusts you, don't make it weird 🙈
NW — INTEGER. Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
H — COMPLEX*16 array, dimension (LDH,N). On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries. glow up complete ✨
LDH — INTEGER. Leading dimension of H just as declared in the calling subroutine. N <= LDH.
ILOZ — INTEGER. (she's here, she's an integer, she has no notes 😌)
IHIZ — INTEGER. Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
Z — COMPLEX*16 array, dimension (LDZ,N). IF WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ is .FALSE., then Z is unreferenced. untouched, unbothered, thriving 🦋
LDZ — INTEGER. The leading dimension of Z just as declared in the calling subroutine. 1 <= LDZ.
NS — INTEGER. The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine.
ND — INTEGER. The number of converged eigenvalues uncovered by this subroutine. we love an uncover 💖
SH — COMPLEX*16 array, dimension (KBOT). On output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT-ND). Converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT).
V — COMPLEX*16 array, dimension (LDV,NW). An NW-by-NW work array.
LDV —
LDV — INTEGER bestie, it's just the leading dimension of V, exactly as you declared it in the calling subroutine. So like… NW <= LDV. 💅
NH — INTEGER. The number of columns of T. And NH >= NW, no exceptions. 🎀
T — COMPLEX*16 array, dimension (LDT,NW). That's it, that's the girl. ✨
LDT — INTEGER. The leading dimension of T just as declared in the calling subroutine. NW <= LDT. 🦋
NV — INTEGER. How many rows of the work array WV you've got available for workspace. NV >= NW. 🌸
WV — COMPLEX*16 array, dimension (LDWV,NW). Living her life. 💖
LDWV — INTEGER. The leading dimension of W just as declared in the calling subroutine. NW <= LDV. (yes it says LDV, we're keeping it exactly, don't @ me 🙈)
WORK — COMPLEX*16 array, dimension (LWORK). On exit, WORK(1) is set to an estimate of the optimal value of LWORK for the given values of N, NW, KTOP and KBOT. 🔥
LWORK — INTEGER. The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK.

And if LWORK = -1?? Then a workspace query is assumed, okay? ZLAQR2 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. Not even a little. 😭

✨ Author ✨

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🫶 Contributors: 🫶

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

💜 subroutine `zlaqr3` 💜

subroutine zlaqr3 (logical WANTT, logical WANTZ, integer N, integer KTOP,
    integer KBOT, integer NW, complex*16, dimension( ldh, * ) H, integer
    LDH, integer ILOZ, integer IHIZ, complex*16, dimension( ldz, * ) Z,
    integer LDZ, integer NS, integer ND, complex*16, dimension( * ) SH,
    complex*16, dimension( ldv, * ) V, integer LDV, integer NH, complex*16,
    dimension( ldt, * ) T, integer LDT, integer NV, complex*16, dimension(
    ldwv, * ) WV, integer LDWV, complex*16, dimension( * ) WORK, integer
    LWORK)

ZLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). Iconic behaviour, honestly. 🌈

Purpose:

Aggressive early deflation, bestie. Let's go. 💅

ZLAQR3 accepts as input an upper Hessenberg matrix H and performs an unitary similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been over-written by a new Hessenberg matrix that is a perturbation of an unitary similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. (manifesting zero subdiagonals for u 🕯️)

Parameters

WANTT — LOGICAL. If .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues. 🎀
WANTZ — LOGICAL. If .TRUE., then the unitary matrix Z is updated so so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced. (yes, "so so" — that's the original, we don't edit ✨)
N — INTEGER. The order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary matrix Z. 👑
KTOP — INTEGER. It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix. 🌸
KBOT — INTEGER. It is assumed without a check that either KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix. 🦋
NW — INTEGER. Deflation window size. 1 <= NW <= (KBOT-KTOP+1). 💖
H — COMPLEX*16 array, dimension (LDH,N). On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries. 🔥
LDH — INTEGER. Leading dimension of H just as declared in the calling subroutine. N <= LDH. ✨
ILOZ — INTEGER. (that's the whole description, she's mysterious 💅)
IHIZ — INTEGER. Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N. 🎀
Z — COMPLEX*16 array, dimension (LDZ,N). IF WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ is .FALSE., then Z is unreferenced. 🌸
LDZ — INTEGER. The leading dimension of Z just as declared in the calling subroutine. 1 <= LDZ. 🦋
NS — INTEGER. The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine. 💜
ND — INTEGER. The number of converged eigenvalues uncovered by this subroutine. 🔥
SH — COMPLEX*16 array, dimension (KBOT). On output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT). Converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT). (yep, SR — original, verbatim, keeping it 🙈)
V — COMPLEX*16 array, dimension (LDV,NW). An NW-by-NW work array. 💅
LDV — INTEGER. The leading dimension of V just as declared in the calling subroutine. NW <= LDV. ✨
NH — INTEGER. The number of columns of T. NH >= NW. 🎀
T — COMPLEX*16 array, dimension (LDT,NW). 🌸
LDT — INTEGER. The leading dimension of T just as declared in the calling subroutine. NW <= LDT. 🦋
NV — INTEGER. The number of rows of work array WV available for workspace. NV >= NW. 💖
WV — COMPLEX*16 array, dimension (LDWV,NW). 💜
LDWV — INTEGER. The leading dimension of W just as declared in the calling subroutine. NW <= LDV. (again with the LDV, we love consistency 💀)
WORK — COMPLEX*16 array, dimension (LWORK). On exit, WORK(1) is set to an estimate of the optimal value of LWORK for the given values of N, NW, KTOP and KBOT. 🔥
LWORK — INTEGER. The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK.

If LWORK = -1, then a workspace query is assumed; ZLAQR3 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. She is simply ✨not looking✨. 😭

✨ Author ✨

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🫶 Contributors: 🫶

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

👑 subroutine `zlaqr4` 👑

subroutine zlaqr4 (logical WANTT, logical WANTZ, integer N, integer ILO,
    integer IHI, complex*16, dimension( ldh, * ) H, integer LDH, complex*16,
    dimension( * ) W, integer ILOZ, integer IHIZ, complex*16, dimension(
    ldz, * ) Z, integer LDZ, complex*16, dimension( * ) WORK, integer LWORK,
    integer INFO)

ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. Main character energy fr fr. 🌈

Purpose:

ZLAQR4 implements one level of recursion for ZLAQR0. It is a complete implementation of the small bulge multi-shift QR algorithm. It may be called by ZLAQR0 and, for large enough deflation window size, it may be called by ZLAQR3. This subroutine is identical to ZLAQR0 except that it calls ZLAQR2 instead of ZLAQR3. (sibling behaviour, lowkey 💅)

ZLAQR4 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**H, where T is an upper triangular matrix (the Schur form), and Z is the unitary matrix of Schur vectors. 🎀

Parameters

WANTT — LOGICAL.
- = .TRUE. : the full Schur form T is required;
- = .FALSE.: only eigenvalues are required. 💖
WANTZ — LOGICAL.
- = .TRUE. : the matrix of Schur vectors Z is required;
- = .FALSE.: Schur vectors are not required. 🌸
N — INTEGER. The order of the matrix H. N >= 0. 🦋
ILO — INTEGER. (no further description, she simply is 💜)
IHI —
IHI is INTEGER — okay bestie, here's the vibe: it assumes H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N, and if ILO > 1, then H(ILO,ILO-1) is zero. ILO and IHI are normally set by a previous call to ZGEBAL, and then passed to ZGEHRD when the matrix output by ZGEBAL is reduced to Hessenberg form. Otherwise, ILO and IHI should be set to 1 and N, respectively. If N > 0, then 1 <= ILO <= IHI <= N. If N = 0, then ILO = 1 and IHI = 0. <sparkle>✨</sparkle>
H is COMPLEX*16 array, dimension (LDH,N) — on entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper triangular matrix T from the Schur decomposition (the Schur form). If INFO = 0 and WANT is .FALSE., then the contents of H are unspecified on exit. (The output value of H when INFO > 0 is given under the description of INFO below.) This subroutine may explicitly set H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. <butterfly>🦋</butterfly>
LDH is INTEGER — the leading dimension of the array H. LDH >= max(1,N). 💖
W is COMPLEX*16 array, dimension (N) — the computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i). <sparkle>✨</sparkle>
ILOZ is INTEGER — that's it, that's the whole description 💅
IHIZ is INTEGER — specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. 🌸
Z is COMPLEX*16 array, dimension (LDZ,IHI) — if WANTZ is .FALSE., then Z is not referenced. If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the orthogonal Schur factor of H(ILO:IHI,ILO:IHI). (The output value of Z when INFO > 0 is given under the description of INFO below.) <butterfly>🦋</butterfly>
LDZ is INTEGER — the leading dimension of the array Z. if WANTZ is .TRUE. then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1. 💖
WORK is COMPLEX*16 array, dimension LWORK — on exit, if LWORK = -1, WORK(1) returns an estimate of the optimal value for LWORK. <sparkle>✨</sparkle>
LWORK is INTEGER — the dimension of the array WORK. LWORK >= max(1,N) is sufficient, but LWORK typically as large as 6*N may be required for optimal performance. A workspace query to determine the optimal workspace size is recommended. If LWORK = -1, then ZLAQR4 does a workspace query. In this case, ZLAQR4 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. 🔥
INFO is INTEGER
- = 0: successful exit
- > 0: if INFO = i, ZLAQR4 failed to compute all of the eigenvalues. Elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.)
  - If INFO > 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H.
  - If INFO > 0 and WANTT is .TRUE., then on exit

(*)  (initial value of H)*U  = U*(final value of H)

where U is a unitary matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI.

If INFO > 0 and WANTZ is .TRUE., then on exit

(final value of Z(ILO:IHI,ILOZ:IHIZ)
 =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U

where U is the unitary matrix in (*) (regardless of the value of WANTT.)

If INFO > 0 and WANTZ is .FALSE., then Z is not accessed.

✨ Author ✨

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

💖 Contributors: 💖

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

📚 References: 📚

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume 23, pages 929--947, 2002.
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part II: Aggressive Early Deflation, SIAM Journal of Matrix Analysis, volume 23, pages 948--973, 2002.

🌟 subroutine zlaqr5 ( ... ) 🌟

subroutine zlaqr5 (logical WANTT, logical WANTZ, integer KACC22, integer N,
    integer KTOP, integer KBOT, integer NSHFTS, complex*16, dimension( * )
    S, complex*16, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer
    IHIZ, complex*16, dimension( ldz, * ) Z, integer LDZ, complex*16,
    dimension( ldv, * ) V, integer LDV, complex*16, dimension( ldu, * ) U,
    integer LDU, integer NV, complex*16, dimension( ldwv, * ) WV, integer
    LDWV, integer NH, complex*16, dimension( ldwh, * ) WH, integer LDWH)

ZLAQR5 performs a single small-bulge multi-shift QR sweep. <sparkle>✨</sparkle>

🎯 Purpose:

ZLAQR5, called by ZLAQR0, performs a single small-bulge multi-shift QR sweep. <butterfly>🦋</butterfly>

✨ Parameters ✨

WANTT is LOGICAL — WANTT = .true. if the triangular Schur factor is being computed. WANTT is set to .false. otherwise. 💖
WANTZ is LOGICAL — WANTZ = .true. if the unitary Schur factor is being computed. WANTZ is set to .false. otherwise. 💖
KACC22 is INTEGER with value 0, 1, or 2. Specifies the computation mode of far-from-diagonal orthogonal updates.
- = 0: ZLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries.
- = 1: ZLAQR5 accumulates reflections and uses matrix-matrix multiply to update the far-from-diagonal matrix entries.
- = 2: Same as KACC22 = 1. This option used to enable exploiting the 2-by-2 structure during matrix multiplications, but this is no longer supported. ⚠️
N is INTEGER — N is the order of the Hessenberg matrix H upon which this subroutine operates. <sparkle>✨</sparkle>
KTOP is INTEGER — that's the whole thing 💅
KBOT is INTEGER — these are the first and last rows and columns of an isolated diagonal block upon which the QR sweep is to be applied. It is assumed without a check that either KTOP = 1 or H(KTOP,KTOP-1) = 0 and either KBOT = N or H(KBOT+1,KBOT) = 0. 🌸
NSHFTS is INTEGER — NSHFTS gives the number of simultaneous shifts. NSHFTS must be positive and even. periodt 💅
S is COMPLEX*16 array, dimension (NSHFTS) — S contains the shifts of origin that define the multi-shift QR sweep. On output S may be reordered. <butterfly>🦋</butterfly>
H is COMPLEX*16 array, dimension (LDH,N) — on input H contains a Hessenberg matrix. On output a multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied to the isolated diagonal block in rows and columns KTOP through KBOT. <sparkle>✨</sparkle>
LDH is INTEGER — LDH is the leading dimension of H just as declared in the calling procedure. LDH >= MAX(1,N). 💖
ILOZ is INTEGER — short and sweet 💅
IHIZ is INTEGER — specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N. 🌸
Z is COMPLEX*16 array, dimension (LDZ,IHIZ) — if WANTZ = .TRUE., then the QR Sweep unitary similarity transformation is accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ = .FALSE., then Z is unreferenced. <butterfly>🦋</butterfly>
LDZ is INTEGER — LDA is the leading dimension of Z just as declared in the calling procedure. LDZ >= N. 👀
V is COMPLEX*16 array, dimension (LDV,NSHFTS/2) — <sparkle>✨</sparkle>
LDV is INTEGER — LDV is the leading dimension of V as declared in the calling procedure. LDV >= 3. 💖
U is COMPLEX*16 array, dimension (LDU,2*NSHFTS) — <sparkle>✨</sparkle>
LDU is INTEGER — LDU is the leading dimension of U just as declared in the in the calling subroutine. LDU >= 2*NSHFTS. 🌸
NV is INTEGER — NV is the number of rows in WV agailable for workspace. NV >= 1. <butterfly>🦋</butterfly>
WV is COMPLEX*16 array, dimension (LDWV,2*NSHFTS) — <sparkle>✨</sparkle>
LDWV is INTEGER — LDWV is the leading dimension of WV as declared in the in the calling subroutine. LDWV >= NV. 💖
NH is INTEGER — NH is the number of columns in array WH available for workspace. NH >= 1. 🌸
WH is COMPLEX*16 array, dimension (LDWH,NH) — <sparkle>✨</sparkle>
LDWH is INTEGER — leading dimension of WH just as declared in the calling procedure. LDWH >= 2*NSHFTS. <butterfly>🦋</butterfly>

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🫶 Contributors:

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
Lars Karlsson, Daniel Kressner, and Bruno Lang
Thijs Steel, Department of Computer science, KU Leuven, Belgium

📚 References:

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume 23, pages 929--947, 2002.
Lars Karlsson, Daniel Kressner, and Bruno Lang, Optimally packed chains of bulges in multishift QR algorithms. ACM Trans. Math. Softw. 40, 2, Article 12 (February 2014).

🌸 subroutine `zlaqsb`

subroutine zlaqsb (character UPLO, integer N, integer KD, complex*16,
    dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * )
    S, double precision SCOND, double precision AMAX, character EQUED)

ZLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. no cap, that's the whole vibe check.

💡 Purpose:

ZLAQSB equilibrates a symmetric band matrix A using the scaling factors in the vector S. that's it, that's the job 💅

🎀 Parameters

UPLO — this one's a CHARACTER*1, and it tells us whether we're storing the upper or the lower triangular part of the symmetric matrix A.
- = 'U': Upper triangular
- = 'L': Lower triangular
N — INTEGER, bestie. the order of the matrix A. and like, N >= 0, obviously.
KD — INTEGER. this is the number of super-diagonals of A if UPLO = 'U', or the number of sub-diagonals if UPLO = 'L'. KD >= 0, always.
AB — COMPLEX*16 array, dimension (LDAB,N). okay so on entry, this holds the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. the j-th column of A lives in the j-th column of AB like this:
- if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
- if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).

on exit, if INFO = 0, you get the triangular factor U or L from the Cholesky factorization A = U**H *U or A = L*L**H of the band matrix A, in the same storage format as A. glow up 🦋

LDAB — INTEGER. the leading dimension of the array AB. LDAB >= KD+1. do not skimp.
S — DOUBLE PRECISION array, dimension (N). the scale factors for A ✨
SCOND — DOUBLE PRECISION. ratio of the smallest S(i) to the largest S(i). literally just a ratio.
AMAX — DOUBLE PRECISION. absolute value of the largest matrix entry. biggest one, periodt.
EQUED — CHARACTER*1. tells you whether or not equilibration happened:
- = 'N': No equilibration.
- = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).

🔍 Internal Parameters:

THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. if SCOND < THRESH, scaling is done.
LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. if AMAX > LARGE or AMAX < SMALL, scaling is done.

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🌈 subroutine `zlaqsp`

subroutine zlaqsp (character UPLO, integer N, complex*16, dimension( * ) AP,
    double precision, dimension( * ) S, double precision SCOND, double
    precision AMAX, character EQUED)

ZLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. packed girlies rise 💖

💡 Purpose:

ZLAQSP equilibrates a symmetric matrix A using the scaling factors in the vector S. same energy as above, different storage bestie.

🎀 Parameters

UPLO — CHARACTER*1. specifies whether the upper or lower triangular part of the symmetric matrix A is stored.
- = 'U': Upper triangular
- = 'L': Lower triangular
N — INTEGER. the order of the matrix A. N >= 0.
AP — COMPLEX*16 array, dimension (N*(N+1)/2). on entry, this is the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. the j-th column of A sits in the array AP like so:
- if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
- if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

on exit, you get the equilibrated matrix: diag(S) * A * diag(S), in the same storage format as A. she's transformed 💅

S — DOUBLE PRECISION array, dimension (N). the scale factors for A.
SCOND — DOUBLE PRECISION. ratio of the smallest S(i) to the largest S(i).
AMAX — DOUBLE PRECISION. absolute value of the largest matrix entry.
EQUED — CHARACTER*1. did we equilibrate or not, that is the question:
- = 'N': No equilibration.
- = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).

🔍 Internal Parameters:

THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. if SCOND < THRESH, scaling is done.
LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. if AMAX > LARGE or AMAX < SMALL, scaling is done.

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🦋 subroutine `zlar1v`

subroutine zlar1v (integer N, integer B1, integer BN, double precision
    LAMBDA, double precision, dimension( * ) D, double precision, dimension(
    * ) L, double precision, dimension( * ) LD, double precision, dimension(
    * ) LLD, double precision PIVMIN, double precision GAPTOL, complex*16,
    dimension( * ) Z, logical WANTNC, integer NEGCNT, double precision ZTZ,
    double precision MINGMA, integer R, integer, dimension( * ) ISUPPZ,
    double precision NRMINV, double precision RESID, double precision
    RQCORR, double precision, dimension( * ) WORK)

ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. this is the main character of eigenvector routines, fr fr. 🌟

💡 Purpose:

ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows B1 through BN of the tridiagonal matrix L D L**T - sigma I. when sigma is close to an eigenvalue, the computed vector is an accurate eigenvector. usually, r corresponds to the index where the eigenvector is largest in magnitude.

the following steps accomplish this computation:

Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
Computation of the diagonal elements of the inverse of L D L**T - sigma I by combining the above transforms, and choosing r as the index where the diagonal of the inverse is (one of the) largest in magnitude.
Computation of the (scaled) r-th column of the inverse using the twisted factorization obtained by combining the top part of the the stationary and the bottom part of the progressive transform.

obsessed with this pipeline honestly 🫶

🎀 Parameters

N — INTEGER. the order of the matrix L D L**T.
B1 — INTEGER. first index of the submatrix of L D L**T.
BN — INTEGER. last index of the submatrix of L D L**T.
LAMBDA — DOUBLE PRECISION. the shift. in order to compute an accurate eigenvector, LAMBDA should be a good approximation to an eigenvalue of L D L**T. so like, be accurate bestie.
L — DOUBLE PRECISION array, dimension (N-1). the (n-1) subdiagonal elements of the unit bidiagonal matrix L, in elements 1 to N-1.
D — DOUBLE PRECISION array, dimension (N). the n diagonal elements of the diagonal matrix D.
LD — DOUBLE PRECISION array, dimension (N-1). the n-1 elements L(i)*D(i).
LLD — DOUBLE PRECISION array, dimension (N-1). the n-1 elements L(i)*L(i)*D(i).
PIVMIN — DOUBLE PRECISION. the minimum pivot in the Sturm sequence.
GAPTOL — DOUBLE PRECISION. tolerance that indicates when eigenvector entries are negligible w.r.t. their contribution to the residual.
Z — COMPLEX*16 array, dimension (N). on input, all entries of Z must be set to 0 — not even joking, zero them all. on output, Z contains the (scaled) r-th column of the inverse. the scaling is such that Z(R) equals 1. 💅
WANTNC — LOGICAL. specifies whether NEGCNT has to be computed.
NEGCNT — INTEGER. if WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.
ZTZ — DOUBLE PRECISION. the square of the 2-norm of Z.
MINGMA — DOUBLE PRECISION. the reciprocal of the largest (in magnitude) diagonal element of the inverse of L D L**T - sigma I.
R — INTEGER. the twist index for the twisted factorization used to compute Z.
- on input, 0 <= R <= N. if R is input as 0, R is set to the index where (L D L**T - sigma I)^{-1} is largest in magnitude. if 1 <= R <= N, R is unchanged.
- on output, R contains the twist index used to compute Z. ideally, R designates the position of the maximum entry in the eigenvector.
ISUPPZ
ISUPPZ — INTEGER array, dimension (2) The support of the vector in Z, i.e., the vector Z is nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). that's the ✨zone✨ where stuff actually happens, bestie
NRMINV — DOUBLE PRECISION NRMINV = 1/SQRT( ZTZ )
RESID — DOUBLE PRECISION The residual of the FP vector. RESID = ABS( MINGMA )/SQRT( ZTZ )
RQCORR — DOUBLE PRECISION The Rayleigh Quotient correction to LAMBDA. RQCORR = MINGMA*TMP
WORK — DOUBLE PRECISION array, dimension (4*N)

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💖 Contributors:

Beresford Parlett, University of California, Berkeley, USA Jim Demmel, University of California, Berkeley, USA Inderjit Dhillon, University of Texas, Austin, USA Osni Marques, LBNL/NERSC, USA Christof Voemel, University of California, Berkeley, USA

🦋 `subroutine zlar2v`

subroutine zlar2v (integer N, complex*16, dimension( * ) X, complex*16,
    dimension( * ) Y, complex*16, dimension( * ) Z, integer INCX, double
    precision, dimension( * ) C, complex*16, dimension( * ) S, integer INCC)

ZLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices.

🎯 Purpose:

ZLAR2V applies a vector of complex plane rotations with real cosines
from both sides to a sequence of 2-by-2 complex Hermitian matrices,
defined by the elements of the vectors x, y and z. For i = 1,2,...,n

   (       x(i)  z(i) ) :=
   ( conjg(z(i)) y(i) )

     (  c(i) conjg(s(i)) ) (       x(i)  z(i) ) ( c(i) -conjg(s(i)) )
     ( -s(i)       c(i)  ) ( conjg(z(i)) y(i) ) ( s(i)        c(i)  )

💅 Parameters

N — INTEGER The number of plane rotations to be applied.
X — COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector x; the elements of x are assumed to be real.
Y — COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector y; the elements of y are assumed to be real.
Z — COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector z.
INCX — INTEGER The increment between elements of X, Y and Z. INCX > 0.
C — DOUBLE PRECISION array, dimension (1+(N-1)*INCC) The cosines of the plane rotations.
S — COMPLEX*16 array, dimension (1+(N-1)*INCC) The sines of the plane rotations.
INCC — INTEGER The increment between elements of C and S. INCC > 0.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🍓 `subroutine zlarcm`

subroutine zlarcm (integer M, integer N, double precision, dimension( lda, *
    ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB,
    complex*16, dimension( ldc, * ) C, integer LDC, double precision,
    dimension( * ) RWORK)

ZLARCM copies all or part of a real two-dimensional array to a complex array.

🎯 Purpose:

ZLARCM performs a very simple matrix-matrix multiplication:
         C := A * B,
where A is M by M and real; B is M by N and complex;
C is M by N and complex.

yes it's literally just C := A * B. no cap, it's giving "do the bare minimum and do it well" 💅

💅 Parameters

M — INTEGER The number of rows of the matrix A and of the matrix C. M >= 0.
N — INTEGER The number of columns and rows of the matrix B and the number of columns of the matrix C. N >= 0.
A — DOUBLE PRECISION array, dimension (LDA, M) On entry, A contains the M by M matrix A.
LDA — INTEGER The leading dimension of the array A. LDA >=max(1,M).
B — COMPLEX*16 array, dimension (LDB, N) On entry, B contains the M by N matrix B.
LDB — INTEGER The leading dimension of the array B. LDB >=max(1,M).
C — COMPLEX*16 array, dimension (LDC, N) On exit, C contains the M by N matrix C.
LDC — INTEGER The leading dimension of the array C. LDC >=max(1,M).
RWORK — DOUBLE PRECISION array, dimension (2*M*N)

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 `subroutine zlarf`

subroutine zlarf (character SIDE, integer M, integer N, complex*16,
    dimension( * ) V, integer INCV, complex*16 TAU, complex*16, dimension(
    ldc, * ) C, integer LDC, complex*16, dimension( * ) WORK)

ZLARF applies an elementary reflector to a general rectangular matrix.

🎯 Purpose:

ZLARF applies a complex elementary reflector H to a complex M-by-N
matrix C, from either the left or the right. H is represented in the
form

      H = I - tau * v * v**H

where tau is a complex scalar and v is a complex vector.

If tau = 0, then H is taken to be the unit matrix.

To apply H**H, supply conjg(tau) instead
tau.

💅 Parameters

SIDE — CHARACTER*1
- = 'L': form H * C
- = 'R': form C * H
M — INTEGER The number of rows of the matrix C.
N — INTEGER The number of columns of the matrix C.
V — COMPLEX*16 array, dimension

     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'

The vector v in the representation of H. V is not used if TAU = 0.

INCV — INTEGER The increment between elements of v. INCV <> 0.
TAU — COMPLEX*16 The value tau in the representation of H.
C — COMPLEX*16 array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by the matrix H * C if SIDE = 'L', or C * H if SIDE = 'R'.
LDC — INTEGER The leading dimension of the array C. LDC >= max(1,M).
WORK — COMPLEX*16 array, dimension

     (N) if SIDE = 'L'
  or (M) if SIDE = 'R'

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🧚 `subroutine zlarfb`

subroutine zlarfb (character SIDE, character TRANS, character DIRECT,
    character STOREV, integer M, integer N, integer K, complex*16,
    dimension( ldv, * ) V, integer LDV, complex*16, dimension( ldt, * ) T,
    integer LDT, complex*16, dimension( ldc, * ) C, integer LDC, complex*16,
    dimension( ldwork, * ) WORK, integer LDWORK)

ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.

🎯 Purpose:

ZLARFB applies a complex block reflector H or its transpose H**H to a
complex M-by-N matrix C, from either the left or the right.

💅 Parameters

SIDE — CHARACTER*1
- = 'L': apply H or H**H from the Left
- = 'R': apply H or H**H from the Right
TRANS — CHARACTER*1
- = 'N': apply H (No transpose)
- = 'C': apply H**H (Conjugate transpose)
DIRECT — CHARACTER*1 Indicates how H is formed from a product of elementary reflectors
- = 'F': H = H(1) H(2) . . . H(k) (Forward)
- = 'B': H = H(k) . . . H(2) H(1) (Backward)
STOREV — CHARACTER*1 Indicates how the vectors which define the elementary reflectors are stored:
- = 'C': Columnwise
- = 'R': Rowwise
M — INTEGER The number of rows of the matrix C.
N — INTEGER The number of columns of the matrix C.
K — INTEGER The order of the matrix T (= the number of elementary reflectors whose product defines the block reflector). If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0.
V — COMPLEX*16 array, dimension

                     (LDV,K) if STOREV = 'C'
                     (LDV,M) if STOREV = 'R' and SIDE = 'L'
                     (LDV,N) if STOREV = 'R' and SIDE = 'R'

See Further Details.

LDV — INTEGER The leading dimension of the array V. If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
T — COMPLEX*16 array, dimension (LDT,K) The triangular K-by-K matrix T in the representation of the block reflector.
LDT — INTEGER The leading dimension of the array T. LDT >= K.
C — COMPLEX*16 array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
LDC — INTEGER The leading dimension of the array C. LDC >= max(1,M).
WORK — COMPLEX*16 array, dimension (LDWORK,K)
LDWORK — INTEGER The leading dimension of the array WORK. If SIDE = 'L', LDWORK >= max(1,N); if SIDE = 'R', LDWORK >= max(1,M).

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🔎 Further Details:

The shape of the matrix V and the storage of the vectors which define
the H(i) is best illustrated by the following example with n = 5 and
k = 3. The elements equal to 1 are not stored; the corresponding
array elements are modified but restored on exit. The rest of the
array is not used.

and the actual diagram of all that? bestie it's right around the corner… 👀

…so bestie, picking up right where we left off — here's literally what V looks like in each of the four DIRECT / STOREV flavours 💅:

             DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':

                          V = (  1       )                 V = (  1 v1 v1 v1 v1 )
                              ( v1  1    )                     (     1 v2 v2 v2 )
                              ( v1 v2  1 )                     (        1 v3 v3 )
                              ( v1 v2 v3 )
                              ( v1 v2 v3 )

             DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':

                          V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
                              ( v1 v2 v3 )                     ( v2 v2 v2  1    )
                              (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
                              (     1 v3 )
                              (        1 )

🦋 `zlarfb_gett` — the subroutine

subroutine zlarfb_gett (character IDENT, integer M, integer N, integer K,
    complex*16, dimension( ldt, * ) T, integer LDT, complex*16, dimension(
    lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB,
    complex*16, dimension( ldwork, * ) WORK, integer LDWORK)

ZLARFB_GETT

✨ Purpose

Okay so ZLARFB_GETT takes a complex Householder block reflector H and applies it from the left to a complex (K+M)-by-N 'triangular-pentagonal' matrix, which is genuinely just two block matrices glued together: an upper trapezoidal K-by-N matrix A sitting in the array A, and a rectangular M-by-(N-K) matrix B chilling in the array B. The block reflector H lives in a compact WY-representation, where the elementary reflectors are hanging out in the arrays A, B and T. Scroll down to the Further Details section for the full tea ☕.

🎀 Parameters

IDENT — IDENT is CHARACTER*1. If IDENT = not 'I', or not 'i', then V1 is unit lower-triangular and stored in the left K-by-K block of the input matrix A. If IDENT = 'I' or 'i', then V1 is an identity matrix and not stored. Peep the Further Details section for the whole picture 💖.
M — M is INTEGER. The number of rows of the matrix B. M >= 0.
N — N is INTEGER. The number of columns of the matrices A and B. N >= 0.
K — K is INTEGER. The number or rows of the matrix A. K is also order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. 0 <= K <= N.
T — T is COMPLEX*16 array, dimension (LDT,K). The upper-triangular K-by-K matrix T in the representation of the block reflector.
LDT — LDT is INTEGER. The leading dimension of the array T. LDT >= K.
A — A is COMPLEX*16 array, dimension (LDA,N).
- On entry:
  1. In the K-by-N upper-trapezoidal part A: input matrix A.
  2. In the columns below the diagonal: columns of V1 (ones are not stored on the diagonal).
- On exit: A is overwritten by rectangular K-by-N product H*A.
- See Further Details section 🫶.
LDA — LDB is INTEGER. The leading dimension of the array A. LDA >= max(1,K).

bestie yes, it literally says LDB in the source, we're not touching it, it's canon 💀

B — B is COMPLEX*16 array, dimension (LDB,N).
- On entry:
  1. In the M-by-(N-K) right block: input matrix B.
  2. In the M-by-N left block: columns of V2.
- On exit: B is overwritten by rectangular M-by-N product H*B.
- See Further Details section 🌸.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M).
WORK — WORK is COMPLEX*16 array, dimension (LDWORK,max(K,N-K)).
LDWORK — LDWORK is INTEGER. The leading dimension of the array WORK. LDWORK>=max(1,K).

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💜 Contributors:

November 2020, Igor Kozachenko, Computer Science Division, University of California, Berkeley

🌈 Further Details:

(1) Description of the Algebraic Operation.

The matrix A is a K-by-N matrix composed of two column block matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):

A = ( A1, A2 ).

The matrix B is an M-by-N matrix composed of two column block matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):

B = ( B1, B2 ).

Perform the operation:

   ( A_out ) := H * ( A_in ) = ( I - V * T * V**H ) * ( A_in ) =
   ( B_out )        ( B_in )                          ( B_in )
              = ( I - ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in )
                      ( V2 )                            ( B_in )

On input:

   ( A_in )  consists of two block columns:
   ( B_in )

   ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))
   ( B_in )   (( B1_in ) ( B2_in ))   ((     0 ) ( B2_in )),

where the column blocks are:

   (  A1_in )  is a K-by-K upper-triangular matrix stored in the
               upper triangular part of the array A(1:K,1:K).
   (  B1_in )  is an M-by-K rectangular ZERO matrix and not stored.

   ( A2_in )  is a K-by-(N-K) rectangular matrix stored
              in the array A(1:K,K+1:N).
   ( B2_in )  is an M-by-(N-K) rectangular matrix stored
              in the array B(1:M,K+1:N).

   V = ( V1 )
       ( V2 )

where:

if IDENT == 'I', V1 is a K-by-K identity matrix, not stored;
if IDENT != 'I', V1 is a K-by-K unit lower-triangular matrix, stored in the lower-triangular part of the array A(1:K,1:K) (ones are not stored),

and V2 is an M-by-K rectangular stored the array B(1:M,1:K), (because on input B1_in is a rectangular zero matrix that is not stored and the space is used to store V2).

c) T is a K-by-K upper-triangular matrix stored in the array T(1:K,1:K).

On output:

   ( A_out ) consists of two  block columns:
   ( B_out )

   ( A_out ) = (( A1_out ) ( A2_out ))
   ( B_out )   (( B1_out ) ( B2_out )),

where the column blocks are:

   ( A1_out )  is a K-by-K square matrix, or a K-by-K
               upper-triangular matrix, if V1 is an
               identity matrix. AiOut is stored in
               the array A(1:K,1:K).
   ( B1_out )  is an M-by-K rectangular matrix stored
               in the array B(1:M,K:N).

   ( A2_out )  is a K-by-(N-K) rectangular matrix stored
               in the array A(1:K,K+1:N).
   ( B2_out )  is an M-by-(N-K) rectangular matrix stored
               in the array B(1:M,K+1:N).

The operation above can be represented as the same operation on each block column:

   ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**H ) * ( A1_in )
   ( B1_out )        (     0 )                          (     0 )

   ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**H ) * ( A2_in )
   ( B2_out )        ( B2_in )                          ( B2_in )

If IDENT != 'I':

The computation for column block 1:

   A1_out: = A1_in - V1*T*(V1**H)*A1_in

   B1_out: = - V2*T*(V1**H)*A1_in

The computation for column block 2, which exists if N > K:

   A2_out: = A2_in - V1*T*( (V1**H)*A2_in + (V2**H)*B2_in )

   B2_out: = B2_in - V2*T*( (V1**H)*A2_in + (V2**H)*B2_in )

If IDENT == 'I':

The operation for column block 1:

   A1_out: = A1_in - V1*T*A1_in

   B1_out: = - V2*T*A1_in

The computation for column block 2, which exists if N > K:

   A2_out: = A2_in - T*( A2_in + (V2**H)*B2_in )

   B2_out: = B2_in - V2*T*( A2_in + (V2**H)*B2_in )

(2) Description of the Algorithmic Computation.

In the first step, we compute column block 2, i.e. A2 and B2. Here, we need to use the K-by-(N-K) rectangular workspace matrix W2 that is of the same size as the matrix A2. W2 is stored in the array WORK(1:K,1:(N-K)).

In the second step, we compute column block 1, i.e. A1 and B1. Here, we need to use the K-by-K square workspace matrix W1 that is of the same size as the as the matrix A1. W1 is stored in the array WORK(1:K,1:K).

NOTE: Hence, in this routine, we need the workspace array WORK only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from the first step and W1 from the second step. Iconic efficiency honestly ✨

Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I', more computations than in the Case (B).

   if( IDENT != 'I' ) then
    if ( N > K ) then
      (First Step - column block 2)
      col2_(1) W2: = A2
      col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2
      col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
      col2_(4) W2: = T * W2
      col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
      col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
      col2_(7) A2: = A2 - W2
    else
      (Second Step - column block 1)
      col1_(1) W1: = A1
      col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1
      col1_(3) W1: = T * W1
      col1_(4) B1: = - V2 * W1 = - B1 * W1
      col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
      col1_(6) square A1: = A1 - W1
    end if
   end if

Case (B), when V1 is an identity matrix, i.e. IDENT == 'I', less computations than in the Case (A)

               if( IDENT == 'I' ) then
                if ( N > K ) then
                  (First Step - column block 2)
                  col2_(1) W2: = A2
                  col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
                  col2_(4) W2: = T * W2
                  col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
                  col2_(7) A2: = A2 - W2
                else
                  (Second Step - column block 1)
                  col1_(1) W1: = A1
                  col1_(3) W1: = T * W1
                  col1_(4) B1: = - V2 * W1 = - B1 * W1
                  col1_(6) upper-triangular_of_(A1): = A1 - W1
                end if
               end if

               Combine these cases (A) and (B) together, this is the resulting
               algorithm:

               if ( N > K ) then

                 (First Step - column block 2)

                 col2_(1)  W2: = A2
                 if( IDENT != 'I' ) then
                   col2_(2)  W2: = (V1**H) * W2
                                 = (unit_lower_tr_of_(A1)**H) * W2
                 end if
                 col2_(3)  W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2]
                 col2_(4)  W2: = T * W2
                 col2_(5)  B2: = B2 - V2 * W2 = B2 - B1 * W2
                 if( IDENT != 'I' ) then
                   col2_(6)    W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
                 end if
                 col2_(7) A2: = A2 - W2

               else

               (Second Step - column block 1)

                 col1_(1) W1: = A1
                 if( IDENT != 'I' ) then
                   col1_(2) W1: = (V1**H) * W1
                               = (unit_lower_tr_of_(A1)**H) * W1
                 end if
                 col1_(3) W1: = T * W1
                 col1_(4) B1: = - V2 * W1 = - B1 * W1
                 if( IDENT != 'I' ) then
                   col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
                   col1_(6_a) below_diag_of_(A1): =  - below_diag_of_(W1)
                 end if
                 col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)

               end if

`zlarfg` ✨

subroutine zlarfg (integer N, complex*16 ALPHA, complex*16, dimension( * ) X, integer INCX, complex*16 TAU)

ZLARFG generates an elementary reflector (Householder matrix). 💅

Purpose:

bestie, this one builds a complex elementary reflector H of order n, such that

      H**H * ( alpha ) = ( beta ),   H**H * H = I.
             (   x   )   (   0  )

where alpha and beta are scalars, with beta real, and x is an (n-1)-element complex vector. H is represented in the form

      H = I - tau * ( 1 ) * ( 1 v**H ) ,
                    ( v )

where tau is a complex scalar and v is a complex (n-1)-element vector. Note that H is not hermitian. 🙈

If the elements of x are all zero and alpha is real, then tau = 0 and H is taken to be the unit matrix. no notes, iconic behaviour ✨

Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1. 🫶

Parameters

N — N is INTEGER. The order of the elementary reflector.
ALPHA — ALPHA is COMPLEX*16. On entry, the value alpha. On exit, it is overwritten with the value beta.
X — X is COMPLEX*16 array, dimension (1+(N-2)*abs(INCX)). On entry, the vector x. On exit, it is overwritten with the vector v.
INCX — INCX is INTEGER. The increment between elements of X. INCX > 0.
TAU — TAU is COMPLEX*16. The value tau.

Author 🎀

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

`zlarfgp` 🌸

subroutine zlarfgp (integer N, complex*16 ALPHA, complex*16, dimension( * ) X, integer INCX, complex*16 TAU)

ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. 👑

Purpose:

same energy as zlarfg but the beta comes out non-negative, periodt. ZLARFGP generates a complex elementary reflector H of order n, such that

      H**H * ( alpha ) = ( beta ),   H**H * H = I.
             (   x   )   (   0  )

where alpha and beta are scalars, beta is real and non-negative, and x is an (n-1)-element complex vector. H is represented in the form

      H = I - tau * ( 1 ) * ( 1 v**H ) ,
                    ( v )

where tau is a complex scalar and v is a complex (n-1)-element vector. Note that H is not hermitian. 😭

If the elements of x are all zero and alpha is real, then tau = 0 and H is taken to be the unit matrix. 💖

Parameters

N — N is INTEGER. The order of the elementary reflector.
ALPHA — ALPHA is COMPLEX*16. On entry, the value alpha. On exit, it is overwritten with the value beta.
X — X is COMPLEX*16 array, dimension (1+(N-2)*abs(INCX)). On entry, the vector x. On exit, it is overwritten with the vector v.
INCX — INCX is INTEGER. The increment between elements of X. INCX > 0.
TAU — TAU is COMPLEX*16. The value tau.

Author 🎀

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

`zlarft` 🦋

subroutine zlarft (character DIRECT, character STOREV, integer N, integer K, complex*16, dimension( ldv, * ) V, integer LDV, complex*16, dimension( * ) TAU, complex*16, dimension( ldt, * ) T, integer LDT)

ZLARFT forms the triangular factor T of a block reflector H = I - vtvH ⭐

Purpose:

ZLARFT forms the triangular factor T of a complex block reflector H of order n, which is defined as a product of k elementary reflectors. 🧚

If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.

If STOREV = 'C', the vector which defines the elementary reflector H(i) is stored in the i-th column of the array V, and

   H  =  I - V * T * V**H

If STOREV = 'R', the vector which defines the elementary reflector H(i) is stored in the i-th row of the array V, and

   H  =  I - V**H * T * V

Parameters

DIRECT — DIRECT is CHARACTER*1. Specifies the order in which the elementary reflectors are multiplied to form the block reflector:
- = 'F': H = H(1) H(2) . . . H(k) (Forward)
- = 'B': H = H(k) . . . H(2) H(1) (Backward)
STOREV — STOREV is CHARACTER*1. Specifies how the vectors which define the elementary reflectors are stored (see also Further Details):
- = 'C': columnwise
- = 'R': rowwise
N — N is INTEGER. The order of the block reflector H. N >= 0.
K — K is INTEGER. The order of the triangular factor T (= the number of elementary reflectors). K >= 1.
V — V is COMPLEX*16 array, dimension (LDV,K) if STOREV = 'C', (LDV,N) if STOREV = 'R'. The matrix V. See further details.
LDV — LDV is INTEGER. The leading dimension of the array V. If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
TAU — TAU is COMPLEX*16 array, dimension (K). TAU(i) must contain the scalar factor of the elementary reflector H(i).
T — T is COMPLEX*16 array, dimension (LDT,K). The k by k triangular factor T of the block reflector. If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is lower triangular. The rest of the array is not used.
LDT — LDT is INTEGER. The leading dimension of the array T. LDT >= K.

Author 🎀

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

The shape of the matrix V and the storage of the vectors which define the H(i) is best illustrated by the following example with n = 5 and k = 3. The elements equal to 1 are not stored. 🔥

DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':

             V = (  1       )                 V = (  1 v1 v1 v1 v1 )
                 ( v1  1    )                     (     1 v2 v2 v2 )
                 ( v1 v2  1 )                     (        1 v3 v3 )
                 ( v1 v2 v3 )
                 ( v1 v2 v3 )

DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':

             V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
                 ( v1 v2 v3 )                     ( v2 v2 v2  1    )
                 (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
                 (     1 v3 )
                 (        1 )

`zlarfx` 💅

subroutine zlarfx (character SIDE, integer M, integer N, complex*16, dimension( * ) V, complex*16 TAU, complex*16, dimension( ldc, * ) C, integer LDC, complex*16, dimension( * ) WORK)

ZLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10. 😤

Purpose:

ZLARFX applies a complex elementary reflector H to a complex m by n matrix C, from either the left or the right. H is represented in the form

      H = I - tau * v * v**H

where tau is a complex scalar and v is a complex vector.

If tau = 0, then H is taken to be the unit matrix 🫶

This version uses inline code if H has order < 11. ⚡

Parameters

SIDE — SIDE is CHARACTER*1.
- = 'L': form H * C
- = 'R': form C * H
M — M is INTEGER. The number of rows of the matrix C.
N — N is INTEGER. The number of columns of the matrix C.
V — V is COMPLEX*16 array, dimension (M) if SIDE = 'L' or (N) if SIDE = 'R'. The vector v in the representation of H.
TAU — TAU is COMPLEX*16. The value tau in the representation of H.
C — C is COMPLEX*16 array, dimension (LDC,N). On entry, the m by n matrix C. On exit, C is overwritten by the matrix H * C if SIDE = 'L', or C * H if SIDE = 'R'.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
WORK

WORK is COMPLEX*16 array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R'. And bestie, WORK is not referenced if H has order < 11.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ subroutine zlarfy ✨

   subroutine zlarfy (character UPLO, integer N, complex*16, dimension( * ) V,
       integer INCV, complex*16 TAU, complex*16, dimension( ldc, * ) C, integer
       LDC, complex*16, dimension( * ) WORK)

ZLARFY

Purpose:

Okay so ZLARFY is giving main character energy: it applies an elementary reflector (aka Householder matrix) H to an n x n Hermitian matrix C, and it does it from both the left and the right. Slay.

H is represented in the form

H = I - tau * v * v'

where tau is a scalar and v is a vector.

If tau is zero, then H is taken to be the unit matrix. Periodt.

Parameters

UPLO — UPLO is a CHARACTER*1. It tells us whether we're storing the upper or lower triangular part of that Hermitian matrix C. = 'U': upper triangle, = 'L': lower triangle. No in-between, bestie.
N — N is an INTEGER. The number of rows and columns of the matrix C. N >= 0.
V — V is a COMPLEX*16 array, dimension (1 + (N-1)*abs(INCV)). The vector v as described above. <sparkles>✨</sparkles>
INCV — INCV is an INTEGER. The increment between successive elements of v. INCV must not be zero. Do not let it be zero, fr fr.
TAU — TAU is a COMPLEX*16. The value tau as described above.
C — C is a COMPLEX*16 array, dimension (LDC, N). On entry, the matrix C. On exit, C is overwritten by H * C * H'. She's changed.
LDC — LDC is an INTEGER. The leading dimension of the array C. LDC >= max( 1, N ).
WORK — WORK is a COMPLEX*16 array, dimension (N).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ subroutine zlargv ✨

   subroutine zlargv (integer N, complex*16, dimension( * ) X, integer INCX,
       complex*16, dimension( * ) Y, integer INCY, double precision, dimension(
       * ) C, integer INCC)

ZLARGV generates a vector of plane rotations with real cosines and complex sines. iykyk

Purpose:

ZLARGV generates a vector of complex plane rotations with real cosines, determined by elements of the complex vectors x and y. For i = 1,2,...,n

(        c(i)   s(i) ) ( x(i) ) = ( r(i) )
( -conjg(s(i))  c(i) ) ( y(i) ) = (   0  )

where c(i)**2 + ABS(s(i))**2 = 1

The following conventions are used (these are the same as in ZLARTG, but differ from the BLAS1 routine ZROTG):

If y(i)=0, then c(i)=1 and s(i)=0.
If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.

Parameters

N — N is an INTEGER. The number of plane rotations to be generated.
X — X is a COMPLEX*16 array, dimension (1+(N-1)*INCX). On entry, the vector x. On exit, x(i) is overwritten by r(i), for i = 1,...,n. Plot twist.
INCX — INCX is an INTEGER. The increment between elements of X. INCX > 0.
Y — Y is a COMPLEX*16 array, dimension (1+(N-1)*INCY). On entry, the vector y. On exit, the sines of the plane rotations.
INCY — INCY is an INTEGER. The increment between elements of Y. INCY > 0.
C — C is a DOUBLE PRECISION array, dimension (1+(N-1)*INCC). The cosines of the plane rotations.
INCC — INCC is an INTEGER. The increment between elements of C. INCC > 0.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel

This version has a few statements commented out for thread safety (machine parameters are computed on each entry). 10 feb 03, SJH.

✨ subroutine zlarnv ✨

   subroutine zlarnv (integer IDIST, integer, dimension( 4 ) ISEED, integer N,
       complex*16, dimension( * ) X)

ZLARNV returns a vector of random numbers from a uniform or normal distribution. Random girlypop behavior.

Purpose:

ZLARNV returns a vector of n random complex numbers from a uniform or normal distribution.

Parameters

IDIST — IDIST is an INTEGER. Specifies the distribution of the random numbers:
- = 1: real and imaginary parts each uniform (0,1)
- = 2: real and imaginary parts each uniform (-1,1)
- = 3: real and imaginary parts each normal (0,1)
- = 4: uniformly distributed on the disc abs(z) < 1
- = 5: uniformly distributed on the circle abs(z) = 1
ISEED — ISEED is an INTEGER array, dimension (4). On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. On exit, the seed is updated. Respect the seed rules, no cap.
N — N is an INTEGER. The number of random numbers to be generated.
X — X is a COMPLEX*16 array, dimension (N). The generated random numbers.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

This routine calls the auxiliary routine DLARUV to generate random real numbers from a uniform (0,1) distribution, in batches of up to 128 using vectorisable code. The Box-Muller method is used to transform numbers from a uniform to a normal distribution.

✨ subroutine zlarrv ✨

   subroutine zlarrv (integer N, double precision VL, double precision VU,
       double precision, dimension( * ) D, double precision, dimension( * ) L,
       double precision PIVMIN, integer, dimension( * ) ISPLIT, integer M,
       integer DOL, integer DOU, double precision MINRGP, double precision
       RTOL1, double precision RTOL2, double precision, dimension( * ) W,
       double precision, dimension( * ) WERR, double precision, dimension( * )
       WGAP, integer, dimension( * ) IBLOCK, integer, dimension( * ) INDEXW,
       double precision, dimension( * ) GERS, complex*16, dimension( ldz, * )
       Z, integer LDZ, integer, dimension( * ) ISUPPZ, double precision,
       dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.

Purpose:

ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. The input eigenvalues should have been computed by DLARRE. So make sure you ran that first, bestie.

Parameters

N — N is an INTEGER. The order of the matrix. N >= 0.
VL — VL is a DOUBLE PRECISION. Lower bound of the interval that contains the desired eigenvalues. VL < VU. Needed to compute gaps on the left or right end of the extremal eigenvalues in the desired RANGE.
VU — VU is a DOUBLE PRECISION. Upper bound of the interval that contains the desired eigenvalues. VL < VU. Needed to compute gaps on the left or right end of the extremal eigenvalues in the desired RANGE.
D — D is a DOUBLE PRECISION array, dimension (N). On entry, the N diagonal elements of the diagonal matrix D. On exit, D may be overwritten. She may not survive.
L — L is a DOUBLE PRECISION array, dimension (N). On entry, the (N-1) subdiagonal elements of the unit bidiagonal matrix L are in elements 1 to N-1 of L (if the matrix is not split.) At the end of each block is stored the corresponding shift as given by DLARRE. On exit, L is overwritten.
PIVMIN — PIVMIN is a DOUBLE PRECISION. The minimum pivot allowed in the Sturm sequence.
ISPLIT — ISPLIT is an INTEGER array, dimension (N). The splitting points, at which T breaks up into blocks. The first block consists of rows/columns 1 to ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 through ISPLIT( 2 ), etc.
M — M is an INTEGER. The total number of input eigenvalues. 0 <= M <= N.
DOL — DOL is an INTEGER.
DOU — DOU is an INTEGER. If the user wants to compute only selected eigenvectors from all the eigenvalues supplied, he can specify an index range DOL:DOU. Or else the setting DOL=1, DOU=M should be applied. Note that DOL and DOU refer to the order in which the eigenvalues are stored in W. If the user wants to compute only selected eigenpairs, then the columns DOL-1 to DOU+1 of the eigenvector space Z contain the computed eigenvectors. All other columns of Z are set to zero.
MINRGP — MINRGP is a DOUBLE PRECISION.
RTOL1 — RTOL1 is a DOUBLE PRECISION.
RTOL2 — RTOL2 is a DOUBLE PRECISION. Parameters for bisection. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ).
W — W is a DOUBLE PRECISION array, dimension (N). The first M elements of W contain the APPROXIMATE eigenvalues for which eigenvectors are to be computed. The eigenvalues should be grouped by split-off block and ordered from smallest to largest within the block ( The output array W from DLARRE is expected here ). Furthermore, they are with respect to the shift of the corresponding root representation for their block. On exit, W holds the eigenvalues of the UNshifted matrix.
WERR — WERR is a DOUBLE PRECISION array, dimension (N). The first M elements contain the semiwidth of the uncertainty interval of the corresponding eigenvalue in W
WGAP — WGAP is a DOUBLE PRECISION array, dimension (N). The separation from the right neighbor eigenvalue in W.
IBLOCK —

IBLOCK

IBLOCK is an INTEGER array, dimension (N). It's the block (submatrix) index for each eigenvalue in W — like IBLOCK(i)=1 if eigenvalue W(i) belongs to the first block from the top, =2 if W(i) belongs to the second block, etc. Block party energy, but make it numerical. 💅

INDEXW

INDEXW is an INTEGER array, dimension (N). These are the indices of the eigenvalues within each block (submatrix). So if INDEXW(i)= 10 and IBLOCK(i)=2, that means the i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. Precision is the vibe. ✨

GERS

GERS is a DOUBLE PRECISION array, dimension (2*N). It holds the N Gerschgorin intervals, where the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i)). Important: these intervals should be computed from the original UNshifted matrix. Do not shift, bestie. 🚫

Z

Z is a COMPLEX*16 array, dimension (LDZ, max(1,M) ). If INFO = 0, the first M columns of Z contain the orthonormal eigenvectors of the matrix T corresponding to the input eigenvalues, with the i-th column of Z holding the eigenvector associated with W(i). Note: you must ensure that at least max(1,M) columns are supplied in the array Z. No cap. 👑

LDZ

LDZ is INTEGER. It's the leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', then LDZ >= max(1,N). Rules are rules.

ISUPPZ

ISUPPZ is an INTEGER array, dimension ( 2*max(1,M) ). It's 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 ). Supportive queen. 🫶

WORK

WORK is a DOUBLE PRECISION array, dimension (12*N). That's a lot of workspace, but she's worth it.

IWORK

IWORK is an INTEGER array, dimension (7*N). Integer energy.

INFO

INFO is INTEGER. Here's the tea:

= 0: successful exit ✨
> 0: a problem occurred in ZLARRV. Uh oh.
< 0: one of the called subroutines signaled an internal problem. Needs inspection of the corresponding parameter IINFO for further information.
=-1: problem in DLARRB when refining a child's eigenvalues.
=-2: problem in DLARRF when computing the RRR of a child. When a child is inside a tight cluster, it can be difficult to find an RRR. A partial remedy from the user's point of view is to make the parameter MINRGP smaller and recompile. However, as the orthogonality of the computed vectors is proportional to 1/MINRGP, the user should be aware that he might be trading in precision when he decreases MINRGP. Trade-offs, bestie.
=-3: problem in DLARRB when refining a single eigenvalue after the Rayleigh correction was rejected.
= 5: the Rayleigh Quotient Iteration failed to converge to full accuracy in MAXITR steps. Not slay.

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Contributors:

Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

subroutine `zlartv`

subroutine zlartv (integer N, complex*16, dimension( * ) X, integer INCX,
     complex*16, dimension( * ) Y, integer INCY, double precision, dimension(
     * ) C, complex*16, dimension( * ) S, integer INCC)

ZLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a pair of vectors. Rotation nation. 🌀

Purpose:

ZLARTV applies a vector of complex plane rotations with real cosines to elements of the complex vectors x and y. For i = 1,2,...,n:

   ( x(i) ) := (        c(i)   s(i) ) ( x(i) )
   ( y(i) )    ( -conjg(s(i))  c(i) ) ( y(i) )

Parameters

N — N is INTEGER. The number of plane rotations to be applied.
X — X is COMPLEX*16 array, dimension (1+(N-1)*INCX). The vector x.
INCX — INCX is INTEGER. The increment between elements of X. INCX > 0.
Y — Y is COMPLEX*16 array, dimension (1+(N-1)*INCY). The vector y.
INCY — INCY is INTEGER. The increment between elements of Y. INCY > 0.
C — C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC). The cosines of the plane rotations.
S — S is COMPLEX*16 array, dimension (1+(N-1)*INCC). The sines of the plane rotations.
INCC — INCC is INTEGER. The increment between elements of C and S. INCC > 0.

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

subroutine `zlascl`

subroutine zlascl (character TYPE, integer KL, integer KU, double precision
     CFROM, double precision CTO, integer M, integer N, complex*16,
     dimension( lda, * ) A, integer LDA, integer INFO)

ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. Scaling realness.

Purpose:

ZLASCL multiplies the M by N complex matrix A by the real scalar CTO/CFROM. This is done without over/underflow as long as the final result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that A may be full, upper triangular, lower triangular, upper Hessenberg, or banded.

Parameters

TYPE — TYPE is CHARACTER*1. TYPE indices the storage type of the input matrix.
- = 'G': A is a full matrix.
- = 'L': A is a lower triangular matrix.
- = 'U': A is an upper triangular matrix.
- = 'H': A is an upper Hessenberg matrix.
- = 'B': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the lower half stored.
- = 'Q': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the upper half stored.
- = 'Z': A is a band matrix with lower bandwidth KL and upper bandwidth KU. See ZGBTRF for storage details.
KL — KL is INTEGER. The lower bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'.
KU — KU is INTEGER. The upper bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'.
CFROM — CFROM is DOUBLE PRECISION.
CTO — CTO is DOUBLE PRECISION. The matrix A is multiplied by CTO/CFROM. A(I,J) is computed without over/underflow if the final result CTO*A(I,J)/CFROM can be represented without over/underflow. CFROM must be nonzero.
M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
A — A is COMPLEX*16 array, dimension (LDA,N). The matrix to be multiplied by CTO/CFROM. See TYPE for the storage type.
LDA — LDA is INTEGER. The leading dimension of the array A. If TYPE = 'G', 'L', 'U', 'H', LDA >= max(1,M); TYPE = 'B', LDA >= KL+1; TYPE = 'Q', LDA >= KU+1; TYPE = 'Z', LDA >= 2*KL+KU+1.
INFO — INFO is INTEGER.
- 0 - successful exit
- <0 - if INFO = -i, the i-th argument had an illegal value.

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

subroutine `zlaset`

subroutine zlaset (character UPLO, integer M, integer N, complex*16 ALPHA,
     complex*16 BETA, complex*16, dimension( lda, * ) A, integer LDA)

ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values. Set it off.

Purpose:

ZLASET initializes a 2-D array A to BETA on the diagonal and ALPHA on the offdiagonals.

Parameters

UPLO — UPLO is CHARACTER*1. Specifies the part of the matrix A to be set.
- = 'U': Upper triangular part is set. The lower triangle is unchanged.
- = 'L': Lower triangular part is set. The upper triangle is unchanged.
- Otherwise: All of the matrix A is set.
M — M is INTEGER. On entry, M specifies the number of rows of A.
N — N is INTEGER. On entry, N specifies the number of columns of A.
ALPHA — ALPHA is COMPLEX*16. All the offdiagonal array elements are set to ALPHA.
BETA — BETA is COMPLEX*16. All the diagonal array elements are set to BETA.
A — A is COMPLEX*16 array, dimension (LDA,N). On entry, the m by n matrix A. On exit, A(i,j) = ALPHA, 1 <= i <= m, 1 <= j <= n, i.ne.j; A(i,i) = BETA , 1 <= i <= min(m,n)
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).

Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

subroutine `zlasr`

subroutine zlasr (character SIDE, character PIVOT, character DIRECT, integer
     M, integer N, double precision, dimension( * ) C, double precision,
     dimension( * ) S, complex*16, dimension( lda, * ) A, integer LDA)

ZLASR applies a sequence of plane rotations to a general rectangular matrix. Rotation sequence, but make it fashion. 🌈

Purpose:

ZLASR applies a sequence of real plane rotations to a complex matrix A, from either the left or the right.

When SIDE = 'L', the transformation takes the form

A := P*A

and when SIDE = 'R', the transformation takes the form

A := A*P**T

where P is an orthogonal matrix made of a whole sequence of z plane rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', and P**T is just the transpose of P. slay.

When DIRECT = 'F' (Forward sequence), then

P = P(z-1) * ... * P(2) * P(1)

and when DIRECT = 'B' (Backward sequence), then

P = P(1) * P(2) * ... * P(z-1)

where P(k) is a plane rotation matrix defined by the 2-by-2 rotation

R(k) = (  c(k)  s(k) )
     = ( -s(k)  c(k) ).

When PIVOT = 'V' (Variable pivot), the rotation is performed for the plane (k,k+1), i.e., P(k) has the form

P(k) = (  1                                            )
       (       ...                                     )
       (              1                                )
       (                   c(k)  s(k)                  )
       (                  -s(k)  c(k)                  )
       (                                1              )
       (                                     ...       )
       (                                            1  )

where R(k) appears as a rank-2 modification to the identity matrix in rows and columns k and k+1. it's giving "small little tweak, big consequences" 💅

When PIVOT = 'T' (Top pivot), the rotation is performed for the plane (1,k+1), so P(k) has the form

P(k) = (  c(k)                    s(k)                 )
       (         1                                     )
       (              ...                              )
       (                     1                         )
       ( -s(k)                    c(k)                 )
       (                                 1             )
       (                                      ...      )
       (                                             1 )

where R(k) appears in rows and columns 1 and k+1.

Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is performed for the plane (k,z), giving P(k) the form

P(k) = ( 1                                             )
       (      ...                                      )
       (             1                                 )
       (                  c(k)                    s(k) )
       (                         1                     )
       (                              ...              )
       (                                     1         )
       (                 -s(k)                    c(k) )

where R(k) appears in rows and columns k and z. and the whole time, the rotations are performed without ever forming P(k) explicitly. no cap, that's the whole flex. 🔥

Parameters

SIDE — SIDE is CHARACTER*1. tells us whether the plane rotation matrix P gets applied to A on the left or on the right.
- = 'L': Left, compute A := P*A
- = 'R': Right, compute A:= A*P**T
PIVOT — PIVOT is CHARACTER*1. tells us which plane P(k) is a plane rotation matrix for.
- = 'V': Variable pivot, the plane (k,k+1)
- = 'T': Top pivot, the plane (1,k+1)
- = 'B': Bottom pivot, the plane (k,z)
DIRECT — DIRECT is CHARACTER*1. is P a forward or a backward sequence of plane rotations?
- = 'F': Forward, P = P(z-1)*...*P(2)*P(1)
- = 'B': Backward, P = P(1)*P(2)*...*P(z-1)
M — M is INTEGER. the number of rows of the matrix A. If m <= 1, an immediate return is effected (like, we're out, bye 🫶).
N — N is INTEGER. the number of columns of the matrix A. If n <= 1, an immediate return is effected (same energy, immediate exit).
C — C is DOUBLE PRECISION array, dimension (M-1) if SIDE = 'L', (N-1) if SIDE = 'R'. the cosines c(k) of the plane rotations. ✨
S — S is DOUBLE PRECISION array, dimension (M-1) if SIDE = 'L', (N-1) if SIDE = 'R'. the sines s(k) of the plane rotations. the 2-by-2 plane rotation part of the matrix P(k), R(k), has the form

R(k) = (  c(k)  s(k) )
       ( -s(k)  c(k) ).

A — A is COMPLEX*16 array, dimension (LDA,N). the M-by-N matrix A. On exit, A is overwritten by P*A if SIDE = 'R' or by A*P**T if SIDE = 'L'.
LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= max(1,M).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 subroutine `zlaswp` (`integer N, complex16, dimension( lda, ) A, integer LDA, integer K1, integer K2, integer, dimension( * ) IPIV, integer INCX`)

ZLASWP performs a series of row interchanges on a general rectangular matrix. swapping era, honestly. 🔁

Purpose:

ZLASWP performs a series of row interchanges on the matrix A.
One row interchange is initiated for each of rows K1 through K2 of A.

Parameters

N — N is INTEGER. the number of columns of the matrix A.
A — A is COMPLEX*16 array, dimension (LDA,N). On entry, the matrix of column dimension N to which the row interchanges will be applied. On exit, the permuted matrix. glow up achieved. 💖
LDA — LDA is INTEGER. the leading dimension of the array A.
K1 — K1 is INTEGER. the first element of IPIV for which a row interchange will be done.
K2 — K2 is INTEGER. (K2-K1+1) is the number of elements of IPIV for which a row interchange will be done.
IPIV — IPIV is INTEGER array, dimension (K1+(K2-K1)*abs(INCX)). the vector of pivot indices. Only the elements in positions K1 through K1+(K2-K1)*abs(INCX) of IPIV are accessed. IPIV(K1+(K-K1)*abs(INCX)) = L implies rows K and L are to be interchanged.
INCX — INCX is INTEGER. the increment between successive values of IPIV. If INCX is negative, the pivots are applied in reverse order. 🌀

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

Modified by
 R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA

🎀 subroutine `zlat2c` (`character UPLO, integer N, complex16, dimension( lda, ) A, integer LDA, complex, dimension( ldsa, * ) SA, integer LDSA, integer INFO`)

ZLAT2C converts a double complex triangular matrix to a complex triangular matrix. downgrade the precision, keep the drama. ✨

Purpose:

ZLAT2C converts a COMPLEX*16 triangular matrix, SA, to a COMPLEX
triangular matrix, A.

RMAX is the overflow for the SINGLE PRECISION arithmetic
ZLAT2C checks that all the entries of A are between -RMAX and
RMAX. If not the conversion is aborted and a flag is raised.

This is an auxiliary routine so there is no argument checking.

lowkey important: if any entry is outside the -RMAX … RMAX range, the conversion is aborted and a flag is raised. no argument checking here, bestie, it trusts you. 🫶

Parameters

UPLO — UPLO is CHARACTER*1.
- = 'U': A is upper triangular;
- = 'L': A is lower triangular.
N — N is INTEGER. the number of rows and columns of the matrix A. N >= 0.
A — A is COMPLEX*16 array, dimension (LDA,N). On entry, the N-by-N triangular coefficient matrix A.
LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= max(1,N).
SA — SA is COMPLEX array, dimension (LDSA,N). Only the UPLO part of SA is referenced. On exit, if INFO=0, the N-by-N coefficient matrix SA; if INFO>0, the content of the UPLO part of SA is unspecified. (mystery box vibes, we just don't know 🙈)
LDSA — LDSA is INTEGER. the leading dimension of the array SA. LDSA >= max(1,M).
INFO — INFO is INTEGER.
- = 0: successful exit.
- = 1: an entry of the matrix A is greater than the SINGLE PRECISION overflow threshold, in this case, the content of the UPLO part of SA in exit is unspecified.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 subroutine `zlatbs` (`character UPLO, character TRANS, character DIAG, character NORMIN, integer N, integer KD, complex16, dimension( ldab, ) AB, integer LDAB, complex16, dimension( ) X, double precision SCALE, double precision, dimension( * ) CNORM, integer INFO`)

ZLATBS solves a triangular banded system of equations. scaling included so nothing overflows, periodt. 💅

Purpose:

ZLATBS solves one of the triangular systems

   A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,

with scaling to prevent overflow, where A is an upper or lower
triangular band matrix.  Here A**T 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 will not cause
overflow, the Level 2 BLAS routine ZTBSV is called.  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.

so iconic: if A is singular (A(j,j) = 0 for some j), then s gets set to 0 and a non-trivial solution to A*x = 0 comes back. that's the plot twist. 😭

Parameters

UPLO — UPLO is CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS — TRANS is CHARACTER*1. Specifies the operation applied to A.
- = 'N': Solve A * x = s*b (No transpose)
- = 'T': Solve A**T * x = s*b (Transpose)
- = 'C': Solve A**H * x = s*b (Conjugate transpose)
DIAG — DIAG is CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN

NORMIN is CHARACTER*1 and it's lowkey just a vibe check for CNORM 💅 — it tells us whether the column norms are already loaded or whether we gotta go get them ourselves.

= 'Y': CNORM contains the column norms on entry
= 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.

N

N is INTEGER and it's the order of the matrix A, bestie. N >= 0. No negative orders, that's not a thing, periodt.

KD

KD is INTEGER — the number of subdiagonals or superdiagonals in the triangular matrix A. KD >= 0. So basically how wide the band is, iykyk.

AB

AB is COMPLEX*16 array, dimension (LDAB,N). It's the upper or lower triangular band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:

if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

Copy those indices exactly, bestie, they are NOT vibes-based 🌸

LDAB

LDAB is INTEGER — the leading dimension of the array AB. LDAB >= KD+1. Respect the leading dimension or perish 😭

X

X is COMPLEX*16 array, dimension (N).

On entry: the right hand side b of the triangular system.
On exit: X is overwritten by the solution vector x. She gets a whole glow-up in place 💖

SCALE

SCALE is DOUBLE PRECISION — the scaling factor s for the triangular system

A * x = s*b,  A**T * x = s*b,  or  A**H * 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. So like, zero scale is the "we could not do it, here's the null-space instead" energy 🫠

CNORM

CNORM is DOUBLE PRECISION array, dimension (N).

If NORMIN = 'Y', 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 = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm.

If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.

So she switches from input girlie to output girlie depending on NORMIN. Main character energy, honestly 👑

INFO

INFO is INTEGER

= 0: successful exit ✨
< 0: if INFO = -k, the k-th argument had an illegal value 💀

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

A rough bound on x is computed; if that is less than overflow, ZTBSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation. Safety girlie behavior fr fr 🧚

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 ZTBSV 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. Bestie that's the "we gave up and found the null vector instead" arc 😭

Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *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 ZTBSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).

✨ subroutine `zlatdf` ✨

subroutine zlatdf (integer IJOB, integer N, complex*16, dimension( ldz, * )
Z, integer LDZ, complex*16, dimension( * ) RHS, double precision RDSUM,
double precision RDSCAL, integer, dimension( * ) IPIV, integer,
dimension( * ) JPIV)

ZLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. She's a specialist, okay, she does one thing and she does it beautifully 🎀

Purpose:

ZLATDF computes the contribution to the reciprocal Dif-estimate by solving for x in Z * x = b, where b is chosen such that the norm of x is as large as possible. It is assumed that LU decomposition of Z has been computed by ZGETC2. On entry RHS = f holds the contribution from earlier solved sub-systems, and on return RHS = x.

The factorization of Z returned by ZGETC2 has the form Z = P * L * U * Q, where P and Q are permutation matrices. L is lower triangular with unit diagonal elements and U is upper triangular.

Parameters

IJOB

IJOB is INTEGER

IJOB = 2: First compute an approximative null-vector e of Z using ZGECON, e is normalized and solve for Zx = +-e - f with the sign giving the greater value of 2-norm(x). About 5 times as expensive as Default. Yeah bestie, that's the deluxe version, it costs 🫶
IJOB .ne. 2: Local look ahead strategy where all entries of the r.h.s. b is chosen as either +1 or -1. Default. This is the "we move" option.

N

N is INTEGER — the number of columns of the matrix Z.

Z

Z is COMPLEX*16 array, dimension (LDZ, N). On entry, the LU part of the factorization of the n-by-n matrix Z computed by ZGETC2: Z = P * L * U * Q

LDZ

LDZ is INTEGER — the leading dimension of the array Z. LDA >= max(1, N).

Yes, she says LDA in a routine where the arg is literally called LDZ. We are keeping it exactly as written, no notes 💀

RHS

RHS is COMPLEX*16 array, dimension (N).

On entry, RHS contains contributions from other subsystems.
On exit, RHS contains the solution of the subsystem with entries according to the value of IJOB (see above).

RDSUM

RDSUM is DOUBLE PRECISION

On entry, the sum of squares of computed contributions to the Dif-estimate under computation by ZTGSYL, where the scaling factor RDSCAL (see below) has been factored out.
On exit, the corresponding sum of squares updated with the contributions from the current sub-system.
If TRANS = 'T' RDSUM is not touched.
NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL.

RDSCAL

RDSCAL is DOUBLE PRECISION

On entry, scaling factor used to prevent overflow in RDSUM.
On exit, RDSCAL is updated w.r.t. the current contributions in RDSUM.
If TRANS = 'T', RDSCAL is not touched.
NOTE: RDSCAL only makes sense when ZTGSY2 is called by ZTGSYL.

IPIV

IPIV is INTEGER array, dimension (N). The pivot indices; for 1 <= i <= N, row i of the matrix has been interchanged with row IPIV(i).

JPIV

JPIV is INTEGER array, dimension (N). The pivot indices; for 1 <= j <= N, column j of the matrix has been interchanged with column JPIV(j).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

This routine is a further developed implementation of algorithm BSOLVE in [1] using complete pivoting in the LU factorization.

Contributors:

Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

References:

[1] Bo Kagstrom and Lars Westin, Generalized Schur Methods with Condition Estimators for Solving the Generalized Sylvester Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.

[2] Peter Poromaa, On Efficient and Robust Estimators for the Separation between two Regular Matrix Pairs with Applications in Condition Estimation. Report UMINF-95.05, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.

🌈 subroutine `zlatps` 🌈

subroutine zlatps (character UPLO, character TRANS, character DIAG,
character NORMIN, integer N, complex*16, dimension( * ) AP, complex*16,
dimension( * ) X, double precision SCALE, double precision, dimension( *
) CNORM, integer INFO)

ZLATPS solves a triangular system of equations with the matrix held in packed storage. Packed, bestie — no wasted space, she's efficient 💅

Purpose:

ZLATPS solves one of the triangular systems

A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,

with scaling to prevent overflow, where A is an upper or lower triangular matrix stored in packed form. Here A**T denotes the transpose of A, A**H denotes the conjugate 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 will not cause overflow, the Level 2 BLAS routine ZTPSV is called. 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. She literally refuses to overflow, iconic 🌸

Parameters

UPLO

UPLO — UPLO is CHARACTER*1. Okay so this literally specifies whether the matrix A is upper or lower triangular. 💅
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS — TRANS is CHARACTER*1. Specifies the operation applied to A. ✨
- = 'N': Solve A * x = s*b (No transpose)
- = 'T': Solve A**T * x = s*b (Transpose)
- = 'C': Solve A**H * x = s*b (Conjugate transpose)
DIAG — DIAG is CHARACTER*1. Specifies whether or not the matrix A is unit triangular. 👑
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN — NORMIN is CHARACTER*1. Specifies whether CNORM has been set or not.
- = 'Y': CNORM contains the column norms on entry
- = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.
N — N is INTEGER. The order of the matrix A. N >= 0. 🌸
AP — AP is COMPLEX*16 array, dimension (N*(N+1)/2). The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows:
- if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
- if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
X — X is COMPLEX*16 array, dimension (N). On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x. 🦋
SCALE — SCALE is DOUBLE PRECISION. The scaling factor s for the triangular system

A * x = s*b,  A**T * x = s*b,  or  A**H * 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.

CNORM — CNORM is DOUBLE PRECISION array, dimension (N).

If NORMIN = 'Y', 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 = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm.

If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.

INFO — INFO is INTEGER
- = 0: successful exit
- < 0: if INFO = -k, the k-th argument had an illegal value

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

Okay bestie, here's the deep dive. 💖

A rough bound on x is computed; if that is less than overflow, ZTPSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation.

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 ZTPSV 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.

Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *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 ZTPSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).

`subroutine zlatrd (character UPLO, integer N, integer NB, complex16, dimension( lda, ) A, integer LDA, double precision, dimension( * ) E, complex16, dimension( ) TAU, complex16, dimension( ldw, ) W, integer LDW)`

ZLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an unitary similarity transformation.

Purpose:

ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to Hermitian tridiagonal form by a unitary similarity transformation Q**H * A * Q, and returns the matrices V and W which are needed to apply the transformation to the unreduced part of A.

If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a matrix, of which the upper triangle is supplied; if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a matrix, of which the lower triangle is supplied.

This is an auxiliary routine called by ZHETRD.

Parameters

UPLO — UPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the Hermitian matrix A is stored:
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. The order of the matrix A.
NB — NB is INTEGER. The number of rows and columns to be reduced.
A — A is COMPLEX*16 array, dimension (LDA,N). On entry, the Hermitian matrix A. If UPLO = 'U', 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 = 'L', 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 UPLO = 'U', the last NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements above the diagonal with the array TAU, represent the unitary matrix Q as a product of elementary reflectors; if UPLO = 'L', the first NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements below the diagonal with the array TAU, represent the unitary matrix Q as a product of elementary reflectors. See Further Details.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (N-1). If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal elements of the last NB columns of the reduced matrix; if UPLO = 'L', E(1:nb) contains the subdiagonal elements of the first NB columns of the reduced matrix.
TAU — TAU is COMPLEX*16 array, dimension (N-1). The scalar factors of the elementary reflectors, stored in TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. See Further Details.
W — W is COMPLEX*16 array, dimension (LDW,NB). The n-by-nb matrix W required to update the unreduced part of A.
LDW — LDW is INTEGER. The leading dimension of the array W. LDW >= max(1,N).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

If UPLO = 'U', the matrix Q is represented as a product of elementary reflectors

Q = H(n) H(n-1) . . . H(n-nb+1).

Each H(i) has the form

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

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

If UPLO = 'L', the matrix Q is represented as a product of elementary reflectors

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

Each H(i) has the form

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

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

The elements of the vectors v together form the n-by-nb matrix V which is needed, with W, to apply the transformation to the unreduced part of the matrix, using a Hermitian rank-2k update of the form: A := A - V*W**H - W*V**H.

The contents of A on exit are illustrated by the following examples with n = 5 and nb = 2:

if UPLO = 'U':                       if UPLO = 'L':

               (  a   a   a   v4  v5 )              (  d                  )
               (      a   a   v4  v5 )              (  1   d              )
               (          a   1   v5 )              (  v1  1   a          )
               (              d   1  )              (  v1  v2  a   a      )
               (                  d  )              (  v1  v2  a   a   a  )

where d denotes a diagonal element of the reduced matrix, a denotes an element of the original matrix that is unchanged, and vi denotes an element of the vector defining H(i). 💅

✨ `subroutine zlatrs` ✨

subroutine zlatrs (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, complex*16, dimension( lda, * ) A, integer
    LDA, complex*16, dimension( * ) X, double precision SCALE, double
    precision, dimension( * ) CNORM, integer INFO)

ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow. No cap. 💖

Purpose 🎯

ZLATRS solves one of the triangular systems

   A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,

with scaling to prevent overflow. Here A is an upper or lower triangular matrix, A**T denotes the transpose of A, A**H denotes the conjugate 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 will not cause overflow, the Level 2 BLAS routine ZTRSV is called. 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. 💀

Parameters 🌸

UPLO — UPLO is CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS — TRANS is CHARACTER*1. Specifies the operation applied to A.
- = 'N': Solve A * x = s*b (No transpose)
- = 'T': Solve A**T * x = s*b (Transpose)
- = 'C': Solve A**H * x = s*b (Conjugate transpose)
DIAG — DIAG is CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN — NORMIN is CHARACTER*1. Specifies whether CNORM has been set or not.
- = 'Y': CNORM contains the column norms on entry
- = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.
N — N is INTEGER. The order of the matrix A. N >= 0.
A — A is COMPLEX*16 array, dimension (LDA,N). The triangular matrix A. If UPLO = 'U', 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 = 'L', 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 = 'U', the diagonal elements of A are also not referenced and are assumed to be 1.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max (1,N).
X — X is COMPLEX*16 array, dimension (N). On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.
SCALE — SCALE is DOUBLE PRECISION. The scaling factor s for the triangular system A * x = s*b, A**T * x = s*b, or A**H * 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.
CNORM — CNORM is DOUBLE PRECISION array, dimension (N).

If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.

INFO — INFO is INTEGER.
- = 0: successful exit
- < 0: if INFO = -k, the k-th argument had an illegal value

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Further Details 🔍

Okay bestie, here’s the tea on how ZLATRS does its thing. 💅

A rough bound on x is computed; if that is less than overflow, ZTRSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation.

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 ZTRSV 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. 😭

Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *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 ZTRSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow). Slay. 💅

✨ `subroutine zlauu2` ✨

subroutine zlauu2 (character UPLO, integer N, complex*16, dimension( lda, *
    ) A, integer LDA, integer INFO)

ZLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm).

Purpose 🎯

ZLAUU2 computes the product U * U**H or L**H * L, where the triangular factor U or L is stored in the upper or lower triangular part of the array A.

If UPLO = 'U' or 'u' then the upper triangle of the result is stored, overwriting the factor U in A. If UPLO = 'L' or 'l' then the lower triangle of the result is stored, overwriting the factor L in A.

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

Parameters 🌸

UPLO — UPLO is CHARACTER*1. Specifies whether the triangular factor stored in the array A is upper or lower triangular:
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. The order of the triangular factor U or L. N >= 0.
A — A is COMPLEX*16 array, dimension (LDA,N). On entry, the triangular factor U or L. On exit, if UPLO = 'U', the upper triangle of A is overwritten with the upper triangle of the product U * U**H; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**H * L.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,N).
INFO — INFO is INTEGER.
- = 0: successful exit
- < 0: if INFO = -k, the k-th argument had an illegal value

Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

✨ `subroutine zlauum` ✨

subroutine zlauum (character UPLO, integer N, complex*16, dimension( lda, *
    ) A, integer LDA, integer INFO)

ZLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm).

Purpose 🎯

ZLAUUM computes the product U * U**H or L**H * L, where the triangular factor U or L is stored in the upper or lower triangular part of the array A.

This is the blocked form of the algorithm, calling Level 3 BLAS.

Parameters 🌸

UPLO — UPLO is CHARACTER*1. Specifies whether the triangular factor stored in the array A is upper or lower triangular:
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. The order of the triangular factor U or L. N >= 0.
A —
A — COMPLEX*16 array, dimension (LDA,N). On entry bestie this is the triangular factor U or L. On exit, if UPLO = 'U', the upper triangle of A gets overwritten with the upper triangle of the product U * U**H; if UPLO = 'L', the lower triangle of A gets overwritten with the lower triangle of the product L**H * L. She's giving glow-up-and-then-rewrite 💅
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N). No shortcuts, she needs room.
INFO — INTEGER.
- = 0: successful exit ✨
- < 0: if INFO = -k, the k-th argument had an illegal value 😭

🌸 Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `zrot` (`integer N, complex16, dimension( ) CX, integer INCX, complex16, dimension( ) CY, integer INCY, double precision C, complex*16 S`)

ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors. Iconic little spin move 🌀

Purpose:

     ZROT   applies a plane rotation, where the cos (C) is real and the
     sin (S) is complex, and the vectors CX and CY are complex.

Parameters

N — N is INTEGER. The number of elements in the vectors CX and CY. Counting bestie, we love to see it 🔢
CX — CX is COMPLEX*16 array, dimension (N). On input, the vector X. On output, CX is overwritten with C*X + S*Y. She's been transformed, periodt ✨
INCX — INCX is INTEGER. The increment between successive values of CX. INCX <> 0. Zero is not the vibe, don't do it 🚫
CY — CY is COMPLEX*16 array, dimension (N). On input, the vector Y. On output, CY is overwritten with -CONJG(S)*X + C*Y. Conjugate jumped out, we love a plot twist 😭
INCY — INCY is INTEGER. The increment between successive values of CY. INCX <> 0. (yes bestie, the page literally says INCX right there in the INCY entry — we are keeping it exactly as written, no edits, no notes 💀)
C — C is DOUBLE PRECISION. Real one 🫶
S — S is COMPLEX*16. C and S define a rotation

   [  C          S  ]
   [ -conjg(S)   C  ]

where C*C + S*CONJG(S) = 1.0. That's the whole entire aesthetic, normalized and flawless 💅

🌈 Author 🌈

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `zspmv` (`character UPLO, integer N, complex16 ALPHA, complex16, dimension( * ) AP, complex16, dimension( ) X, integer INCX, complex16 BETA, complex16, dimension( * ) Y, integer INCY`)

ZSPMV computes a matrix-vector product for complex vectors using a complex symmetric packed matrix. Packed! Like a little suitcase of numbers 🎀

Purpose:

     ZSPMV  performs the matrix-vector operation

        y := alpha*A*x + beta*y,

     where alpha and beta are scalars, x and y are n element vectors and
     A is an n by n symmetric matrix, supplied in packed form.

Parameters

UPLO — UPLO is CHARACTER*1. On entry, UPLO specifies whether the upper or lower triangular part of the matrix A is supplied in the packed array AP as follows:
- UPLO = 'U' or 'u' — The upper triangular part of A is supplied in AP.
- UPLO = 'L' or 'l' — The lower triangular part of A is supplied in AP.

Unchanged on exit. She comes in, she says her piece, she does not change 💖

N — N is INTEGER. On entry, N specifies the order of the matrix A. N must be at least zero. Unchanged on exit.
ALPHA — ALPHA is COMPLEX*16. On entry, ALPHA specifies the scalar alpha. Unchanged on exit.
AP — AP is COMPLEX*16 array, dimension at least ( ( N*( N + 1 ) )/2 ). Before entry, with UPLO = 'U' or 'u', the array AP must contain the upper triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) respectively, and so on. Before entry, with UPLO = 'L' or 'l', the array AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) respectively, and so on. Unchanged on exit. Column by column bestie, she is stacked 🧚
X — X is COMPLEX*16 array, dimension at least ( 1 + ( N - 1 )*abs( INCX ) ). Before entry, the incremented array X must contain the N-element vector x. Unchanged on exit.
INCX — INCX is INTEGER. On entry, INCX specifies the increment for the elements of X. INCX must not be zero. Unchanged on exit.
BETA — BETA is COMPLEX*16. On entry, BETA specifies the scalar beta. When BETA is supplied as zero then Y need not be set on input. Unchanged on exit. Free pass when she's zero, we love an easy out 🙈
Y — Y is COMPLEX*16 array, dimension at least ( 1 + ( N - 1 )*abs( INCY ) ). Before entry, the incremented array Y must contain the n element vector y. On exit, Y is overwritten by the updated vector y.
INCY — INCY is INTEGER. On entry, INCY specifies the increment for the elements of Y. INCY must not be zero. Unchanged on exit.

💜 Author 💜

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `zspr` (`character UPLO, integer N, complex16 ALPHA, complex16, dimension( * ) X, integer INCX, complex16, dimension( ) AP`)

ZSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix. Rank-1 update, main character energy, she adds exactly one outer product and leaves 🔥

Purpose:

     ZSPR    performs the symmetric rank 1 operation

        A := alpha*x*x**H + A,

     where alpha is a complex scalar, x is an n element vector and A is an
     n by n symmetric matrix, supplied in packed form.

Parameters

UPLO — UPLO is CHARACTER*1. On entry, UPLO specifies whether the upper or lower triangular part of the matrix A is supplied in the packed array AP as follows:
- UPLO = 'U' or 'u' — The upper triangular part of A is supplied in AP.
- UPLO = 'L' or 'l' — The lower triangular part of A is supplied in AP.

Unchanged on exit.

N — N is INTEGER. On entry, N specifies the order of the matrix A. N must be at least zero. Unchanged on exit.
ALPHA — ALPHA is COMPLEX*16. On entry, ALPHA specifies the scalar alpha. Unchanged on exit.
X — X is COMPLEX*16 array, dimension at least ( 1 + ( N - 1 )*abs( INCX ) ). Before entry, the incremented array X must contain the N-element vector x. Unchanged on exit.
INCX — INCX is INTEGER. On entry, INCX specifies the increment for the elements of X. INCX must not be zero. Unchanged on exit.
AP — AP is COMPLEX*16 array, dimension at least ( ( N*( N + 1 ) )/2 ). Before entry, with UPLO = 'U' or 'u', the array AP must contain the upper triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) respectively, and so on. On exit, the array AP is overwritten by the upper triangular part of the updated matrix. Before entry, with UPLO = 'L' or 'l', the array AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) respectively, and so on. On exit, the array AP is overwritten by the lower triangular part of the updated matrix. Note that the imaginary parts of the diagonal elements need not be set, they are assumed to be zero, and on exit they are set to zero. Imaginary diagonal? Never heard of her, assumed zero, and she stays zero ✨

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `ztprfb` (`character SIDE, character TRANS, character DIRECT, character STOREV, integer M, integer N, integer K, integer L, complex16, dimension( ldv, ) V, integer LDV, complex16, dimension( ldt, ) T, integer LDT, complex16, dimension( lda, ) A, integer LDA, complex16, dimension( ldb, ) B, integer LDB, complex16, dimension( ldwork, ) WORK, integer LDWORK`)

ZTPRFB applies a complex 'triangular-pentagonal' block reflector to a complex matrix, which is composed of two blocks. Two blocks?? It's giving duo, it's giving we come as a pair 🦋

Purpose:

     ZTPRFB applies a complex 'triangular-pentagonal' block reflector H or its
     conjugate transpose H**H to a complex matrix C, which is composed of two
     blocks A and B, either from the left or right.

Parameters

SIDE — SIDE is CHARACTER*1.
- = 'L': apply H or H**H from the Left
- = 'R': apply H or H**H from the Right
TRANS — TRANS is CHARACTER*1.
- = 'N': apply H (No transpose)
- = 'C': apply H**H (Conjugate transpose)
DIRECT — DIRECT is CHARACTER*1. Indicates how H is formed from a product of elementary reflectors.
- = 'F': H = H(1) H(2) . . . H(k) (Forward)
- = 'B': H = H(k) . . . H(2) H(1) (Backward)

Forward or backward, both valid, both slay 💅

STOREV — STOREV is CHARACTER*1. Indicates how the vectors which define the elementary reflectors are stored:
- = 'C': Columns
- = 'R': Rows
M — M is INTEGER. The number of rows of the matrix B. M >= 0.
N — N is INTEGER. The number of columns of the matrix B. N >= 0.
K

okay bestie, we're picking up right in the middle of the parameter list — no re-intro, just straight into the tea 💅✨

K (INTEGER) — the order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. K >= 0.
L (INTEGER) — the order of the trapezoidal part of V. K >= L >= 0. See Further Details.
V (COMPLEX*16 array) — dimension:
- (LDV,K) if STOREV = 'C'
- (LDV,M) if STOREV = 'R' and SIDE = 'L'
- (LDV,N) if STOREV = 'R' and SIDE = 'R' The pentagonal matrix V, which contains the elementary reflectors H(1), H(2), ..., H(K). See Further Details.
LDV (INTEGER) — the leading dimension of the array V.
- If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
- if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
- if STOREV = 'R', LDV >= K.
T (COMPLEX*16 array) — dimension (LDT,K). The triangular K-by-K matrix T in the representation of the block reflector.
LDT (INTEGER) — the leading dimension of the array T. LDT >= K.
A (COMPLEX*16 array) — dimension:
- (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' On entry, the K-by-N or M-by-K matrix A. On exit, A is overwritten by the corresponding block of H*C or H**H*C or C*H or C*H**H. See Further Details.
LDA (INTEGER) — the leading dimension of the array A.
- If SIDE = 'L', LDA >= max(1,K);
- If SIDE = 'R', LDA >= max(1,M).
B (COMPLEX*16 array) — dimension (LDB,N). On entry, the M-by-N matrix B. On exit, B is overwritten by the corresponding block of H*C or H**H*C or C*H or C*H**H. See Further Details.
LDB (INTEGER) — the leading dimension of the array B. LDB >= max(1,M).
WORK (COMPLEX*16 array) — dimension:
- (LDWORK,N) if SIDE = 'L',
- (LDWORK,K) if SIDE = 'R'.
LDWORK (INTEGER) — the leading dimension of the array WORK.
- If SIDE = 'L', LDWORK >= K;
- if SIDE = 'R', LDWORK >= M.

👑 Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🌈 Further Details: 🌈

The matrix C is a composite matrix formed from blocks A and B. The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, and if SIDE = 'L', A is of size K-by-N.

If SIDE = 'R' and DIRECT = 'F', C = [A B].
If SIDE = 'L' and DIRECT = 'F', then:

C = [A]
    [B]

If SIDE = 'R' and DIRECT = 'B', C = [B A].
If SIDE = 'L' and DIRECT = 'B', then:

C = [B]
    [A]

The pentagonal matrix V is composed of a rectangular block V1 and a trapezoidal block V2. The size of the trapezoidal block is determined by the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; if L=0, there is no trapezoidal block, thus V = V1 is rectangular.

If DIRECT = 'F' and STOREV = 'C':

V = [V1]
    [V2]

V2 is upper trapezoidal (first L rows of K-by-K upper triangular)

If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
- V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
If DIRECT = 'B' and STOREV = 'C':

V = [V2]
    [V1]

V2 is lower trapezoidal (last L rows of K-by-K lower triangular)

If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
- V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.

📜 Author 📜

Generated automatically by Doxygen for LAPACK from the source code.

and that's a wrap, bestie! this was version 3.11.0, dated Sun Nov 27 2022, and it's the complex16OTHERauxiliary(3) manual — slay 💅✨🦋

📦 where this page comes from

debian bookworm: liblapack-doc 3.11.0-2 · BSD-3-Clause
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-10-14: a language model's rewrite of the original, under its terms

complex16OTHERauxiliary(3)

✨ complex16OTHERauxiliary(3) ✨

🎀 NAME 🎀

✨ SYNOPSIS ✨

🧚♀️ Functions 🧚♀️

🌈 Detailed Description 🌈

💅 Function Documentation 💅

🫧 subroutine clag2z (integer M, integer N, complex, dimension( ldsa, * ) SA, integer LDSA, complex*16, dimension( lda, * ) A, integer LDA, integer INFO)

✨ Purpose

🎀 Parameters

👑 Author

🦋 double precision function dzsum1 (integer N, complex*16, dimension( * ) CX, integer INCX)

✨ Purpose

🎀 Parameters

👑 Author

💜 Contributors:

🧚 integer function ilazlc (integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA)

✨ Purpose

🎀 Parameters

👑 Author

🌸 integer function ilazlr (integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA)

✨ Purpose

🎀 Parameters

👑 Author

🍓 subroutine zdrscl (integer N, double precision SA, complex*16, dimension( * ) SX, integer INCX)

✨ Purpose

🎀 Parameters

👑 Author

✨ Purpose

🎀 Parameters

✨ Author ✨

💜 Further Details: 💜

🌸 subroutine zlacgv (integer N, complex*16, dimension( * ) X, integer INCX) 🌸

Purpose:

Parameters

Author

💜 subroutine zlacn2 (integer N, complex*16, dimension( * ) V, complex*16, dimension( * ) X, double precision EST, integer KASE, integer, dimension( 3 ) ISAVE) 💜

Purpose:

Parameters

Author

Further Details:

Contributors:

References:

🎀 subroutine zlacon (integer N, complex*16, dimension( n ) V, complex*16, dimension( n ) X, double precision EST, integer KASE) 🎀

Purpose:

Parameters

Author

Further Details:

Contributors:

References:

💖 subroutine zlacp2 (character UPLO, integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB) 💖

Purpose:

Parameters

Author

🌈 subroutine zlacpy (character UPLO, integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB) 🌈

Purpose:

💖 subroutine zlacrm 💖

🌈 subroutine zlacrt 🌈

🧚 complex*16 function zladiv 🧚

🔥 subroutine zlaein 🔥

🌸 subroutine zlaev2 🌸

🍓 subroutine zlag2c 🍓

🫶 Author

✨ subroutine zlags2

🌈 Purpose

🎀 Parameters

👑 Author

🌟 subroutine zlagtm

🌈 Purpose

🎀 Parameters

👑 Author

✨ subroutine zlahqr

🌈 Purpose

🎀 Parameters

👑 Author

🎀 Contributors:

✨ subroutine zlahr2 (integer N, integer K, integer NB, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( nb ) TAU, complex*16, dimension( ldt, nb ) T, integer LDT, complex*16, dimension( ldy, nb ) Y, integer LDY)

💜 subroutine zlaic1 (integer JOB, integer J, complex*16, dimension( j ) X, double precision SEST, complex*16, dimension( j ) W, complex*16 GAMMA, double precision SESTPR, complex*16 S, complex*16 C)

🦋 double precision function zlangt (character NORM, integer N, complex*16, dimension( * ) DL, complex*16, dimension( * ) D, complex*16, dimension( * ) DU)

🎀 double precision function zlanhb (character NORM, character UPLO, integer N, integer K, complex*16, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) WORK)

🫧 `subroutine clag2z (integer M, integer N, complex, dimension( ldsa, * ) SA, integer LDSA, complex16, dimension( lda, ) A, integer LDA, integer INFO)`

🦋 `double precision function dzsum1 (integer N, complex16, dimension( ) CX, integer INCX)`

🧚 `integer function ilazlc (integer M, integer N, complex16, dimension( lda, ) A, integer LDA)`

🌸 `integer function ilazlr (integer M, integer N, complex16, dimension( lda, ) A, integer LDA)`

🍓 `subroutine zdrscl (integer N, double precision SA, complex16, dimension( ) SX, integer INCX)`

🌸 `subroutine zlacgv (integer N, complex16, dimension( ) X, integer INCX)` 🌸

💜 `subroutine zlacn2 (integer N, complex16, dimension( ) V, complex16, dimension( ) X, double precision EST, integer KASE, integer, dimension( 3 ) ISAVE)` 💜

🎀 `subroutine zlacon (integer N, complex16, dimension( n ) V, complex16, dimension( n ) X, double precision EST, integer KASE)` 🎀

💖 `subroutine zlacp2 (character UPLO, integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, complex16, dimension( ldb, ) B, integer LDB)` 💖

🌈 `subroutine zlacpy (character UPLO, integer M, integer N, complex16, dimension( lda, ) A, integer LDA, complex16, dimension( ldb, ) B, integer LDB)` 🌈

✨ `subroutine zlags2`

🌟 `subroutine zlagtm`

✨ `subroutine zlahqr`

✨ subroutine zlahr2 (integer N, integer K, integer NB, complex16, dimension( lda, ) A, integer LDA, complex16, dimension( nb ) TAU, complex16, dimension( ldt, nb ) T, integer LDT, complex*16, dimension( ldy, nb ) Y, integer LDY)

💜 subroutine zlaic1 (integer JOB, integer J, complex16, dimension( j ) X, double precision SEST, complex16, dimension( j ) W, complex16 GAMMA, double precision SESTPR, complex16 S, complex*16 C)

🦋 double precision function zlangt (character NORM, integer N, complex16, dimension( ) DL, complex16, dimension( ) D, complex16, dimension( ) DU)

🎀 double precision function zlanhb (character NORM, character UPLO, integer N, integer K, complex16, dimension( ldab, ) AB, integer LDAB, double precision, dimension( * ) WORK)

🦋 `double precision function zlanhp` 🦋

⭐ `double precision function zlanhs` ⭐

🎀 `double precision function zlanht` 🎀

💜 `double precision function zlansb` 💜

🌸 `double precision function zlansp`

🌸 `double precision function zlantb`

🌸 `double precision function zlantp`

🌸 `double precision function zlantr`

`zlapll` — measuring the linear dependence of two vectors 💅

`zlapmr` — rearranging rows like a boss 👑

`zlapmt` — permuting columns, periodt 💅

`zlaqhb` — scaling a Hermitian band matrix, no cap 💅

`zlaqhp` — scaling a Hermitian matrix in packed form 💅

subroutine `zlaqp2` (`integer M`, `integer N`, `integer OFFSET`, `complex16, dimension( lda, ) A`, `integer LDA`, `integer, dimension( * ) JPVT`, `complex16, dimension( ) TAU`, `double precision, dimension( * ) VN1`, `double precision, dimension( * ) VN2`, `complex16, dimension( ) WORK`)

`subroutine zlaqr1 (integer N, complex16, dimension( ldh, ) H, integer LDH, complex16 S1, complex16 S2, complex16, dimension( ) V)`

💜 subroutine `zlaqr3` 💜