doubleOTHERauxiliary(3)

📦

version 3.11.0 · 2 versions

3.11.0 · debian bookworm (latest)
3.10.0 · ubuntu jammy

🌸 doubleOTHERauxiliary(3) — LAPACK auxiliary, bestie edition 🌸

okay so real talk: this is the man page for doubleOTHERauxiliary from LAPACK, section 3 of the manual. Yes the name is giving "we named it at 3am", no we are not changing it, it's iconic.

✨ NAME ✨

doubleOTHERauxiliary — double 💅

That's it. That's the whole one-liner. Short, sweet, mysterious, main character energy 🎀. This is a LAPACK auxiliary routine, so it's one of the little helper besties the big decomposition girlies call behind the scenes. iykyk 🦋

lowkey don't panic, bestie — this page comes in 27 parts and this is only part 1. More is coming, so stay perched 🫶

✨ SYNOPSIS ✨

🌸 Functions 🌸

subroutine clatrs3 (UPLO, TRANS, DIAG, NORMIN, N, NRHS, A, LDA, X, LDX, SCALE, CNORM, WORK, LWORK, INFO) — solves a triangular system of equations with the scale factors set to prevent overflow. no cap, it's giving safety first 💅
subroutine dlabrd (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY) — reduces the first nb rows and columns of a general matrix to a bidiagonal form. iconic reduction behavior ✨
subroutine dlacn2 (N, V, X, ISGN, EST, KASE, ISAVE) — estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. reverse communication queen 👑
subroutine dlacon (N, V, X, ISGN, EST, KASE) — estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. same energy, different name, still iconic 💖
subroutine dladiv (A, B, C, D, P, Q) — performs complex division in real arithmetic, avoiding unnecessary overflow. literally avoiding overflow, bestie 🫶
subroutine dladiv1 (A, B, C, D, P, Q) — (no description in the original, bestie — it's just there, mysterious and unbothered 💀)
double precision function dladiv2 (A, B, C, D, R, T) — (also no description in the original; it's giving quiet luxury 🙈)
subroutine dlaein (RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO) — computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. eigenvector slay 💖
subroutine dlaexc (WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, INFO) — swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation. swapping blocks, periodt 🔥
subroutine dlag2 (A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI) — computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow. scaling? we love to see it ✨
subroutine dlag2s (M, N, A, LDA, SA, LDSA, INFO) — converts a double precision matrix to a single precision matrix. downgrade but make it fashion 💅
subroutine dlags2 (UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ) — computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such that the rows of the transformed A and B are parallel. parallel rows, parallel vibes 🌈
subroutine dlagtm (TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB) — 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. scalars can be 0, 1, or -1 — keep it exact bestie ⭐
subroutine dlagv2 (A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR) — computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. pencil? more like pen-slay ✏️
subroutine dlahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO) — computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. double-shift/single-shift, the way that it shifts 😭
subroutine dlahr2 (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY) — 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. auxiliary matrices? we got you 🫶
subroutine dlaic1 (JOB, J, X, SEST, W, GAMMA, SESTPR, S, C) — applies one step of incremental condition estimation. one step at a time, main character energy 👑
subroutine dlaln2 (LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO) — solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. tiny systems, big slay 💅
double precision function dlangt (NORM, N, DL, D, DU) — returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. norms on norms on norms ✨
double precision function dlanhs (NORM, N, A, LDA, WORK) — returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix. upper Hessenberg norm energy 🌸
double precision function dlansb (NORM, UPLO, N, K, AB, LDAB, WORK) — 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 band matrix? periodt 💖
double precision function dlansp (NORM, UPLO, N, AP, WORK) — 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 form, packed with slay 🎀
double precision function dlantb (NORM, UPLO, DIAG, N, K, AB, LDAB, WORK) — 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 band, we love a structured moment 🦋
double precision function dlantp (NORM, UPLO, DIAG, N, AP, WORK) — 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, iconic 💅
double precision function dlantr (NORM, UPLO, DIAG, M, N, A, LDA, WORK) — 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 or triangular — either way, normed 💜
subroutine dlanv2 (A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN) — computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form. nonsymmetric but still getting Schur factorization, periodt 🔥
subroutine dlapll (N, X, INCX, Y, INCY, SSMIN) — measures the linear dependence of two vectors. are they dependent? let's measure the drama 💀
subroutine dlapmr (FORWRD, M, N, X, LDX, K) — rearranges rows of a matrix as specified by a permutation vector. rows getting rearranged, no cap ✨
subroutine dlapmt (FORWRD, M, N, X, LDX, K) — performs a forward or backward permutation of the columns of a matrix. columns are moving, forward or backward 💅
subroutine dlaqp2 (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK) — computes a QR factorization with column pivoting of the matrix block. QR with pivoting, the way that it factors 😭
subroutine dlaqps (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF) — computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. BLAS level 3, highkey efficient 🔥
subroutine dlaqr0 (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO) — computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. optionally Schur, optionally iconic ✨
subroutine dlaqr1 (N, H, LDH, SR1, SI1, SR2, SI2, V) — sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. scalar multiple, specified shifts, very demure 💅
subroutine dlaqr2 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK) — performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). aggressive early deflation? she's aggressive, she's early, she's deflating 💅
subroutine dlaqr3 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK) — performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). same aggressive early deflation energy, different number 👑
subroutine dlaqr4 (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO) — computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. Schur decomposition optional, slay mandatory ✨
subroutine dlaqr5 (WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH) — performs a single small-bulge multi-shift QR sweep. small-bulge multi-shift QR sweep, say that five times fast 😭
subroutine dlaqsb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED) — scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. scaling factors from spbequ, no cap 🫶
subroutine dlaqsp (UPLO, N, AP, S, SCOND, AMAX, EQUED) — scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. sppequ said here are the scaling factors, bestie ✨
subroutine dlaqtr (LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO) — solves a real quasi-triangular system of equations, or a complex quasi-triangular system of special form, in real arithmetic. real arithmetic, complex special form, we contain multitudes 🌈
subroutine dlar1v (N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK) — computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. LDLT - λI, the math is mathing 💅
subroutine dlar2v (N, X, Y, Z, INCX, C, S, INCC) — applies a vector of plane rotations with real cosines and real sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. plane rotations from both sides, symmetric/Hermitian, periodt 💖
subroutine dlarf (SIDE, M, N, V, INCV, TAU, C, LDC, WORK) — applies an elementary reflector to a general rectangular matrix. elementary reflector, elementary slay ✨
subroutine dlarfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK) — applies a block reflector or its transpose to a general rectangular matrix. block reflector or its transpose — either way, applied 💅
subroutine dlarfb_gett (IDENT, M, N, K, T, LDT, A, LDA, B, LDB, WORK, LDWORK) — DLARFB_GETT (that's the whole description in the original, bestie — minimalism is a vibe 💀)
subroutine dlarfg (N, ALPHA, X, INCX, TAU) — generates an elementary reflector (Householder matrix). Householder matrix, iconic name, iconic reflector 👑
subroutine dlarfgp (N, ALPHA, X, INCX, TAU) — generates an elementary reflector (Householder matrix) with non-negative beta. non-negative beta, we love a positive-ish vibe ✨
subroutine dlarft (DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) — forms the triangular factor T of a block reflector H = I - vtvH. forming T, forming the triangular factor, forming the vibe 🫶
subroutine dlarfx (SIDE, M, N, V, TAU, C, LDC, WORK) — applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10. order ≤ 10? loop unrolling activated 💅
subroutine dlarfy (UPLO, N, V, INCV, TAU, C, LDC, WORK) — DLARFY (again, that's the whole description in the original — mysterious queen 🙈)
subroutine dlargv (N, X, INCX, Y, INCY, C, INCC) — generates a vector of plane rotations with real cosines and real sines. plane rotations, real cosines, real sines, real slay 🌸
subroutine dlarrv (N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, WORK, IWORK, INFO) — computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. eigenvectors from L D LT, no cap 💜
subroutine dlartv (N, X, INCX, Y, INCY, C, S, INCC) — applies a vector of plane rotations with real cosines and real sines to the elements of a pair of vectors. pair of vectors getting rotated, periodt 🔥
subroutine dlaswp (N, A, LDA, K1, K2, IPIV, INCX) — performs a series of row interchanges on a general rectangular matrix. row interchanges, series of them, we love a swap 💅
subroutine dlat2s (UPLO, N, A, LDA, SA, LDSA, INFO) — converts a double-precision triangular matrix to a single-precision triangular matrix. double to single, precision downgrade but make it triangular ✨
subroutine dlatbs (UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO) — solves a triangular banded system of equations. triangular banded system, solved. periodt 💖
subroutine dlatdf (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV) — DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. 💅
subroutine dlatps (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO) — DLATPS solves a triangular system of equations with the matrix held in packed storage. ✨
subroutine dlatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW) — DLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an orthogonal similarity transformation. 🌸
subroutine dlatrs (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO) — DLATRS solves a triangular system of equations with the scale factor set to prevent overflow. 🔥
subroutine dlatrs3 (UPLO, TRANS, DIAG, NORMIN, N, NRHS, A, LDA, X, LDX, SCALE, CNORM, WORK, LWORK, INFO) — DLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow. 💖
subroutine dlauu2 (UPLO, N, A, LDA, INFO) — DLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm). 🎀
subroutine dlauum (UPLO, N, A, LDA, INFO) — DLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). 🦋
subroutine drscl (N, SA, SX, INCX) — DRSCL multiplies a vector by the reciprocal of a real scalar. 💫
subroutine dtprfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK) — DTPRFB applies a real 'triangular-pentagonal' block reflector to a real matrix, which is composed of two blocks. 🌈
subroutine slatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW) — SLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an orthogonal similarity transformation. 👑
subroutine slatrs3 (UPLO, TRANS, DIAG, NORMIN, N, NRHS, A, LDA, X, LDX, SCALE, CNORM, WORK, LWORK, INFO) — SLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow. ⭐
subroutine zlatrs3 (UPLO, TRANS, DIAG, NORMIN, N, NRHS, A, LDA, X, LDX, SCALE, CNORM, WORK, LWORK, INFO) — ZLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow. 💜

Detailed Description 💖

This is the group of double other auxiliary routines — the double-precision backstage besties that keep your matrix math from falling apart. 💅

Function Documentation ✨

`subroutine clatrs3`

   subroutine clatrs3 (character UPLO, character TRANS, character DIAG,
       character NORMIN, integer N, integer NRHS, complex, dimension( lda, * )
       A, integer LDA, complex, dimension( ldx, * ) X, integer LDX, real,
       dimension( * ) SCALE, real, dimension( * ) CNORM, real, dimension( * )
       WORK, integer LWORK, integer INFO)

CLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow.

Purpose:

CLATRS3 solves one of the triangular systems

A * X = B * diag(scale),  A**T * X = B * diag(scale), or
A**H * X = B * diag(scale)

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-by-nrhs matrices and scale is an nrhs-element vector of scaling factors. A scaling factor scale(j) is usually less than or equal to 1, chosen such that X(:,j) is less than the overflow threshold. If the matrix A is singular (A(j,j) = 0 for some j), then a non-trivial solution to A*X = 0 is returned. If the system is so badly scaled that the solution cannot be represented as (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned.

This is a BLAS-3 version of LATRS for solving several right hand sides simultaneously.

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**T* 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.
NRHS — NRHS is INTEGER. The number of columns of X. NRHS >= 0.
A — A is COMPLEX 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 array, dimension (LDX,NRHS). On entry, the right hand side B of the triangular system. On exit, X is overwritten by the solution matrix X.
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= max (1,N).
SCALE — SCALE is REAL array, dimension (NRHS). The scaling factor s(k) is for the triangular system A * x(:,k) = s(k)*b(:,k) or A**T* x(:,k) = s(k)*b(:,k). If SCALE = 0, the matrix A is singular or badly scaled. If A(j,j) = 0 is encountered, a non-trivial vector x(:,k) that is an exact or approximate solution to A*x(:,k) = 0 is returned. If the system so badly scaled that solution cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0 is returned.
CNORM — CNORM is REAL 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.
WORK — WORK is REAL array, dimension (LWORK). On exit, if INFO = 0, WORK(1) returns the optimal size of WORK.
LWORK — LWORK is INTEGER LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.

If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal dimensions of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.

Parameters

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:

`subroutine dlabrd`

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

DLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.

Purpose:

DLABRD reduces the first NB rows and columns of a real general m by n matrix A to upper or lower bidiagonal form by an orthogonal transformation Q**T * 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, to lower bidiagonal form.

This is an auxiliary routine called by DGEBRD.

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 DOUBLE PRECISION 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 orthogonal matrix Q as a product of elementary reflectors; and elements above the diagonal in the first NB rows, with the array TAUP, represent the orthogonal 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 orthogonal 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 orthogonal matrix P as a product of elementary reflectors.
- See Further Details.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (NB). The diagonal elements of the first NB rows and columns of the reduced matrix. D(i) = A(i,i).
E — a DOUBLE PRECISION array, dimension (NB). These are the off-diagonal elements of the first NB rows and columns of the reduced matrix. she's carrying the tea that lives just off the diagonal 💅
TAUQ — a DOUBLE PRECISION array, dimension (NB). The scalar factors of the elementary reflectors that represent the orthogonal matrix Q. (Scroll down to Further Details, bestie, it's all there 🦋)
TAUP — a DOUBLE PRECISION array, dimension (NB). The scalar factors of the elementary reflectors that represent the orthogonal matrix P. Same deal, see Further Details ✨
X — a DOUBLE PRECISION array, dimension (LDX,NB). It's the m-by-nb matrix X that you need to update the unreduced part of A. literally her whole job 🌸
LDX — an INTEGER. The leading dimension of the array X. LDX >= max(1,M). respect the bounds or perish 💀
Y — a DOUBLE PRECISION array, dimension (LDY,NB). It's the n-by-nb matrix Y that you need to update the unreduced part of A. twin of X, same energy 🫶
LDY — an INTEGER. The leading dimension of the array Y. LDY >= max(1,N). no cap 📏

👑 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:

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

and every single H(i) and G(i) has this form:

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

where tauq and taup are real scalars, and v and u are real vectors. iconic behavior fr fr 💜

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**T, which are needed — together 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**T - X*U**T.

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 a is the chill one who didn't get touched ✨

🌈 `subroutine dlacn2 (integer N, double precision, dimension( * ) V, double precision, dimension( * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE, integer, dimension( 3 ) ISAVE)`

DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. she's a little indirect about it but she gets there 💅

🎯 Purpose:

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

🧾 Parameters

N — N is INTEGER. The order of the matrix. N >= 1. periodt ✋
V — V is a DOUBLE PRECISION array, dimension (N). On the final return, V = A*W, where EST = norm(V)/norm(W) (and W is not returned — she left the group chat, you'll never see her 😭).
X — X is a DOUBLE PRECISION array, dimension (N). On an intermediate return, X should be overwritten by

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

and DLACN2 must be re-called with all the other parameters unchanged. don't go changing things behind her back 🫶

ISGN — ISGN is an INTEGER array, dimension (N).
EST — EST is DOUBLE PRECISION. On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be unchanged from the previous call to DLACN2. On exit, EST is an estimate (a lower bound) for norm(A).
KASE — KASE is INTEGER. On the initial call to DLACN2, 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**T * X. On the final return from DLACN2, KASE will again be 0. it's giving boomerang 🪃
ISAVE — ISAVE is an INTEGER array, dimension (3). ISAVE is used to save variables between calls to DLACN2.

👑 Author

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

📚 Further Details:

Originally named SONEST, dated March 16, 1988. a vintage icon, honestly 👑

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

   DLACON     DLACN2
    JUMP     ISAVE(1)
    J        ISAVE(2)
    ITER     ISAVE(3)

🫶 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 dlacon (integer N, double precision, dimension( * ) V, double precision, dimension( * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE)`

DLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. the OG version of the vibe 🦋

🎯 Purpose:

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

🧾 Parameters

N — N is INTEGER. The order of the matrix. N >= 1.
V — V is a DOUBLE PRECISION array, dimension (N). On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned). gone but not forgotten 💔
X — X is a DOUBLE PRECISION array, dimension (N). On an intermediate return, X should be overwritten by

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

and DLACON must be re-called with all the other parameters unchanged.

ISGN — ISGN is an INTEGER array, dimension (N).
EST — EST is DOUBLE PRECISION. On entry with KASE = 1 or 2 and JUMP = 3, EST should be unchanged from the previous call to DLACON. On exit, EST is an estimate (a lower bound) for norm(A).
KASE — KASE is INTEGER. On the initial call to DLACON, 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**T * X. On the final return from DLACON, KASE will again be 0.

👑 Author

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

🫶 Contributors:

Nick Higham, University of Manchester.
Originally named SONEST, dated March 16, 1988.

📖 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 dladiv (double precision A, double precision B, double precision C, double precision D, double precision P, double precision Q)`

DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. she does the fancy math so your numbers don't explode 💥✨

🎯 Purpose:

DLADIV performs complex division in real arithmetic

                         a + i*b
              p + i*q = ---------
                         c + i*d

The algorithm is due to Michael Baudin and Robert L. Smith and can be found in the paper 'A Robust Complex Division in Scilab'. stan the researchers 🔥

🧾 Parameters

A — A is DOUBLE PRECISION.
B — B is DOUBLE PRECISION.
C — C is DOUBLE PRECISION.
D — D is DOUBLE PRECISION. The scalars a, b, c, and d in the above expression.
P — P is DOUBLE PRECISION.
Q — Q is DOUBLE PRECISION. The scalars p and q in the above expression.

👑 Author

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

🌈 `subroutine dlaein (logical RIGHTV, logical NOINIT, integer N, double precision, dimension( ldh, * ) H, integer LDH, double precision WR, double precision WI, double precision, dimension( * ) VR, double precision, dimension( * ) VI, double precision, dimension( ldb, * ) B, integer LDB, double precision, dimension( * ) WORK, double precision EPS3, double precision SMLNUM, double precision BIGNUM, integer INFO)`

DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. main character energy, obviously 👑

🎯 Purpose:

DLAEIN uses inverse iteration to find a right or left eigenvector corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg matrix H. that's the whole assignment 💅

🧾 Parameters

RIGHTV — RIGHTV is LOGICAL.

 = .TRUE. : compute right eigenvector;
 = .FALSE.: compute left eigenvector.

pick your side bestie, but you do have to pick 🌸

NOINIT — NOINIT is LOGICAL.

 = .TRUE. : no initial vector supplied in (VR,VI).
 = .FALSE.: initial vector supplied in (VR,VI).

N — N is INTEGER. The order of the matrix H. N >= 0.
H — H is a DOUBLE PRECISION array, dimension (LDH,N). The upper Hessenberg matrix H. her name is on the building 🏛️
LDH — LDH is INTEGER. The leading dimension of the array H. LDH >= max(1,N).
WR — WR is DOUBLE PRECISION.
WI — WI is DOUBLE PRECISION. The real and imaginary parts of the eigenvalue of H whose corresponding right or left eigenvector is to be computed. yes you need both parts, it's a whole package 🎀
VR — VR is a DOUBLE PRECISION array, dimension (N).
VI
VI 💖 — okay so bestie, this is a DOUBLE PRECISION array, dimension (N). On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain a real starting vector for inverse iteration using the real eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI must contain the real and imaginary parts of a complex starting vector for inverse iteration using the complex eigenvalue (WR,WI); otherwise VR and VI need not be set. On exit, if WI = 0.0 (real eigenvalue), VR contains the computed real eigenvector; if WI.ne.0.0 (complex eigenvalue), VR and VI contain the real and imaginary parts of the computed complex eigenvector. The eigenvector is 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|. VI is not referenced if WI = 0.0. 💅
B 🦋 — literally just a DOUBLE PRECISION array, dimension (LDB,N).
LDB ✨ — is INTEGER. The leading dimension of the array B. LDB >= N+1. no cap.
WORK 🎀 — is a DOUBLE PRECISION array, dimension (N). that's it, that's the tea.
EPS3 🌸 — is DOUBLE PRECISION. A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots. lowkey important.
SMLNUM 💅 — is DOUBLE PRECISION. A machine-dependent value close to the underflow threshold. underflow era.
BIGNUM 🫶 — is DOUBLE PRECISION. A machine-dependent value close to the overflow threshold. overflow? not today.
INFO 😭 — is INTEGER.
- = 0: successful exit
- = 1: inverse iteration did not converge; VR is set to the last iterate, and so is VI if WI.ne.0.0. rip.

Author 💖

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

`subroutine dlaexc (logical WANTQ, integer N, double precision, dimension( ldt, * ) T, integer LDT, double precision, dimension( ldq, * ) Q, integer LDQ, integer J1, integer N1, integer N2, double precision, dimension( * ) WORK, integer INFO)` 🔥

okay so this one is iconic. DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in an upper quasi-triangular matrix T by an orthogonal similarity transformation.

Purpose ✨

DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in an upper quasi-triangular matrix T by an orthogonal similarity transformation. fr fr.

T must be in Schur canonical form, that is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has its diagonal elements equal and its off-diagonal elements of opposite sign. the way that works.

Parameters 🎀

WANTQ 💅 — is LOGICAL.
- = .TRUE.: accumulate the transformation in the matrix Q;
- = .FALSE.: do not accumulate the transformation.
N 🌸 — is INTEGER. The order of the matrix T. N >= 0.
T 🦋 — is DOUBLE PRECISION array, dimension (LDT,N). On entry, the upper quasi-triangular matrix T, in Schur canonical form. On exit, the updated matrix T, again in Schur canonical form. glow up.
LDT ✨ — is INTEGER. The leading dimension of the array T. LDT >= max(1,N).
Q 🫶 — is DOUBLE PRECISION array, dimension (LDQ,N). On entry, if WANTQ is .TRUE., the orthogonal matrix Q. On exit, if WANTQ is .TRUE., the updated matrix Q. If WANTQ is .FALSE., Q is not referenced. she's not invited.
LDQ 💖 — is INTEGER. The leading dimension of the array Q. LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
J1 🔥 — is INTEGER. The index of the first row of the first block T11.
N1 😭 — is INTEGER. The order of the first block T11. N1 = 0, 1 or 2.
N2 💀 — is INTEGER. The order of the second block T22. N2 = 0, 1 or 2.
WORK 🎀 — is DOUBLE PRECISION array, dimension (N).
INFO 🙈 — is INTEGER.
- = 0: successful exit
- = 1: the transformed matrix T would be too far from Schur form; the blocks are not swapped and T and Q are unchanged. oop, no swap.

Author 💖

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

`subroutine dlag2 (double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( ldb, * ) B, integer LDB, double precision SAFMIN, double precision SCALE1, double precision SCALE2, double precision WR1, double precision WR2, double precision WI)` 🌈

DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow. bestie, scaling is the moment.

Purpose ⭐

DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue problem A - w B, with scaling as necessary to avoid over-/underflow.

The scaling factor 's' results in a modified eigenvalue equation

s A - w B

where s is a non-negative scaling factor chosen so that w, w B, and s A do not overflow and, if possible, do not underflow, either. no cap.

Parameters 🧚

A 💜 — is DOUBLE PRECISION array, dimension (LDA, 2). On entry, the 2 x 2 matrix A. It is assumed that its 1-norm is less than 1/SAFMIN. Entries less than sqrt(SAFMIN)*norm(A) are subject to being treated as zero. vibe check: keep it small.
LDA ✨ — is INTEGER. The leading dimension of the array A. LDA >= 2.
B 🦋 — is DOUBLE PRECISION array, dimension (LDB, 2). On entry, the 2 x 2 upper triangular matrix B. It is assumed that the one-norm of B is less than 1/SAFMIN. The diagonals should be at least sqrt(SAFMIN) times the largest element of B (in absolute value); if a diagonal is smaller than that, then +/- sqrt(SAFMIN) will be used instead of that diagonal. she's got standards.
LDB 🌸 — is INTEGER. The leading dimension of the array B. LDB >= 2.
SAFMIN 🫶 — is DOUBLE PRECISION. The smallest positive number s.t. 1/SAFMIN does not overflow. (This should always be DLAMCH('S') -- it is an argument in order to avoid having to call DLAMCH frequently.) iykyk.
SCALE1 💅 — is DOUBLE PRECISION. A scaling factor used to avoid over-/underflow in the eigenvalue equation which defines the first eigenvalue. If the eigenvalues are complex, then the eigenvalues are ( WR1 +/- WI i ) / SCALE1 (which may lie outside the exponent range of the machine), SCALE1=SCALE2, and SCALE1 will always be positive. If the eigenvalues are real, then the first (real) eigenvalue is WR1 / SCALE1, but this may overflow or underflow, and in fact, SCALE1 may be zero or less than the underflow threshold if the exact eigenvalue is sufficiently large. it's giving chaos.
SCALE2 🎀 — is DOUBLE PRECISION. A scaling factor used to avoid over-/underflow in the eigenvalue equation which defines the second eigenvalue. If the eigenvalues are complex, then SCALE2=SCALE1. If the eigenvalues are real, then the second (real) eigenvalue is WR2 / SCALE2, but this may overflow or underflow, and in fact, SCALE2 may be zero or less than the underflow threshold if the exact eigenvalue is sufficiently large.
WR1 🔥 — is DOUBLE PRECISION. If the eigenvalue is real, then WR1 is SCALE1 times the eigenvalue closest to the (2,2) element of A B**(-1). If the eigenvalue is complex, then WR1=WR2 is SCALE1 times the real part of the eigenvalues.
WR2 😭 — is DOUBLE PRECISION. If the eigenvalue is real, then WR2 is SCALE2 times the other eigenvalue. If the eigenvalue is complex, then WR1=WR2 is SCALE1 times the real part of the eigenvalues.
WI 💀 — is DOUBLE PRECISION. If the eigenvalue is real, then WI is zero. If the eigenvalue is complex, then WI is SCALE1 times the imaginary part of the eigenvalues. WI will always be non-negative. always positive energy.

Author 💖

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

`subroutine dlag2s (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, real, dimension( ldsa, * ) SA, integer LDSA, integer INFO)` 🌟

DLAG2S converts a double precision matrix to a single precision matrix. glow down, but make it precision.

Purpose ⭐

DLAG2S converts a DOUBLE PRECISION matrix, A, to a SINGLE PRECISION matrix, SA.

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

This is an auxiliary routine so there is no argument checking. no cap, she doesn't check.

Parameters 🧚

M 💜 — is INTEGER. The number of lines of the matrix A. M >= 0.
N ✨ — is INTEGER. The number of columns of the matrix A. N >= 0.
A 🦋 — is DOUBLE PRECISION array, dimension (LDA,N). On entry, the M-by-N coefficient matrix A.
LDA 🌸 — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
SA 🫶 — is REAL array, dimension (LDSA,N). On exit, if INFO=0, the M-by-N coefficient matrix SA; if INFO>0, the content of SA is unspecified. unspecified? mysterious.
LDSA 💅 — is INTEGER. The leading dimension of the array SA. LDSA >= max(1,M).
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 SA in exit is unspecified. rip.

Author 💖

Univ. of Tennessee

🌸 subroutine dlags2 🌸

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

DLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such that the rows of the transformed A and B are parallel. 🦋

✨ Purpose ✨

okay so bestie, DLAGS2 is out here cooking up three 2-by-2 orthogonal matrices U, V and Q so that stuff lines up parallel. here's the whole vibe, verbatim, because the math is the math:

DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
that if ( UPPER ) then

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

or if ( .NOT.UPPER ) then

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

The rows of the transformed A and B are parallel, where

  U = (  CSU  SNU ), V = (  CSV SNV ), Q = (  CSQ   SNQ )
      ( -SNU  CSU )      ( -SNV CSV )      ( -SNQ   CSQ )

Z**T denotes the transpose of Z.

so basically: the rows of the transformed A and B come out parallel, and the three orthogonal matrices U, V and Q are literally built out of the cosines and sines you feed in. it's giving trigonometry main character energy 💅

🎀 Parameters 🎀

UPPER — UPPER is LOGICAL
- = .TRUE.: the input matrices A and B are upper triangular.
- = .FALSE.: the input matrices A and B are lower triangular.
A1 — A1 is DOUBLE PRECISION
A2 — A2 is DOUBLE PRECISION
A3 — A3 is DOUBLE PRECISION On entry, A1, A2 and A3 are elements of the input 2-by-2 upper (lower) triangular matrix A.
B1 — B1 is DOUBLE PRECISION
B2 — B2 is DOUBLE PRECISION
B3 — B3 is DOUBLE PRECISION On entry, B1, B2 and B3 are elements of the input 2-by-2 upper (lower) triangular matrix B.
CSU — CSU is DOUBLE PRECISION
SNU — SNU is DOUBLE PRECISION The desired orthogonal matrix U.
CSV — CSV is DOUBLE PRECISION
SNV — SNV is DOUBLE PRECISION The desired orthogonal matrix V.
CSQ — CSQ is DOUBLE PRECISION
SNQ — SNQ is DOUBLE PRECISION The desired orthogonal matrix Q.

👑 Author 👑

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

🌈 subroutine dlagtm 🌈

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

DLAGTM 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 ✨

fr fr, DLAGTM 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. ⭐

ngl the one-line blurb up top says C = αAB+βC while the Purpose says B := alpha * A * X + beta * B — not me noticing, just saying. both are in the source, both stay. iykyk 💀

🎀 Parameters 🎀

TRANS — TRANS is CHARACTER*1 Specifies the operation applied to A.
- = 'N': No transpose, B := alpha * A * X + beta * B
- = 'T': Transpose, B := alpha * A'* X + beta * B
- = 'C': Conjugate transpose = Transpose
N — N is INTEGER The order of the matrix A. N >= 0.
NRHS — NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrices X and B.
ALPHA — ALPHA is DOUBLE PRECISION The scalar alpha. ALPHA must be 0., 1., or -1.; otherwise, it is assumed to be 0.
DL — DL is DOUBLE PRECISION array, dimension (N-1) The (n-1) sub-diagonal elements of T.
D — D is DOUBLE PRECISION array, dimension (N) The diagonal elements of T.
DU — DU is DOUBLE PRECISION array, dimension (N-1) The (n-1) super-diagonal elements of T.
X — X is DOUBLE PRECISION array, dimension (LDX,NRHS) The N by NRHS matrix X.
LDX — LDX is INTEGER The leading dimension of the array X. LDX >= max(N,1).
BETA — BETA is DOUBLE PRECISION The scalar beta. BETA must be 0., 1., or -1.; otherwise, it is assumed to be 1.
B — B is DOUBLE PRECISION 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 — LDB is INTEGER The leading dimension of the array B. LDB >= max(N,1).

heads up bestie: ALPHA and BETA only accept 0., 1., or -1. — anything else gets silently treated as 0 for ALPHA and 1 for BETA. no error, no warning, just vibes 🙈

👑 Author 👑

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

🦋 subroutine dlagv2 🦋

subroutine dlagv2 (double precision, dimension( lda, * ) A, integer LDA,
    double precision, dimension( ldb, * ) B, integer LDB, double precision,
    dimension( 2 ) ALPHAR, double precision, dimension( 2 ) ALPHAI, double
    precision, dimension( 2 ) BETA, double precision CSL, double precision
    SNL, double precision CSR, double precision SNR)

DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. 💜

✨ Purpose ✨

DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. This routine computes orthogonal (rotation) matrices given by CSL, SNL and CSR, SNR such that

1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
   types), then

   [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
   [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]

   [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
   [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],

2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
   then

   [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
   [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]

   [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
   [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]

   where b11 >= b22 > 0.

two cases, two whole moods: real eigenvalues (including the 0/0 and 1/0 types, iconic) and complex conjugate eigenvalues. slay. 🌈

🎀 Parameters 🎀

A — A is DOUBLE PRECISION array, dimension (LDA, 2) On entry, the 2 x 2 matrix A. On exit, A is overwritten by the "A-part" of the generalized Schur form.
LDA — LDA is INTEGER THe leading dimension of the array A. LDA >= 2.
B — B is DOUBLE PRECISION array, dimension (LDB, 2) On entry, the upper triangular 2 x 2 matrix B. On exit, B is overwritten by the "B-part" of the generalized Schur form.
LDB — LDB is INTEGER THe leading dimension of the array B. LDB >= 2.
ALPHAR — ALPHAR is DOUBLE PRECISION array, dimension (2)
ALPHAI — ALPHAI is DOUBLE PRECISION array, dimension (2)
BETA — BETA is DOUBLE PRECISION array, dimension (2) (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may be zero.
CSL — CSL is DOUBLE PRECISION The cosine of the left rotation matrix.
SNL — SNL is DOUBLE PRECISION The sine of the left rotation matrix.
CSR — CSR is DOUBLE PRECISION The cosine of the right rotation matrix.
SNR — SNR is DOUBLE PRECISION The sine of the right rotation matrix.

🫶 yes, the source really does spell it "THe leading dimension" twice. we keep it. that's part of the lore now.

👑 Author 👑

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

💖 Contributors 💖

Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

🌟 subroutine dlahqr 🌟

subroutine dlahqr (logical WANTT, logical WANTZ, integer N, integer ILO,
    integer IHI, double precision, dimension( ldh, * ) H, integer LDH,
    double precision, dimension( * ) WR, double precision, dimension( * )
    WI, integer ILOZ, integer IHIZ, double precision, dimension( ldz, * ) Z,
    integer LDZ, integer INFO)

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

✨ Purpose ✨

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

so it's the little helper that DHSEQR calls in to do the cleanup work on rows and columns ILO to IHI. auxiliary but make it essential ✨

🎀 Parameters 🎀

WANTT — WANTT is LOGICAL
- = .TRUE. : the full Schur form T is required;
- = .FALSE.: only eigenvalues are required.
WANTZ — WANTZ is LOGICAL
- = .TRUE. : the matrix of Schur vectors Z is required;
- = .FALSE.: Schur vectors are not required.
N — N is INTEGER The order of the matrix H. N >= 0.
ILO — ILO is INTEGER
IHI — IHI is INTEGER It is assumed that H is already upper quasi-triangular in rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). DLAHQR 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 — okay bestie, H is where this part runs out of runway 🌸 the full 411 on H is rolling in with the next chunk, so hold that thought.
H — DOUBLE PRECISION array, dimension (LDH,N). On entry, it's the upper Hessenberg matrix H. On exit, if INFO is zero and WANTT is .TRUE., H comes back upper quasi-triangular in rows and columns ILO:IHI, and any 2-by-2 diagonal blocks are in standard form 💅. If INFO is zero and WANTT is .FALSE., the contents of H on exit are unspecified (aka: don't look at it, it's giving nothing, mind your business 🙈). The output state of H when INFO is nonzero? It's down below under the INFO description, bestie, keep scrolling 👇
LDH — INTEGER. It's the leading dimension of the array H. Requirement: LDH >= max(1,N).
WR — DOUBLE PRECISION array, dimension (N). That's it, that's the whole note 🎀
WI — DOUBLE PRECISION array, dimension (N). So WR and WI are the real and imaginary parts, respectively, of the computed eigenvalues ILO to IHI, and they're stored in the corresponding elements of WR and WI. If two eigenvalues get computed as a complex conjugate pair, they're stored in consecutive elements of WR and WI — like the i-th and (i+1)th — with WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i), and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). Math ate that fr fr 🔥
ILOZ — INTEGER.
IHIZ — INTEGER. These two specify the rows of Z that transformations get applied to, but only if WANTZ is .TRUE.. Constraints: 1 <= ILOZ <= ILO and IHI <= IHIZ <= N.
Z — DOUBLE PRECISION array, dimension (LDZ,N). If WANTZ is .TRUE., then on entry Z must contain the current matrix Z of transformations accumulated by DHSEQR, 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 literally not referenced. Not even glanced at. It's invisible bestie 👻
LDZ — INTEGER. The leading dimension of the array Z. Requirement: LDZ >= max(1,N).
INFO — INTEGER. The vibes check 👑:
- = 0: successful exit ✨
- > 0: If INFO = i, DLAHQR failed to compute all the eigenvalues ILO to IHI in a total of 30 iterations per eigenvalue; elements i+1:ihi of WR and WI 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.

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.

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 yes, that holds regardless of the value of WANTT 🫶

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

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 DLAHQR 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).

🦋 subroutine dlahr2

subroutine dlahr2 (integer N, integer K, integer NB, double precision,
    dimension( lda, * ) A, integer LDA, double precision, dimension( nb )
    TAU, double precision, dimension( ldt, nb ) T, integer LDT, double
    precision, dimension( ldy, nb ) Y, integer LDY)

DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and it hands back the auxiliary matrices you need to apply that transformation to the unreduced part of A. It's the little helper doing the heavy lifting, and honestly? iconic behavior 💖

Purpose:

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

This is an auxiliary routine called by DGEHRD.

Parameters

N — INTEGER. The order of the matrix A.
K — INTEGER. The offset for the reduction. Elements below the k-th subdiagonal in the first NB columns are reduced to zero. K < N.
NB — INTEGER. The number of columns to be reduced.
A — DOUBLE PRECISION 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, together 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 — INTEGER. The leading dimension of the array A. Requirement: LDA >= max(1,N).
TAU — DOUBLE PRECISION array, dimension (NB). The scalar factors of the elementary reflectors. See Further Details.
T — DOUBLE PRECISION array, dimension (LDT,NB). The upper triangular matrix T.
LDT — INTEGER. The leading dimension of the array T. Requirement: LDT >= NB.
Y — DOUBLE PRECISION array, dimension (LDY,NB). The n-by-nb matrix Y.
LDY — INTEGER. The leading dimension of the array Y. Requirement: 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**T

where tau is a real scalar, and v is a real 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**T) * (A - Y*V**T).

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 DLAHRD 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 DLAHRD routine. (This subroutine is not backward compatible with LAPACK-3.0's DLAHRD.) ⚠️ read that twice bestie, the old code will not recognize the new you 💀

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 dlaic1

subroutine dlaic1 (integer JOB, integer J, double precision, dimension( j )
    X, double precision SEST, double precision, dimension( j ) W, double
    precision GAMMA, double precision SESTPR, double precision S, double
    precision C)

DLAIC1 applies one step of incremental condition estimation. One step. Baby steps. But make it math 🧚

Purpose:

DLAIC1 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 DLAIC1 computes sestpr, s, c such that the vector

                [ s*x ]
         xhat = [  c  ]

is an approximate singular vector of

                [ L       0  ]
         Lhat = [ w**T 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]**T and sestpr**2 is an eigenpair of the system

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

where alpha = x**T*w.

Parameters

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

💫 Author 💫

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 `dlaln2` 🦋

subroutine dlaln2 (logical LTRANS, integer NA, integer NW, double precision
    SMIN, double precision CA, double precision, dimension( lda, * ) A,
    integer LDA, double precision D1, double precision D2, double precision,
    dimension( ldb, * ) B, integer LDB, double precision WR, double
    precision WI, double precision, dimension( ldx, * ) X, integer LDX,
    double precision SCALE, double precision XNORM, integer INFO)

okay so DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. It's giving tiny but mighty. 💅

🌸 Purpose 🌸

DLALN2 solves a system of the form

(ca A - w D ) X = s B

(ca A**T - w D) X = s B

with possible scaling ('s') and perturbation of A. (A**T means A-transpose, iykyk.) 🫶

A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA real diagonal matrix, w is a real or complex value, and X and B are NA x 1 matrices — real if w is real, complex if w is complex. And NA may be 1 or 2. That's it, that's the whole vibe. 💀

If w is complex, X and B are represented as NA x 2 matrices, the first column of each being the real part and the second being the imaginary part. Two columns, two moods. 🦋

's' is a scaling factor (<= 1), computed by DLALN2, which is so chosen that X can be computed without overflow. X is further scaled if necessary to assure that norm(ca A - w D)*norm(X) is less than overflow. Protective queen behavior. 👑

If both singular values of (ca A - w D) are less than SMIN, SMIN*identity will be used instead of (ca A - w D). If only one singular value is less than SMIN, one element of (ca A - w D) will be perturbed enough to make the smallest singular value roughly SMIN. If both singular values are at least SMIN, (ca A - w D) will not be perturbed. In any case, the perturbation will be at most some small multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values are computed by infinity-norm approximations, and thus will only be correct to a factor of 2 or so. Factor of 2?? close enough, we love an approximation era. 😭

Note: all input quantities are assumed to be smaller than overflow by a reasonable factor. (See BIGNUM.) fr fr, don't test it.

✨ Parameters ✨

LTRANS — LTRANS is LOGICAL. =.TRUE.: A-transpose will be used. =.FALSE.: A will be used (not transposed.) Pick your fighter. 💅
NA — NA is INTEGER. The size of the matrix A. It may (only) be 1 or 2. That's the whole range, bestie, no in-between. 🙈
NW — NW is INTEGER. 1 if 'w' is real, 2 if 'w' is complex. It may only be 1 or 2. 🎀
SMIN — SMIN is DOUBLE PRECISION. The desired lower bound on the singular values of A. This should be a safe distance away from underflow or overflow, say, between (underflow/machine precision) and (machine precision * overflow ). (See BIGNUM and ULP.) a safe distance, like emotionally. 🌸
CA — CA is DOUBLE PRECISION. The coefficient c, which A is multiplied by. 🧚
A — A is DOUBLE PRECISION array, dimension (LDA,NA). The NA x NA matrix A. The main character herself. 👑
LDA — LDA is INTEGER. The leading dimension of A. It must be at least NA. ⭐
D1 — D1 is DOUBLE PRECISION. The 1,1 element in the diagonal matrix D. 💜
D2 — D2 is DOUBLE PRECISION. The 2,2 element in the diagonal matrix D. Not used if NA=1. She gets the night off. 💤
B — B is DOUBLE PRECISION array, dimension (LDB,NW). The NA x NW matrix B (right-hand side). If NW=2 ('w' is complex), column 1 contains the real part of B and column 2 contains the imaginary part. 💖
LDB — LDB is INTEGER. The leading dimension of B. It must be at least NA. 🌈
WR — WR is DOUBLE PRECISION. The real part of the scalar 'w'. ✨
WI — WI is DOUBLE PRECISION. The imaginary part of the scalar 'w'. Not used if NW=1. 🎀
X — X is DOUBLE PRECISION array, dimension (LDX,NW). The NA x NW matrix X (unknowns), as computed by DLALN2. If NW=2 ('w' is complex), on exit, column 1 will contain the real part of X and column 2 will contain the imaginary part. the glow-up happens on exit. 💅
LDX — LDX is INTEGER. The leading dimension of X. It must be at least NA. 🦋
SCALE — SCALE is DOUBLE PRECISION. The scale factor that B must be multiplied by to insure that overflow does not occur when computing X. Thus, (ca A - w D) X will be SCALE*B, not B (ignoring perturbations of A.) It will be at most 1. at most 1, never overachieving, respect. 😌
XNORM — XNORM is DOUBLE PRECISION. The infinity-norm of X, when X is regarded as an NA x NW real matrix. 🍓
INFO — INFO is INTEGER. An error flag. It will be set to zero if no error occurs, a negative number if an argument is in error, or a positive number if ca A - w D had to be perturbed. The possible values are:
- = 0: No error occurred, and (ca A - w D) did not have to be perturbed. 🎉
- = 1: (ca A - w D) had to be perturbed to make its smallest (or only) singular value greater than SMIN. 😭

NOTE: In the interests of speed, this routine does not check the inputs for errors. not me finding out the hard way 💀

💫 Author 💫

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 `dlangt` 🦋

double precision function dlangt (character NORM, integer N, double
    precision, dimension( * ) DL, double precision, dimension( * ) D, double
    precision, dimension( * ) DU)

DLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. tridiagonal? we stan a three-line queen. 👑

🌸 Purpose 🌸

DLANGT returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real tridiagonal matrix A. ✨

✨ Returns ✨

DLANGT = ( 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. it's giving rebel, not a real norm, we still love it. 💅

✨ Parameters ✨

NORM — NORM is CHARACTER*1. Specifies the value to be returned in DLANGT as described above. Pick a letter, any letter (from the list). 🎀
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, DLANGT is set to zero. zero matrix, zero drama. 🌸
DL — DL is DOUBLE PRECISION array, dimension (N-1). The (n-1) sub-diagonal elements of A. 🦋
D — D is DOUBLE PRECISION array, dimension (N). The diagonal elements of A. The main character row. 👑
DU — DU is DOUBLE PRECISION 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.

🦋 `dlanhs` 🦋

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

DLANHS 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, rise up. 💖

🌸 Purpose 🌸

DLANHS 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 ✨

DLANHS = ( 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 DLANHS as described above. 🎀
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, DLANHS is set to zero. 🌸
A — A is DOUBLE PRECISION array, dimension (LDA,N). The n by n upper Hessenberg matrix A; the part of A below the first sub-diagonal is not referenced. if you're below the sub-diagonal, you're not invited, sorry bestie. 🙈
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. she only shows up for the 'I' era. 💜

💫 Author 💫

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 `dlansb` 🦋

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

DLANSB 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 AND banded?? iconic. 👑

🌸 Purpose 🌸

DLANSB 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 ✨

DLANSB

okay so if you asked for the biggest absolute value, the 1-norm, the Frobenius norm or the infinity norm, here's the little menu it picks from — for DLANSB it's giving:

DLANSB = ( 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 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). ✨ And bestie, note this: max(abs(A(i,j))) is not a consistent matrix norm. That's just facts, not a vibe.

Parameters

NORM — NORM is CHARACTER*1. It says which value gets handed back to you in DLANSB, as laid out above. 💅
UPLO — UPLO is CHARACTER*1. Tells it whether the upper or lower triangular part of the band matrix A is the one you're supplying.
- = 'U': upper triangular part is supplied
- = 'L': lower triangular part is supplied
N — N is INTEGER. The order of the matrix A. N >= 0. And when N = 0, DLANSB is set to zero. (Zero. Gone. Nothing. 🫶)
K — K is INTEGER. The number of super-diagonals or sub-diagonals of the band matrix A. K >= 0.
AB — AB is DOUBLE PRECISION array, dimension (LDAB,N). The upper or lower triangle of the symmetric 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).
LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= K+1.
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 dlansp (character NORM, character UPLO, integer N, double precision, dimension( * ) AP, double precision, dimension( * ) WORK)`

DLANSP 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. Like your suitcase. 🧳

Purpose:

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

Returns

DLANSP =

DLANSP = ( 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 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). And again, no cap: max(abs(A(i,j))) is not a consistent matrix norm. 💀

Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in DLANSP as described above. 👑
UPLO — UPLO is CHARACTER*1. Says 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 — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, DLANSP is set to zero. (iconic zero behaviour fr fr)
AP — AP is DOUBLE PRECISION 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 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.
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 dlantb (character NORM, character UPLO, character DIAG, integer N, integer K, double precision, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) WORK)`

DLANTB 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 band edition, bestie. 💜

Purpose:

DLANTB  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

DLANTB =

DLANTB = ( 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 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). And once more for the people in the back: max(abs(A(i,j))) is not a consistent matrix norm. 😭

Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in DLANTB as described above.
UPLO — UPLO is CHARACTER*1. Tells it whether the matrix A is upper or lower triangular.
- = 'U': upper triangular
- = 'L': lower triangular
DIAG — DIAG is CHARACTER*1. Says whether or not the matrix A is unit triangular.
- = 'N': non-unit triangular
- = 'U': unit triangular
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, DLANTB is set to zero. ✨
K — K is 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 — AB is DOUBLE PRECISION 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).

Heads up: 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. So don't stress about writing them, they're just ~vibes~ assumed to be 1. 🎀

LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= K+1.
WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced. (Note the difference, bestie — it's only 'I' here, not the whole '1'/'O' crew. iykyk 🫶)

Author

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

`double precision function dlantp (character NORM, character UPLO, character DIAG, integer N, double precision, dimension( * ) AP, double precision, dimension( * ) WORK)`

DLANTP 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 AND packed. The crossover episode nobody asked for but we got anyway. 🌈

Purpose:

DLANTP  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

DLANTP =

DLANTP = ( 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 is the one norm of a matrix (maximum column sum), normI is the infinity norm (maximum row sum), and normF is the Frobenius norm (square root of sum of squares). And yes, one more time: max(abs(A(i,j))) is not a consistent matrix norm. 🍓

Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in DLANTP as described above.
UPLO — UPLO is CHARACTER*1. Tells it whether the matrix A is upper or lower triangular.
- = 'U': upper triangular
- = 'L': lower triangular
DIAG — DIAG is CHARACTER*1. Says whether or not the matrix A is unit triangular.
- = 'N': non-unit triangular
- = 'U': unit triangular
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, DLANTP is set to zero. 💅
AP — AP is DOUBLE PRECISION 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.

And same energy as before: 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. They're literally just placeholders, don't even look at them. 🙈

WORK — WORK is 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.

🦋 `dlantr` 💜

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

okay so bestie, DLANTR hands you back the 1-norm, or the Frobenius norm, or the infinity norm, or just the element with the largest absolute value of a trapezoidal or triangular matrix. four-for-one deal, we love to see it ✨

🌸 Purpose 🌸

DLANTR 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 ✨

DLANTR = ( 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, no cap 🙈

🎀 Parameters 🎀

NORM — NORM is CHARACTER*1. Tells you which value gets returned in DLANTR, exactly as described up there ⬆️
UPLO — UPLO is CHARACTER*1. Says whether the matrix A is upper or lower trapezoidal.
- = 'U': Upper trapezoidal
- = 'L': Lower trapezoidal
- ngl, A is triangular instead of trapezoidal if M = N 💅
DIAG — DIAG is CHARACTER*1. Says whether or not the matrix A has a 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, DLANTR 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, DLANTR is set to zero 🫶
A — A is DOUBLE PRECISION 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. she only shows up for the infinity norm, main character energy 🌈

👑 Author 👑

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

🔥 `dlanv2` 🔥

subroutine dlanv2 (double precision A, double precision B, double precision
    C, double precision D, double precision RT1R, double precision RT1I,
    double precision RT2R, double precision RT2I, double precision CS,
    double precision SN)

DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form. tiny matrix, big brain energy 🧚

🌟 Purpose 🌟

DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form:

     [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
     [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]

where either

CC = 0 so that AA and DD are real eigenvalues of the matrix, or
AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex conjugate eigenvalues.

💖 Parameters 💖

A — A is DOUBLE PRECISION
B — B is DOUBLE PRECISION
C — C is DOUBLE PRECISION
D — D is DOUBLE PRECISION. On entry, the elements of the input matrix. On exit, they are overwritten by the elements of the standardised Schur form. bye original values, it was real 😭
RT1R — RT1R is DOUBLE PRECISION
RT1I — RT1I is DOUBLE PRECISION
RT2R — RT2R is DOUBLE PRECISION
RT2I — RT2I is DOUBLE PRECISION. The real and imaginary parts of the eigenvalues. If the eigenvalues are a complex conjugate pair, RT1I > 0.
CS — CS is DOUBLE PRECISION
SN — SN is DOUBLE PRECISION. Parameters of the rotation matrix. spinny spinny 💫

👑 Author 👑

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

🦋 Further Details: 🦋

Modified by V. Sima, Research Institute for Informatics, Bucharest, Romania, to reduce the risk of cancellation errors, when computing real eigenvalues, and to ensure, if possible, that abs(RT1R) >= abs(RT2R). protective queen behaviour, we stan 🫶

🌈 `dlapll` 🌈

subroutine dlapll (integer N, double precision, dimension( * ) X, integer
    INCX, double precision, dimension( * ) Y, integer INCY, double precision
    SSMIN)

DLAPLL measures the linear dependence of two vectors. are y'all actually different or just vibing the same? vibe check incoming ⭐

🌸 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. literally a whole detective arc for two vectors 💀

🎀 Parameters 🎀

N — N is INTEGER. The length of the vectors X and Y.
X — X is DOUBLE PRECISION array, dimension (1+(N-1)*INCX). On entry, X contains the N-vector X. On exit, X is overwritten. sorry babe, she's changed now 🙈
INCX — INCX is INTEGER. The increment between successive elements of X. INCX > 0.
Y — Y is DOUBLE PRECISION array, dimension (1+(N-1)*INCY). On entry, Y contains the N-vector Y. On exit, Y is overwritten. same fate fr fr 😭
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.

💅 `dlapmr` 💅

subroutine dlapmr (logical FORWRD, integer M, integer N, double precision,
    dimension( ldx, * ) X, integer LDX, integer, dimension( * ) K)

DLAPMR rearranges rows of a matrix as specified by a permutation vector. row shuffle era ✨

🌟 Purpose 🌟

DLAPMR 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 DOUBLE PRECISION 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. she borrows your stuff and gives it back, that's the respectful icon behaviour 🫶

👑 Author 👑

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

🦋 `dlapmt` 🦋

subroutine dlapmt (logical FORWRD, integer M, integer N, double precision,
    dimension( ldx, * ) X, integer LDX, integer, dimension( * ) K)

DLAPMT performs a forward or backward permutation of the columns of a matrix. DLAPMR's sibling but she does columns, not rows, and yes they're both iconic in their own lane 👑

🌈 Purpose 🌈

DLAPMT 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 DOUBLE PRECISION 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. again, borrowed and returned pristine, queen shit 💅

👑 Author 👑

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

🌸 `dlaqp2` 🌸

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

DLAQP2 computes a QR factorization with column pivoting of the matrix block. yes bestie, pivoting, the whole bit 🔥

🌟 Purpose 🌟

(this is where we pick up in the next part, don't you dare stop reading here) 💖

🧮 `dlaqp2` — the little pivoting QR bestie

What it does: DLAQP2 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. Like, it gets shuffled around, but it does not get the glow-up. 🙈

✨ Parameters ✨

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.
OFFSET — OFFSET is INTEGER. The number of rows of the matrix A that must be pivoted but no factorized. OFFSET >= 0. (Yes bestie, the original really does say "no factorized" and I am not touching it. 💅)
A — A is DOUBLE PRECISION 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 — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
JPVT — JPVT is 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 — TAU is DOUBLE PRECISION array, dimension (min(M,N)). The scalar factors of the elementary reflectors.
VN1 — VN1 is DOUBLE PRECISION array, dimension (N). The vector with the partial column norms.
VN2 — VN2 is DOUBLE PRECISION array, dimension (N). The vector with the exact column norms.
WORK — WORK is DOUBLE PRECISION 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 `dlaqps` 🌈

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

DLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. BLAS-3, bestie — we love a queen who batches her work. ✨

🦋 Purpose:

DLAQPS computes a step of QR factorization with column pivoting of a real 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.

Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. (Again. It's giving "invited to the party but not on the list." 🎀)

✨ Parameters ✨

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
OFFSET — OFFSET is INTEGER. The number of rows of A that have been factorized in previous steps.
NB — NB is INTEGER. The number of columns to factorize.
KB — KB is INTEGER. The number of columns actually factorized.
A — A is DOUBLE PRECISION 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 — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
JPVT — JPVT is INTEGER array, dimension (N). JPVT(I) = K <==> Column K of the full matrix A has been permuted into position I in AP.
TAU — TAU is DOUBLE PRECISION array, dimension (KB). The scalar factors of the elementary reflectors.
VN1 — VN1 is DOUBLE PRECISION array, dimension (N). The vector with the partial column norms.
VN2 — VN2 is DOUBLE PRECISION array, dimension (N). The vector with the exact column norms.
AUXV — AUXV is DOUBLE PRECISION array, dimension (NB). Auxiliary vector.
F — F is DOUBLE PRECISION array, dimension (LDF,NB). Matrix F**T = L*Y**T*A.
LDF — LDF is 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 `dlaqr0` 💜

subroutine dlaqr0 (logical WANTT, logical WANTZ, integer N, integer ILO,
    integer IHI, double precision, dimension( ldh, * ) H, integer LDH,
    double precision, dimension( * ) WR, double precision, dimension( * )
    WI, integer ILOZ, integer IHIZ, double precision, dimension( ldz, * ) Z,
    integer LDZ, double precision, dimension( * ) WORK, integer LWORK,
    integer INFO)

DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. EIGENVALUES AND SCHUR FORM, the full main-character arc. 👑

🌸 Purpose:

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

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

✨ Parameters ✨

WANTT — WANTT is LOGICAL.
- = .TRUE. : the full Schur form T is required;
- = .FALSE.: only eigenvalues are required.
WANTZ — WANTZ is LOGICAL.
- = .TRUE. : the matrix of Schur vectors Z is required;
- = .FALSE.: Schur vectors are not required.
N — N is INTEGER. The order of the matrix H. N >= 0.
ILO — ILO is INTEGER.
IHI — IHI is 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 DGEBAL, and then passed to DGEHRD when the matrix output by DGEBAL 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 — H is DOUBLE PRECISION 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 quasi-triangular matrix T from the Schur decomposition (the Schur form); 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT 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 is INTEGER. The leading dimension of the array H. LDH >= max(1,N).
WR — WR is DOUBLE PRECISION array, dimension (IHI).
WI — WI is DOUBLE PRECISION array, dimension (IHI).
- The real and imaginary parts, respectively, of the computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI) and WI(ILO:IHI). If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
ILOZ — ILOZ is INTEGER.
IHIZ — 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 — Z is DOUBLE PRECISION 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.)
LDZ
LDZ — INTEGER. This is the leading dimension of the array Z. If WANTZ is .TRUE., then LDZ >= MAX(1,IHIZ). Otherwise? Chill, LDZ >= 1 is all you need. 💅
WORK — DOUBLE PRECISION array, dimension LWORK. On exit, if LWORK = -1, then WORK(1) hands you back an estimate of the optimal value for LWORK. ✨
LWORK — INTEGER. The dimension of the array WORK. Okay so LWORK >= max(1,N) is sufficient, but bestie, LWORK typically as large as 6*N may be required for optimal performance 😤. A workspace query to determine the optimal workspace size is recommended — like, do it.

If LWORK = -1, then DLAQR0 does a workspace query. In this case, DLAQR0 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. (not a single one of them gets touched, iconic)

INFO — INTEGER
- = 0: successful exit 🎉
- > 0: if INFO = i, DLAQR0 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 an orthogonal matrix. The final value of H is upper Hessenberg and quasi-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 orthogonal matrix in (*) (regardless of the value of WANTT.) no cap, it really does not care what WANTT is

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

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

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 subroutine dlaqr1 (integer N, double precision, dimension( ldh, * ) H, integer LDH, double precision SR1, double precision SI1, double precision SR2, double precision SI2, double precision, dimension( * ) V)

DLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. 🧚

✨ Purpose

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

(*)  K = (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)

scaling to avoid overflows and most underflows. It is assumed that either

    1) sr1 = sr2 and si1 = -si2
or
    2) si1 = si2 = 0.

This is useful for starting double implicit shift bulges in the QR algorithm. 🔥

💜 Parameters

N — INTEGER. Order of the matrix H. N must be either 2 or 3. (just those two, bestie, no in-between)
H — DOUBLE PRECISION 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
SR1 — DOUBLE PRECISION
SI1 — DOUBLE PRECISION
SR2 — DOUBLE PRECISION
SI2 — DOUBLE PRECISION. The shifts in (*). ✨
V — DOUBLE PRECISION array, dimension (N). A scalar multiple of the first column of the matrix K in (*).

👑 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 dlaqr2 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT, integer NW, double precision, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ, double precision, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, double precision, dimension( * ) SR, double precision, dimension( * ) SI, double precision, dimension( ldv, * ) V, integer LDV, integer NH, double precision, dimension( ldt, * ) T, integer LDT, integer NV, double precision, dimension( ldwv, * ) WV, integer LDWV, double precision, dimension( * ) WORK, integer LWORK)

DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). 😤✨

✨ Purpose

DLAQR2 is identical to DLAQR3 except that it avoids recursion by calling DLAHQR instead of DLAQR4.

Aggressive early deflation:

This subroutine accepts as input an upper Hessenberg matrix H and performs an orthogonal 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 orthogonal similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. manifesting zero subdiagonal entries fr fr 🙏

💜 Parameters

WANTT — LOGICAL. If .TRUE., then the Hessenberg matrix H is fully updated so that the quasi-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 orthogonal matrix Z is updated so so that the orthogonal 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 orthogonal 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. (yes, "without a check" — it is simply trusting you, don't betray it)
NW — INTEGER. Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
H — DOUBLE PRECISION 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 an orthogonal 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
IHIZ — INTEGER. Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
Z — DOUBLE PRECISION array, dimension (LDZ,N). IF WANTZ is .TRUE., then on output, the orthogonal 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. 🎀
SR — DOUBLE PRECISION array, dimension (KBOT)
SI — DOUBLE PRECISION array, dimension (KBOT). On output, the real and imaginary parts of approximate eigenvalues that may be used for shifts are stored in SR(KBOT-ND-NS+1) through SR(KBOT-ND) and SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively. The real and imaginary parts of converged eigenvalues are stored in SR(KBOT-ND+1) through SR(KBOT) and SI(KBOT-ND+1) through SI(KBOT), respectively.
V — DOUBLE PRECISION 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 — DOUBLE PRECISION array, dimension (LDT,NW)
LDT
LDT — bestie, LDT is INTEGER, and it's the leading dimension of T, just as declared in the calling subroutine. So like: NW <= LDT 💖
NV — NV is INTEGER. It's the number of rows of work array WV that are actually available for workspace. The rule is NV >= NW. slay ✨
WV — WV is a DOUBLE PRECISION array, dimension (LDWV,NW) 🦋
LDWV — LDWV is INTEGER. It's the leading dimension of W just as declared in the calling subroutine. And the constraint as written is NW <= LDV 👑 (yes, exactly like that, we're relaying it character for character, no editing bestie)
WORK — WORK is a DOUBLE PRECISION 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. it's giving "here's your little hint, go off" 🎀
LWORK — LWORK is INTEGER. It's the dimension of the work array WORK. Here's the tea: LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK. 💅

And now the plot twist, the part you actually need to know: if LWORK = -1, then a workspace query is assumed. DLAQR2 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. That estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. And — huge — neither H nor Z are accessed. so nothing gets touched, it's just vibes and math 🌸

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 `dlaqr3` 💅

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

DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).

Purpose: 🎯

Aggressive early deflation:

DLAQR3 accepts as input an upper Hessenberg matrix H and performs an orthogonal 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 orthogonal similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. ✨ (manifesting zero subdiagonal entries for you bestie 🫶)

Parameters 🎀

WANTT — WANTT is LOGICAL. If .TRUE., then the Hessenberg matrix H is fully updated so that the quasi-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. lowkey minimalist mode 💅
WANTZ — WANTZ is LOGICAL. If .TRUE., then the orthogonal matrix Z is updated so so that the orthogonal Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced. (and yes, the original really says "so so" — iconic, we're keeping it 💖)
N — N is INTEGER. The order of the matrix H and (if WANTZ is .TRUE.) the order of the orthogonal matrix Z. 👑
KTOP — KTOP is 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 — KBOT is 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. (no check bestie, we're trusting you 🙈)
NW — NW is INTEGER. Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
H — H is a DOUBLE PRECISION 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 an orthogonal similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries. 🌸
LDH — LDH is INTEGER. Leading dimension of H just as declared in the calling subroutine. N <= LDH
ILOZ — ILOZ is INTEGER 💀 (that's it, that's the whole description, we love a mysterious queen)
IHIZ — IHIZ is INTEGER. Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
Z — Z is a DOUBLE PRECISION array, dimension (LDZ,N). IF WANTZ is .TRUE., then on output, the orthogonal 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 — LDZ is INTEGER. The leading dimension of Z just as declared in the calling subroutine. 1 <= LDZ.
NS — NS is INTEGER. The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine.
ND — ND is INTEGER. The number of converged eigenvalues uncovered by this subroutine. 💖
SR — SR is a DOUBLE PRECISION array, dimension (KBOT)
SI — SI is a DOUBLE PRECISION array, dimension (KBOT). On output, the real and imaginary parts of approximate eigenvalues that may be used for shifts are stored in SR(KBOT-ND-NS+1) through SR(KBOT-ND) and SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively. The real and imaginary parts of converged eigenvalues are stored in SR(KBOT-ND+1) through SR(KBOT) and SI(KBOT-ND+1) through SI(KBOT), respectively. ✨
V — V is a DOUBLE PRECISION array, dimension (LDV,NW). An NW-by-NW work array. 🎀
LDV — LDV is INTEGER. The leading dimension of V just as declared in the calling subroutine. NW <= LDV
NH — NH is INTEGER. The number of columns of T. NH >= NW.
T — T is a DOUBLE PRECISION array, dimension (LDT,NW)
LDT — LDT is INTEGER. The leading dimension of T just as declared in the calling subroutine. NW <= LDT
NV — NV is INTEGER. The number of rows of work array WV available for workspace. NV >= NW.
WV — WV is a DOUBLE PRECISION array, dimension (LDWV,NW)
LDWV — LDWV is INTEGER. The leading dimension of W just as declared in the calling subroutine. NW <= LDV 👑 (again, exactly as written, we don't edit the queen)
WORK — WORK is a DOUBLE PRECISION 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 — LWORK is INTEGER. The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK. 💅

And the iconic plot twist again: if LWORK = -1, then a workspace query is assumed; DLAQR3 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. nothing touched, nothing harmed, just a little estimate bestie 🌸

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 `dlaqr4` 🦋

subroutine dlaqr4 (logical WANTT, logical WANTZ, integer N, integer ILO,
    integer IHI, double precision, dimension( ldh, * ) H, integer LDH,
    double precision, dimension( * ) WR, double precision, dimension( * )
    WI, integer ILOZ, integer IHIZ, double precision, dimension( ldz, * ) Z,
    integer LDZ, double precision, dimension( * ) WORK, integer LWORK,
    integer INFO)

DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.

Purpose: 🎯

DLAQR4 implements one level of recursion for DLAQR0. It is a complete implementation of the small bulge multi-shift QR algorithm. It may be called by DLAQR0 and, for large enough deflation window size, it may be called by DLAQR3. This subroutine is identical to DLAQR0 except that it calls DLAQR2 instead of DLAQR3. 👑

DLAQR4 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur form), and Z is the orthogonal matrix of Schur vectors. ✨

Parameters 🎀

WANTT — WANTT is LOGICAL. = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.
WANTZ — WANTZ is LOGICAL. = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.
N — N is INTEGER. The order of the matrix H. N >= 0. 💖
ILO — ILO is INTEGER
IHI —
IHI is INTEGER — and the vibe here is that we're assuming H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N, and, if ILO > 1, that H(ILO,ILO-1) is zero. ILO and IHI are normally set by a previous call to DGEBAL, and then passed to DGEHRD when the matrix output by DGEBAL gets reduced to Hessenberg form. If you're not doing that whole pipeline, just set them to 1 and N respectively. Also: if N > 0, then 1 <= ILO <= IHI <= N. And if N = 0, then ILO = 1 and IHI = 0. periodt ✨
H — H is DOUBLE PRECISION array, dimension (LDH,N). On entry, she's the upper Hessenberg matrix H 💅. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper quasi-triangular matrix T from the Schur decomposition (aka the Schur form, iconic), and the 2-by-2 diagonal blocks (those are the complex conjugate pairs of eigenvalues, bestie) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is .FALSE., then ngl the contents of H are just unspecified on exit. (Wanna know the output value of H when INFO > 0? It's down under the description of INFO, keep scrolling 🙈.)

Heads up: 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 is INTEGER — the leading dimension of the array H. LDH >= max(1,N).
WR is DOUBLE PRECISION array, dimension (IHI).
WI is DOUBLE PRECISION array, dimension (IHI). The real and imaginary parts, respectively, of the computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI) and WI(ILO:IHI). If two eigenvalues get computed as a complex conjugate pair, they're stored in consecutive elements of WR and WI — say the i-th and (i+1)th — with WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). math girlies stay winning 🧮💖
ILOZ is INTEGER.
IHIZ is INTEGER — together these two specify the rows of Z that transformations have to be applied to if WANTZ is .TRUE.. Constraints, no cap: 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
Z — Z is DOUBLE PRECISION array, dimension (LDZ,IHI). If WANTZ is .FALSE., then Z is not referenced, she's off the clock 🫶. 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 lives under the description of INFO below, iykyk.)
LDZ is INTEGER — the leading dimension of the array Z. If WANTZ is .TRUE. then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
WORK — WORK is DOUBLE PRECISION array, dimension LWORK. On exit, if LWORK = -1, WORK(1) returns an estimate of the optimal value for LWORK.
LWORK is 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, so a workspace query to determine the optimal workspace size is recommended, fr fr.

If LWORK = -1, then DLAQR4 does a workspace query. In this case DLAQR4 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, DLAQR4 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, bestie, 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 orthogonal 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 orthogonal 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 dlaqr5

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

DLAQR5 performs a single small-bulge multi-shift QR sweep.

✨ Purpose

DLAQR5, called by DLAQR0, performs a single small-bulge multi-shift QR sweep. that's it, that's the whole job description, and she does it flawlessly 💅

🎀 Parameters

WANTT is LOGICAL — WANTT = .true. if the quasi-triangular Schur factor is being computed. WANTT is set to .false. otherwise.
WANTZ is LOGICAL — WANTZ = .true. if the orthogonal Schur factor is being computed. WANTZ is set to .false. otherwise.
KACC22 is INTEGER with value 0, 1, or 2. It specifies the computation mode of far-from-diagonal orthogonal updates. 👀
- = 0: DLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries.
- = 1: DLAQR5 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. rip to that feature, she's retired 💀
N is INTEGER — N is the order of the Hessenberg matrix H upon which this subroutine operates.
KTOP is INTEGER.
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. no odd numbers allowed, she's picky like that ✨
SR is DOUBLE PRECISION array, dimension (NSHFTS).
SI is DOUBLE PRECISION array, dimension (NSHFTS). SR contains the real parts and SI contains the imaginary parts of the NSHFTS shifts of origin that define the multi-shift QR sweep. On output SR and SI may be reordered.
H is DOUBLE PRECISION 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.
LDH is INTEGER — LDH is the leading dimension of H just as declared in the calling procedure. LDH >= MAX(1,N).
ILOZ is INTEGER.
IHIZ is INTEGER — these specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
Z is DOUBLE PRECISION array, dimension (LDZ,IHIZ). If WANTZ = .TRUE., then the QR Sweep orthogonal similarity transformation is accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ = .FALSE., then Z is unreferenced, she said not today 💅
LDZ is INTEGER — LDA is the leading dimension of Z just as declared in the calling procedure. LDZ >= N.
V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2).
LDV
LDV — an INTEGER, and she's the leading dimension of V, exactly as you declared it back in the calling procedure. The rule is LDV >= 3, no exceptions bestie ✨
U — a DOUBLE PRECISION array, dimension (LDU,2*NSHFTS). she's a whole 2D situation 🦋
LDU — an INTEGER. LDU is the leading dimension of U, just as you declared it in the calling subroutine. And it needs LDU >= 2*NSHFTS. periodt 💅
NV — an INTEGER. NV is the number of rows in WV available for workspace, and NV >= 1. give her at least one row, she's not asking for much 🥺
WV — a DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS).
LDWV — an INTEGER. LDWV is the leading dimension of WV as declared in the calling subroutine, and LDWV >= NV. the way that these leading dimensions just keep coming 😭
NH — an INTEGER. NH is the number of columns in the array WH available for workspace, with NH >= 1.
WH — a DOUBLE PRECISION array, dimension (LDWH,NH).
LDWH — an INTEGER. Leading dimension of WH, just as declared in the calling procedure, and LDWH >= 2*NSHFTS. last one of the workspace girlies, we made it 🎉

🌈 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 dlaqsb` 🎀

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

okay so DLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. she does not compute the factors herself, she just applies them. iconic delegation 💅

Purpose ✨

DLAQSB equilibrates a symmetric 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 DOUBLE PRECISION 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**T*U or A = L*L**T of the band matrix A, in the same storage format as A. so she swaps the storage in place, rude but efficient 💀

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.

🌸 `subroutine dlaqsp` 🌸

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

DLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. packed storage bestie, she's living her compact little life 🎀

Purpose ✨

DLAQSP equilibrates a symmetric 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.
AP — AP is DOUBLE PRECISION array, dimension (N*(N+1)/2).

On entry, 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.

On exit, the equilibrated matrix: diag(S) * A * diag(S), in the same storage format as A. yes she rewrites herself in place, main character energy 👑

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.

🧚 `subroutine dlaqtr` 🧚

subroutine dlaqtr (logical LTRAN, logical LREAL, integer N, double
    precision, dimension( ldt, * ) T, integer LDT, double precision,
    dimension( * ) B, double precision W, double precision SCALE, double
    precision, dimension( * ) X, double precision, dimension( * ) WORK,
    integer INFO)

DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of special form, in real arithmetic. complex math done entirely in real arithmetic?? it's giving chaos but make it elegant, honestly obsessed 😭✨

Purpose ✨

DLAQTR solves the real quasi-triangular system

             op(T)*p = scale*c,               if LREAL = .TRUE.

or the complex quasi-triangular systems

           op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.

in real arithmetic, where T is upper quasi-triangular.

If LREAL = .FALSE., then the first diagonal block of T must be 1 by 1, B is the specially structured matrix

               B = [ b(1) b(2) ... b(n) ]
                   [       w            ]
                   [           w        ]
                   [              .     ]
                   [                 w  ]

op(A) = A or A**T, A**T denotes the transpose of matrix A.

On input, X = [ c ].  On output, X = [ p ].
              [ d ]                  [ q ]

This subroutine is designed for the condition number estimation in routine DTRSNA. so she's a supporting character but she ate 🫶

Parameters 🦋

LTRAN — LTRAN is LOGICAL. On entry, LTRAN specifies the option of conjugate transpose:
- = .FALSE., op(T+i*B) = T+i*B,
- = .TRUE., op(T+i*B) = (T+i*B)**T.
LREAL — LREAL is LOGICAL. On entry, LREAL specifies the input matrix structure:
- = .FALSE., the input is complex
- = .TRUE., the input is real
N — N is INTEGER. On entry, N specifies the order of T+i*B. N >= 0.
T — T is DOUBLE PRECISION array, dimension (LDT,N). On entry, T contains a matrix in Schur canonical form. If LREAL = .FALSE., then the first diagonal block of T must be 1 by 1.
LDT — LDT is INTEGER. The leading dimension of the matrix T. LDT >= max(1,N).
B — B is DOUBLE PRECISION array, dimension (N). On entry, B contains the elements to form the matrix B as described above. If LREAL = .TRUE., B is not referenced. she's literally not even looking at you bestie 🙈
W — W is DOUBLE PRECISION. On entry, W is the diagonal element of the matrix B. If LREAL = .TRUE., W is not referenced.
SCALE — SCALE is DOUBLE PRECISION. On exit, SCALE is the scale factor.
X — X is DOUBLE PRECISION array, dimension (2*N). On entry, X contains the right hand side of the system. On exit, X is overwritten by the solution.
WORK — WORK is DOUBLE PRECISION array, dimension (N)
INFO

INFO is INTEGER — and bestie, on the way out it gets set to one of these:

0: successful exit. slay, we move. 💅
1: the some diagonal 1 by 1 block has been perturbed by a small number SMIN to keep nonsingularity. (iykyk — a lil nudge so the thing doesn't fall apart)
2: the some diagonal 2 by 2 block has been perturbed by a small number in DLALN2 to keep nonsingularity.

NOTE: In the interests of speed, this routine does not check the inputs for errors. literally no vibe check, no seatbelt, she just GOES 🏃‍♀️💨 so pls pls pls know what you're doing.

Author 🫶

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

`subroutine dlar1v`

subroutine dlar1v (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, double
    precision, 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)

DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. iconic behaviour honestly 🌸

Purpose:

DLAR1V computes the (scaled) r-th column of the inverse of the sumbmatrix 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:

(a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
(b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
(c) 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.
(d) 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.

Parameters

N — N is INTEGER. The order of the matrix L D L**T. 📐
B1 — B1 is INTEGER. First index of the submatrix of L D L**T.
BN — BN is INTEGER. Last index of the submatrix of L D L**T.
LAMBDA — LAMBDA is 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. (emphasis on good, bestie, she's picky ✨)
L — L is DOUBLE PRECISION array, dimension (N-1). The (n-1) subdiagonal elements of the unit bidiagonal matrix L, in elements 1 to N-1.
D — D is DOUBLE PRECISION array, dimension (N). The n diagonal elements of the diagonal matrix D.
LD — LD is DOUBLE PRECISION array, dimension (N-1). The n-1 elements L(i)*D(i).
LLD — LLD is DOUBLE PRECISION array, dimension (N-1). The n-1 elements L(i)*L(i)*D(i).
PIVMIN — PIVMIN is DOUBLE PRECISION. The minimum pivot in the Sturm sequence.
GAPTOL — GAPTOL is DOUBLE PRECISION. Tolerance that indicates when eigenvector entries are negligible w.r.t. their contribution to the residual.
Z — Z is DOUBLE PRECISION array, dimension (N). On input, all entries of Z must be set to 0. On output, Z contains the (scaled) r-th column of the inverse. The scaling is such that Z(R) equals 1. 🎀
WANTNC — WANTNC is LOGICAL. Specifies whether NEGCNT has to be computed.
NEGCNT — NEGCNT is 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 — ZTZ is DOUBLE PRECISION. The square of the 2-norm of Z.
MINGMA — MINGMA is DOUBLE PRECISION. The reciprocal of the largest (in magnitude) diagonal element of the inverse of L D L**T - sigma I.
R — R is 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 is 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 ).
NRMINV — NRMINV is DOUBLE PRECISION. NRMINV = 1/SQRT( ZTZ )
RESID — RESID is DOUBLE PRECISION. The residual of the FP vector. RESID = ABS( MINGMA )/SQRT( ZTZ )
RQCORR — RQCORR is DOUBLE PRECISION. The Rayleigh Quotient correction to LAMBDA. RQCORR = MINGMA*TMP
WORK — WORK is 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 dlar2v`

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

DLAR2V applies a vector of plane rotations with real cosines and real sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. both sides?? she's ambidextrous 💅✨

Purpose:

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

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

Parameters

N — N is INTEGER. The number of plane rotations to be applied.
X — X is DOUBLE PRECISION array, dimension (1+(N-1)*INCX). The vector x.
Y — Y is DOUBLE PRECISION array, dimension (1+(N-1)*INCX). The vector y.
Z — Z is DOUBLE PRECISION array, dimension (1+(N-1)*INCX). The vector z.
INCX — INCX is INTEGER. The increment between elements of X, Y and Z. INCX > 0.
C — C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC). The cosines of the plane rotations.
S — S is DOUBLE PRECISION 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 dlarf`

subroutine dlarf (character SIDE, integer M, integer N, double precision,
    dimension( * ) V, integer INCV, double precision TAU, double precision,
    dimension( ldc, * ) C, integer LDC, double precision, dimension( * )
    WORK)

DLARF applies an elementary reflector to a general rectangular matrix. 🪞✨

Purpose:

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

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

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

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

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 DOUBLE PRECISION 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 — INCV is INTEGER. The increment between elements of v. INCV <> 0.
TAU — TAU is DOUBLE PRECISION. The value tau in the representation of H.
C — C is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dlarfb`

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

DLARFB applies a block reflector or its transpose to a general rectangular matrix. 🔥 block reflector?? she said BLOCK, we love a queen who does things in bulk.

Purpose:

<span>okay so hold that thought bestie, the purpose section is pulling up right now… 🦋</span>

…so like, DLARFB takes a real block reflector H — or its transpose H**T — and applies it to a real m-by-n matrix C, from either the left or the right. That's the whole vibe. 💅

✨ Parameters ✨

SIDE — SIDE is CHARACTER*1
- = 'L': apply H or H**T from the Left 🫲
- = 'R': apply H or H**T from the Right 🫱
TRANS — TRANS is CHARACTER*1
- = 'N': apply H (No transpose)
- = 'T': apply H**T (Transpose)
DIRECT — DIRECT is CHARACTER*1. Tells you how H is formed from a product of elementary reflectors, bestie:
- = 'F': H = H(1) H(2) . . . H(k) (Forward)
- = 'B': H = H(k) . . . H(2) H(1) (Backward)
STOREV — STOREV is CHARACTER*1. Tells you how the vectors that define the elementary reflectors are stored:
- = 'C': Columnwise
- = 'R': Rowwise
M — M is INTEGER. The number of rows of the matrix C. 📏
N — N is INTEGER. The number of columns of the matrix C. 📐
K — K is INTEGER. The order of the matrix T (aka the number of elementary reflectors whose product defines the block reflector). The rules, no cap:
- If SIDE = 'L', then M >= K >= 0;
- if SIDE = 'R', then N >= K >= 0.
V — V is DOUBLE PRECISION array, dimension

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

The matrix V. See Further Details, bestie, there's a whole diagram. 🎀

LDV — LDV is INTEGER. The leading dimension of the array V. The rules:
- If STOREV = 'C' and SIDE = 'L', then LDV >= max(1,M);
- if STOREV = 'C' and SIDE = 'R', then LDV >= max(1,N);
- if STOREV = 'R', then LDV >= K.
T — T is DOUBLE PRECISION array, dimension (LDT,K). The 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.
C — C is DOUBLE PRECISION array, dimension (LDC,N).
- On entry: the m-by-n matrix C.
- On exit: C is overwritten by H*C or H**T*C or C*H or C*H**T. She's getting replaced, periodt. 💀
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
WORK — WORK is DOUBLE PRECISION array, dimension (LDWORK,K).
LDWORK — LDWORK is INTEGER. The leading dimension of the array WORK. The rules:
- If SIDE = 'L', then LDWORK >= max(1,N);
- if SIDE = 'R', then 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. Lowkey iconic. ✨

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 )

💖 subroutine dlarfb_gett 💖

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

DLARFB_GETT 🧚

🎀 Purpose: 🎀

DLARFB_GETT applies a real Householder block reflector H from the left to a real (K+M)-by-N 'triangular-pentagonal' matrix composed of two block matrices: an upper trapezoidal K-by-N matrix A stored in the array A, and a rectangular M-by-(N-K) matrix B, stored in the array B. The block reflector H is stored in a compact WY-representation, where the elementary reflectors are in the arrays A, B and T. See the Further Details section. 🫶

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

See the Further Details section.

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 the order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. 0 <= K <= N.
T — T is DOUBLE PRECISION 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 DOUBLE PRECISION 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 the Further Details section.

LDA — LDB is INTEGER. The leading dimension of the array A. LDA >= max(1,K).

okay so yes, the manpage literally says LDB there instead of LDA — we do not fix it, we do not touch it, we report it exactly as it is. Not my typo, not your typo, it's just the lore. 🙈

B — B is DOUBLE PRECISION 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 the Further Details section.

LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M).
WORK — WORK is DOUBLE PRECISION 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**T ) * ( A_in ) =
   ( B_out )        ( B_in )                          ( B_in )
              = ( I - ( V1 ) * T * ( V1**T, V2**T ) ) * ( A_in )
                      ( V2 )                            ( B_in )

On input:

( A_in )
( B_in )

consists of two block columns:

( 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 )
( B_out )

consists of two block columns:

( 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**T ) * ( A1_in )
                  ( B1_out )        (     0 )                          (     0 )

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

okay so like, if IDENT != 'I', the computation for column block 1 is giving:

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

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

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

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

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

but if IDENT == 'I', it's a whole different vibe 💖. The operation for column block 1:

                  A1_out: = A1_in - V1*T**A1_in

                  B1_out: = - V2*T**A1_in

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

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

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

✨ (2) Description of the Algorithmic Computation ✨

so like, step one: we compute column block 2, aka A2 and B2. For that we need the K-by-(N-K) rectangular workspace matrix W2, same size as A2. W2 lives in WORK(1:K,1:(N-K)).

step two: we compute column block 1, aka A1 and B1. Now we need the K-by-K square workspace matrix W1, same size as A1. W1 is stored in 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.

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**T) * W2 = (unit_lower_tr_of_(A1)**T) * W2
                  col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W1 = (unit_lower_tr_of_(A1)**T) * 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**T) * B2 = W2 + (B1**T) * 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**T) * W2
                                 = (unit_lower_tr_of_(A1)**T) * W2
                 end if
                 col2_(3)  W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W1
                               = (unit_lower_tr_of_(A1)**T) * 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

✨ subroutine dlarfg ✨

   subroutine dlarfg (integer N, double precision ALPHA, double precision,
       dimension( * ) X, integer INCX, double precision TAU)

DLARFG generates an elementary reflector (Householder matrix). periodt 💅

Purpose:

DLARFG generates a real elementary reflector H of order n, such that

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

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

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

where tau is a real scalar and v is a real (n-1)-element vector.

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

Otherwise 1 <= tau <= 2.

Parameters

N — N is INTEGER. The order of the elementary reflector.
ALPHA — ALPHA is DOUBLE PRECISION. On entry, the value alpha. On exit, it is overwritten with the value beta.
X — X is DOUBLE PRECISION 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 DOUBLE PRECISION. The value tau.

Author

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

✨ subroutine dlarfgp ✨

   subroutine dlarfgp (integer N, double precision ALPHA, double precision,
       dimension( * ) X, integer INCX, double precision TAU)

DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. so extra ✨

Purpose:

DLARFGP generates a real elementary reflector H of order n, such that

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

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

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

where tau is a real scalar and v is a real (n-1)-element vector.

If the elements of x are all zero, 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 DOUBLE PRECISION. On entry, the value alpha. On exit, it is overwritten with the value beta.
X — X is DOUBLE PRECISION 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 DOUBLE PRECISION. The value tau.

Author

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

✨ subroutine dlarft ✨

   subroutine dlarft (character DIRECT, character STOREV, integer N, integer K,
       double precision, dimension( ldv, * ) V, integer LDV, double precision,
       dimension( * ) TAU, double precision, dimension( ldt, * ) T, integer
       LDT)

DLARFT forms the triangular factor T of a block reflector H = I - vtvH. iconic 🌈

Purpose:

DLARFT forms the triangular factor T of a real 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**T

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**T * 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (K). TAU(i) must contain the scalar factor of the elementary reflector H(i).
T — T is DOUBLE PRECISION 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. Periodt.

✨ Author ✨

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 Further Details: 🦋

Okay so 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. Iykyk.

             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 )

🌸 subroutine dlarfx 🌸

   subroutine dlarfx (character SIDE, integer M, integer N, double precision,
       dimension( * ) V, double precision TAU, double precision, dimension(
       ldc, * ) C, integer LDC, double precision, dimension( * ) WORK)

Bestie, DLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10. It's giving optimized.

✨ Purpose ✨

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

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

where tau is a real scalar and v is a real 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 DOUBLE PRECISION array, dimension (M) if SIDE = 'L' or (N) if SIDE = 'R' The vector v in the representation of H.
TAU — TAU is DOUBLE PRECISION The value tau in the representation of H.
C — C is DOUBLE PRECISION 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 >= (1,M).
WORK — WORK is DOUBLE PRECISION array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R' WORK is not referenced if H has order < 11.

💖 Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 subroutine dlarfy 🦋

   subroutine dlarfy (character UPLO, integer N, double precision, dimension( *
       ) V, integer INCV, double precision TAU, double precision, dimension(
       ldc, * ) C, integer LDC, double precision, dimension( * ) WORK)

DLARFY — that's it, that's the name. No cap.

✨ Purpose ✨

DLARFY applies an elementary reflector, or Householder matrix, H, to an n x n symmetric matrix C, from both the left and the right.

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.

🎀 Parameters 🎀

UPLO — UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix C is stored. = 'U': Upper triangle = 'L': Lower triangle
N — N is INTEGER The number of rows and columns of the matrix C. N >= 0.
V — V is DOUBLE PRECISION array, dimension (1 + (N-1)*abs(INCV)) The vector v as described above.
INCV — INCV is INTEGER The increment between successive elements of v. INCV must not be zero.
TAU — TAU is DOUBLE PRECISION The value tau as described above.
C — C is DOUBLE PRECISION array, dimension (LDC, N) On entry, the matrix C. On exit, C is overwritten by H * C * H'.
LDC — LDC is INTEGER The leading dimension of the array C. LDC >= max( 1, N ).
WORK — WORK is DOUBLE PRECISION array, dimension (N)

💖 Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 subroutine dlargv 🌈

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

DLARGV generates a vector of plane rotations with real cosines and real sines. Periodt.

✨ Purpose ✨

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

   (  c(i)  s(i) ) ( x(i) ) = ( a(i) )
   ( -s(i)  c(i) ) ( y(i) ) = (   0  )

🎀 Parameters 🎀

N — N is INTEGER The number of plane rotations to be generated.
X — X is DOUBLE PRECISION array, dimension (1+(N-1)*INCX) On entry, the vector x. On exit, x(i) is overwritten by a(i), for i = 1,...,n.
INCX — INCX is INTEGER The increment between elements of X. INCX > 0.
Y — Y is DOUBLE PRECISION array, dimension (1+(N-1)*INCY) On entry, the vector y. On exit, the sines of the plane rotations.
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.
INCC — INCC is INTEGER The increment between elements of C. INCC > 0.

💖 Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

👑 subroutine dlarrv 👑

   subroutine dlarrv (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, double precision, dimension( ldz,
       * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ, double precision,
       dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

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

✨ Purpose ✨

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

🎀 Parameters 🎀

N — N is INTEGER The order of the matrix. N >= 0.
VL — VL is 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 DOUBLE PRECISION Upper bound of the interval that contains the desired eigenvalues. VL < VU. Note: VU is currently not used by this implementation of DLARRV, VU is passed to DLARRV because it could be used compute gaps on the right end of the extremal eigenvalues. However, with not much initial accuracy in LAMBDA and VU, the formula can lead to an overestimation of the right gap and thus to inadequately early RQI 'convergence'. This is currently prevented this by forcing a small right gap. And so it turns out that VU is currently not used by this implementation of DLARRV.
D — D is DOUBLE PRECISION array, dimension (N) On entry, the N diagonal elements of the diagonal matrix D. On exit, D may be overwritten.
L — L is 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 DOUBLE PRECISION The minimum pivot allowed in the Sturm sequence.
ISPLIT — ISPLIT is 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 INTEGER The total number of input eigenvalues. 0 <= M <= N.
DOL — DOL is INTEGER
DOU — DOU is 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 DOUBLE PRECISION
RTOL1
RTOL1 is DOUBLE PRECISION
RTOL2 is DOUBLE PRECISION — Parameters for bisection. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
W — 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 — DOUBLE PRECISION array, dimension (N) — The first M elements contain the semiwidth of the uncertainty interval of the corresponding eigenvalue in W
WGAP — DOUBLE PRECISION array, dimension (N) — The separation from the right neighbor eigenvalue in W.
IBLOCK — INTEGER array, dimension (N) — The indices of the blocks (submatrices) associated with the corresponding eigenvalues in W; 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.
INDEXW — INTEGER array, dimension (N) — The indices of the eigenvalues within each block (submatrix); for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
GERS — DOUBLE PRECISION array, dimension (2*N) — The N Gerschgorin intervals (the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should be computed from the original UNshifted matrix.
Z — DOUBLE PRECISION 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: the user must ensure that at least max(1,M) columns are supplied in the array Z.
LDZ is INTEGER — The leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', LDZ >= max(1,N).
ISUPPZ — INTEGER array, dimension ( 2*max(1,M) ) — The support of the eigenvectors in Z, i.e., the indices indicating the nonzero elements in Z. The I-th eigenvector is nonzero only in elements ISUPPZ( 2*I-1 ) through ISUPPZ( 2*I ).
WORK — DOUBLE PRECISION array, dimension (12*N)
IWORK — INTEGER array, dimension (7*N)
INFO is INTEGER
- = 0: successful exit
- > 0: A problem occurred in DLARRV.
- < 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 they might be trading in precision when they decrease MINRGP.
- =-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.

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 `dlartv` (`integer N`, `double precision, dimension( * ) X`, `integer INCX`, `double precision, dimension( * ) Y`, `integer INCY`, `double precision, dimension( * ) C`, `double precision, dimension( * ) S`, `integer INCC`)

DLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair of vectors. Rotation energy, fr fr. 💫

Purpose 💫

DLARTV applies a vector of real plane rotations to elements of the real vectors x and y. For i = 1,2,...,n

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

Parameters 📋

N is INTEGER — The number of plane rotations to be applied.
X — DOUBLE PRECISION array, dimension (1+(N-1)*INCX) — The vector x.
INCX is INTEGER — The increment between elements of X. INCX > 0.
Y — DOUBLE PRECISION array, dimension (1+(N-1)*INCY) — The vector y.
INCY is INTEGER — The increment between elements of Y. INCY > 0.
C — DOUBLE PRECISION array, dimension (1+(N-1)*INCC) — The cosines of the plane rotations.
S — DOUBLE PRECISION array, dimension (1+(N-1)*INCC) — The sines of the plane rotations.
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 `dlaswp` (`integer N`, `double precision, dimension( lda, * ) A`, `integer LDA`, `integer K1`, `integer K2`, `integer, dimension( * ) IPIV`, `integer INCX`)

DLASWP performs a series of row interchanges on a general rectangular matrix. It's giving pivot energy, no cap. 🔥

Purpose 💫

DLASWP 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 is INTEGER — The number of columns of the matrix A.
A — DOUBLE PRECISION 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.
LDA is INTEGER — The leading dimension of the array A.
K1 is INTEGER — The first element of IPIV for which a row interchange will be done.
K2 is INTEGER — (K2-K1+1) is the number of elements of IPIV for which a row interchange will be done.
IPIV — 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 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 `dlat2s` (`character UPLO`, `integer N`, `double precision, dimension( lda, * ) A`, `integer LDA`, `real, dimension( ldsa, * ) SA`, `integer LDSA`, `integer INFO`)

DLAT2S converts a double-precision triangular matrix to a single-precision triangular matrix. Precision glow-down, but make it safe. 🦋

Purpose 💫

DLAT2S converts a DOUBLE PRECISION triangular matrix, SA, to a SINGLE PRECISION triangular matrix, A.

RMAX is the overflow for the SINGLE PRECISION arithmetic DLAS2S 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.

Parameters 📋

UPLO is CHARACTER*1
- = 'U': A is upper triangular;
- = 'L': A is lower triangular.
N is INTEGER — The number of rows and columns of the matrix A. N >= 0.
A — DOUBLE PRECISION array, dimension (LDA,N) — On entry, the N-by-N triangular coefficient matrix A.
LDA is INTEGER — The leading dimension of the array A. LDA >= max(1,N).
SA — REAL 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.
LDSA is INTEGER — The leading dimension of the array SA. LDSA >= max(1,M).
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.

🦋 `dlatbs` — triangular banded solve, but make it safe

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

DLATBS solves a triangular banded system of equations. no cap, that's literally the whole job description 💅

💖 Purpose

DLATBS solves one of the triangular systems

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

bestie if A is singular it literally just hands you back a non-trivial solution to A*x = 0 and sets s = 0. that's not a bug that's the whole design. 🫶

🎀 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**T* x = s*b (Conjugate transpose = 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.
KD — KD is INTEGER. The number of subdiagonals or superdiagonals in the triangular matrix A. KD >= 0.
AB — AB is DOUBLE PRECISION array, dimension (LDAB,N). 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).

LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= KD+1.
X — X is DOUBLE PRECISION 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  or  A**T* 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:

A rough bound on x is computed; if that is less than overflow, DTBSV 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 DTBSV 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. The basic algorithm for A upper triangular is

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

We simultaneously compute two bounds

     G(j) = bound on ( b(i) - A[1:i-1,i]**T * 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 DTBSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).

🌈 `dlatdf` — the reciprocal Dif-estimate contributor

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

DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. iykyk 💜

💫 Purpose

DLATDF uses the LU factorization of the n-by-n matrix Z computed by DGETC2 and computes a contribution to the reciprocal Dif-estimate by solving Z * x = b for x, and choosing the r.h.s. b such that the norm of x is as large as possible. On entry RHS = b holds the contribution from earlier solved sub-systems, and on return RHS = x.

The factorization of Z returned by DGETC2 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 DGECON, 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.
- IJOB .ne. 2: Local look ahead strategy where all entries of the r.h.s. b is chosen as either +1 or -1 (Default).
N — N is INTEGER. The number of columns of the matrix Z.
Z — Z is DOUBLE PRECISION array, dimension (LDZ, N). On entry, the LU part of the factorization of the n-by-n matrix Z computed by DGETC2: Z = P * L * U * Q
LDZ — LDZ is INTEGER. The leading dimension of the array Z. LDA >= max(1, N).
RHS — RHS is DOUBLE PRECISION 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 DTGSYL, 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 DTGSY2 is called by STGSYL.

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 DTGSY2 is called by DTGSYL.

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 IMINF-95.05, Departement of Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.

🦋 `dlatps` 🦋

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

DLATPS solves a triangular system of equations with the matrix held in packed storage. cute lil compressed matrix, she's frugal 💅

Purpose ✨

DLATPS solves one of the triangular systems

A *x = s*b  or  A**T*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, 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 DTPSV 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. no cap, that's the whole vibe 🌈

Parameters 🌸

UPLO — CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS — 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**T* x = s*b (Conjugate transpose = Transpose)
DIAG — CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN — 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 — INTEGER. The order of the matrix A. N >= 0.
AP — DOUBLE PRECISION 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 — DOUBLE PRECISION array, dimension (N). On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x. bye b, hello x 👋
SCALE — DOUBLE PRECISION. The scaling factor s for the triangular system

A * x = s*b  or  A**T* 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 — 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 — 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, DTPSV 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 DTPSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).

Similarly, a row-wise scheme is used to solve A**T*x = b. The basic algorithm for A upper triangular is

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

We simultaneously compute two bounds

G(j) = bound on ( b(i) - A[1:i-1,i]**T * 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 DTPSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow). math girlies stay winning ✨

🎀 `dlatrd` 🎀

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

DLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an orthogonal similarity transformation. she's a transformer fr fr 🦋

Purpose ✨

DLATRD reduces NB rows and columns of a real symmetric matrix A to symmetric tridiagonal form by an orthogonal similarity transformation Q**T * 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', DLATRD reduces the last NB rows and columns of a matrix, of which the upper triangle is supplied; if UPLO = 'L', DLATRD reduces the first NB rows and columns of a matrix, of which the lower triangle is supplied.

This is an auxiliary routine called by DSYTRD. just a lil helper 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.
NB — INTEGER. The number of rows and columns to be reduced.
A — DOUBLE PRECISION array, dimension (LDA,N). On entry, the symmetric 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 orthogonal 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 orthogonal matrix Q as a product of elementary reflectors.

See Further Details. 👀

LDA — INTEGER. The leading dimension of the array A. LDA >= (1,N).
E — 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 — DOUBLE PRECISION 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 — a DOUBLE PRECISION array, dimension (LDW,NB). This is the n-by-nb matrix W, and bestie it's the whole reason we're here: you need it to update the unreduced part of A. No W, no update. Periodt. 💅
LDW — INTEGER. The leading dimension of the array W. Constraint: LDW >= max(1,N). Respect it or the vibes (and your memory) get wrecked. 🫶

✨ Author ✨

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

(yeah, the whole squad is credited, iconic behavior honestly 🎀)

💖 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**T

where tau is a real scalar, and v is a real 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**T

where tau is a real scalar, and v is a real 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 symmetric rank-2k update of the form: A := A - V*W**T - W*V**T.

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

(ngl staring at those little diagrams until it clicks is a rite of passage 🌸)

🦋 subroutine dlatrs

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

DLATRS solves a triangular system of equations with the scale factor set to prevent overflow.

Purpose:

DLATRS solves one of the triangular systems

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

with scaling to prevent overflow. Here A is an upper or lower triangular matrix, 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 DTRSV 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 — CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS — 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**T* x = s*b (Conjugate transpose = Transpose)
DIAG — CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN — 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 — INTEGER. The order of the matrix A. Constraint: N >= 0.
A — DOUBLE PRECISION 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 — INTEGER. The leading dimension of the array A. Constraint: LDA >= max (1,N).
X — DOUBLE PRECISION 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 — DOUBLE PRECISION. The scaling factor s for the triangular system

A * x = s*b  or  A**T* 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 — 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 — 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, DTRSV 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 DTRSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).

Similarly, a row-wise scheme is used to solve A**T*x = b. The basic algorithm for A upper triangular is

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

We simultaneously compute two bounds

G(j) = bound on ( b(i) - A[1:i-1,i]**T * 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 DTRSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).

🌈 subroutine dlatrs3

subroutine dlatrs3 (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, integer NRHS, double precision, dimension(
    lda, * ) A, integer LDA, double precision, dimension( ldx, * ) X,
    integer LDX, double precision, dimension( * ) SCALE, double precision,
    dimension( * ) CNORM, double precision, dimension( * ) WORK, integer
    LWORK, integer INFO)

DLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow.

Purpose:

DLATRS3 solves one of the triangular systems

A * X = B * diag(scale)  or  A**T * X = B * diag(scale)

with scaling to prevent overflow. Here A is an upper or lower triangular matrix, A**T denotes the transpose of A. X and B are n by nrhs matrices and scale is an nrhs element vector of scaling factors. A scaling factor scale(j) is usually less than or equal to 1, chosen such that X(:,j) is less than the overflow threshold. If the matrix A is singular (A(j,j) = 0 for some j), then a non-trivial solution to A*X = 0 is returned. If the system is so badly scaled that the solution cannot be represented as (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned.

This is a BLAS-3 version of LATRS for solving several right hand sides simultaneously. 🔥

Parameters

UPLO — CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS —
TRANS — this one's a CHARACTER*1 and it picks which operation gets applied to A, bestie. The options are:
- = 'N': Solve A * x = s*b (No transpose)
- = 'T': Solve A**T* x = s*b (Transpose)
- = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose)
DIAG — CHARACTER*1, and it's the vibe check on whether A is a unit triangular matrix or not:
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN — CHARACTER*1, and it's literally just "has CNORM 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 — INTEGER, the order of the matrix A. N >= 0. Not negotiable. 🙅
NRHS — INTEGER, the number of columns of X. NRHS >= 0.
A — DOUBLE PRECISION array, dimension (LDA,N). It's the triangular matrix A, and here's the tea: 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. And if DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1. They said "ignore that part" and meant it. 💅
LDA — INTEGER, the leading dimension of the array A. LDA >= max (1,N).
X — DOUBLE PRECISION array, dimension (LDX,NRHS). On entry, it's the right hand side B of the triangular system. On exit, X is overwritten by the solution matrix X. glow up moment ✨
LDX — INTEGER, the leading dimension of the array X. LDX >= max (1,N).
SCALE — DOUBLE PRECISION array, dimension (NRHS). The scaling factor s(k) is for the triangular system A * x(:,k) = s(k)*b(:,k) or A**T* x(:,k) = s(k)*b(:,k). If SCALE = 0, the matrix A is singular or badly scaled. If A(j,j) = 0 is encountered, a non-trivial vector x(:,k) that is an exact or approximate solution to A*x(:,k) = 0 is returned. If the system so badly scaled that solution cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0 is returned. Yeah. It's giving "we tried our best". 😭
CNORM — 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.

WORK — DOUBLE PRECISION array, dimension (LWORK). On exit, if INFO = 0, WORK(1) returns the optimal size of WORK.
LWORK — LWORK is INTEGER LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.

If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal dimensions of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. It's the "just tell me how much room I need, I'm not moving in yet" mode. 🫶

Parameters 👑

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.

Further Details:

That's the whole section, bestie. Nothing under this heading in the source — no extra lore. Moving on.

✨ subroutine `dlauu2` ✨

subroutine dlauu2 (character UPLO, integer N, double precision, dimension( lda, * ) A, integer LDA, integer INFO)

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

Purpose:

DLAUU2 computes the product U * U**T or L**T * 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. It's the small-batch, hand-crafted artisanal version. 🧚

Parameters

UPLO — CHARACTER*1. Specifies whether the triangular factor stored in the array A is upper or lower triangular:
- = 'U': Upper triangular
- = 'L': Lower triangular
N — INTEGER. The order of the triangular factor U or L. N >= 0.
A — DOUBLE PRECISION 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**T; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**T * L.
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N).
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 `dlauum` ✨

subroutine dlauum (character UPLO, integer N, double precision, dimension( lda, * ) A, integer LDA, integer INFO)

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

Purpose:

DLAUUM computes the product U * U**T or L**T * 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 blocked form of the algorithm, calling Level 3 BLAS. Same energy, but it eats bigger bites at a time. 🔥

Parameters

UPLO — CHARACTER*1. Specifies whether the triangular factor stored in the array A is upper or lower triangular:
- = 'U': Upper triangular
- = 'L': Lower triangular
N — INTEGER. The order of the triangular factor U or L. N >= 0.
A — DOUBLE PRECISION 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**T; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**T * L.
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N).
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 `drscl` 🦋

subroutine drscl (integer N, double precision SA, double precision, dimension( * ) SX, integer INCX)

DRSCL multiplies a vector by the reciprocal of a real scalar. Math but make it fashion.

Purpose:

DRSCL multiplies an n-element real vector x by the real scalar 1/a.

This is done without overflow or underflow as long as the final result x/a does not overflow or underflow. Iconic behaviour honestly. 💅

Parameters

N — INTEGER. The number of components of the vector x.
SA — 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. 👀 read that again, bestie.
SX — DOUBLE PRECISION array, dimension (1+(N-1)*abs(INCX)). The n-element vector x.
INCX — 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 `dtprfb` 🌈

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

DTPRFB applies a real 'triangular-pentagonal' block reflector to a real matrix, which is composed of two blocks. Yeah, two blocks. It's a whole thing and we're here for it. 🎀

Purpose:

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

Parameters

SIDE — CHARACTER*1
- = 'L': apply H or H**T from the Left
- = 'R': apply H or H**T from the Right
TRANS — CHARACTER*1
- = 'N': apply H (No transpose)
- = 'T': apply H**T (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': Columns
- = 'R': Rows
M — INTEGER. The number of rows of the matrix B. M >= 0.
N — N is INTEGER. it's the number of columns of the matrix B, and like… N >= 0, bestie. no negatives here 💅
K — K is INTEGER. this is the order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. K >= 0. yes queen, zero is allowed.
L — L is INTEGER. the order of the trapezoidal part of V. and the rule is K >= L >= 0. go read Further Details if you're confused, i got you 🫶
V — V is DOUBLE PRECISION array, dimension
- (LDV,K) if STOREV = 'C'
- (LDV,M) if STOREV = 'R' and SIDE = 'L'
- (LDV,N) if STOREV = 'R' and SIDE = 'R'

this is the pentagonal matrix V, and it's carrying the elementary reflectors H(1), H(2), ..., H(K). See Further Details. it's giving main character energy honestly ✨

LDV — LDV is INTEGER. the leading dimension of the array V. so:
- if STOREV = 'C' and SIDE = 'L', then LDV >= max(1,M);
- if STOREV = 'C' and SIDE = 'R', then LDV >= max(1,N);
- if STOREV = 'R', then LDV >= K.
T — T is DOUBLE PRECISION array, dimension (LDT,K). it's the triangular K-by-K matrix T in the representation of the block reflector. she's triangular and she knows it 👑
LDT — LDT is INTEGER. the leading dimension of the array T. and LDT >= K, periodt.
A — A is DOUBLE PRECISION array, dimension
- (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'

on entry, it's the K-by-N or M-by-K matrix A. on exit, A is overwritten by the corresponding block of H*C or H**T*C or C*H or C*H**T. See Further Details. yes bestie, overwritten, she's gone, don't come crying 😭

LDA — LDA is INTEGER. the leading dimension of the array A.
- If SIDE = 'L', LDA >= max(1,K);
- If SIDE = 'R', LDA >= max(1,M).
B — B is DOUBLE PRECISION 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**T*C or C*H or C*H**T. See Further Details. same fate as A, sorry babe 🥀
LDB — LDB is INTEGER. the leading dimension of the array B. and LDB >= max(1,M).
WORK — WORK is DOUBLE PRECISION array, dimension
- (LDWORK,N) if SIDE = 'L',
- (LDWORK,K) if SIDE = 'R'.
LDWORK — LDWORK is 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', C = [A]
                                    [B].

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

If SIDE = 'L' and DIRECT = 'B', 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.

🧚 subroutine slatrd

subroutine slatrd (character UPLO, integer N, integer NB, real, dimension(
    lda, * ) A, integer LDA, real, dimension( * ) E, real, dimension( * )
    TAU, real, dimension( ldw, * ) W, integer LDW)

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

✨ Purpose

SLATRD reduces NB rows and columns of a real symmetric matrix A to symmetric tridiagonal form by an orthogonal similarity transformation Q**T * 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', SLATRD reduces the last NB rows and columns of a matrix, of which the upper triangle is supplied; if UPLO = 'L', SLATRD reduces the first NB rows and columns of a matrix, of which the lower triangle is supplied.

This is an auxiliary routine called by SSYTRD. yeah she's the support act but she's lowkey the moment 💅

🌸 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.
NB — NB is INTEGER. the number of rows and columns to be reduced.
A — A is REAL array, dimension (LDA,N). on entry, the symmetric 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 orthogonal 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 orthogonal matrix Q as a product of elementary reflectors.

See Further Details.

LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= (1,N).
E — E is REAL 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 REAL 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 REAL 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**T

where tau is a real scalar, and v is a real 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**T

where tau is a real scalar, and v is a real 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 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 slatrs3

subroutine slatrs3 (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, integer NRHS, real, dimension( lda, * ) A,
    integer LDA, real, dimension( ldx, * ) X, integer LDX, real, dimension(
    * ) SCALE, real, dimension( * ) CNORM, real, dimension( * ) WORK,
    integer LWORK, integer INFO)

SLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow.

✨ Purpose

SLATRS3 solves one of the triangular systems

   A * X = B * diag(scale)  or  A**T * X = B * diag(scale)

...with scaling to prevent overflow, bestie. So like, here A is an upper or lower triangular matrix, A**T denotes the transpose of A. X and B are n by nrhs matrices and scale is an nrhs-element vector of scaling factors. A scaling factor scale(j) is usually less than or equal to 1, chosen such that X(:,j) is less than the overflow threshold. If the matrix A is singular (A(j,j) = 0 for some j), then a non-trivial solution to A*X = 0 is returned. If the system is so badly scaled that the solution cannot be represented as (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned. 💅

And like, this is a BLAS-3 version of LATRS for solving several right hand sides simultaneously. Periodt. ✨

🌸 Parameters 🌸

UPLO — UPLO is CHARACTER*1. Tells you 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**T* x = s*b (Conjugate transpose = 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.
NRHS — NRHS is INTEGER. The number of columns of X. NRHS >= 0.
A — A is REAL 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 REAL array, dimension (LDX,NRHS). On entry, the right hand side B of the triangular system. On exit, X is overwritten by the solution matrix X.
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= max (1,N).
SCALE — SCALE is REAL array, dimension (NRHS). The scaling factor s(k) is for the triangular system A * x(:,k) = s(k)*b(:,k) or A**T* x(:,k) = s(k)*b(:,k). If SCALE = 0, the matrix A is singular or badly scaled. If A(j,j) = 0 is encountered, a non-trivial vector x(:,k) that is an exact or approximate solution to A*x(:,k) = 0 is returned. If the system so badly scaled that solution cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0 is returned.
CNORM — CNORM is REAL 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.
WORK — WORK is REAL array, dimension (LWORK). On exit, if INFO = 0, WORK(1) returns the optimal size of WORK.
LWORK — LWORK is INTEGER. LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.

If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal dimensions of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. 👑

Parameters

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:

🦋 subroutine zlatrs3 🦋

subroutine zlatrs3 (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, integer NRHS, complex*16, dimension( lda, * )
    A, integer LDA, complex*16, dimension( ldx, * ) X, integer LDX, double
    precision, dimension( * ) SCALE, double precision, dimension( * ) CNORM,
    double precision, dimension( * ) WORK, integer LWORK, integer INFO)

ZLATRS3 solves a triangular system of equations with the scale factors set to prevent overflow. 💅✨

🌈 Purpose: 🌈

ZLATRS3 solves one of the triangular systems

A * X = B * diag(scale),  A**T * X = B * diag(scale), or
A**H * X = B * diag(scale)

with scaling to prevent overflow, bestie. So like, 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-by-nrhs matrices and scale is an nrhs-element vector of scaling factors. A scaling factor scale(j) is usually less than or equal to 1, chosen such that X(:,j) is less than the overflow threshold. If the matrix A is singular (A(j,j) = 0 for some j), then a non-trivial solution to A*X = 0 is returned. If the system is so badly scaled that the solution cannot be represented as (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned. 😭

And like, this is a BLAS-3 version of LATRS for solving several right hand sides simultaneously. Iconic. ✨

🎀 Parameters 🎀

UPLO — UPLO is CHARACTER*1. Tells you 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**T* 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.
NRHS — NRHS is INTEGER. The number of columns of X. NRHS >= 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 (LDX,NRHS). On entry, the right hand side B of the triangular system. On exit, X is overwritten by the solution matrix X.
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= max (1,N).
SCALE — SCALE is DOUBLE PRECISION array, dimension (NRHS). The scaling factor s(k) is for the triangular system A * x(:,k) = s(k)*b(:,k) or A**T* x(:,k) = s(k)*b(:,k). If SCALE = 0, the matrix A is singular or badly scaled. If A(j,j) = 0 is encountered, a non-trivial vector x(:,k) that is an exact or approximate solution to A*x(:,k) = 0 is returned. If the system so badly scaled that solution cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0 is returned.
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.
WORK — WORK is DOUBLE PRECISION array, dimension (LWORK). On exit, if INFO = 0, WORK(1) returns the optimal size of WORK.
LWORK — LWORK is INTEGER. LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.

If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal dimensions of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. 👑

Parameters

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

…and last but never least, the closing credits give a lil shoutout to NAG Ltd. 💖 yup, they're in the building, they're part of the story, iconic behavior.

🌸 Further Details 🌸

okay so this heading is literally here ✨ and the actual details are the ones we already spilled above, bestie — nothing new gets dropped under this label, it's just the vibe check at the end of the function docs. no cap, that's the whole section. 🫶

🎀 Author 🎀

Generated automatically by Doxygen for LAPACK from the source code. 🤖✨

so like… nobody hand-wrote this bad boy, Doxygen just read the LAPACK source and went "i got you, bestie" and generated the whole page. machine-made, still ate. 💅

lowkey a whole mood that this entire man page is auto-generated and still this detailed. main character energy from a doc tool, fr fr. 😭👑

and that's a wrap on doubleOTHERauxiliary(3) babes!! 🦋🌸 this page came to you from LAPACK, version 3.11.0, dated Sun Nov 27 2022 — screenshot it, cite it, love it, and go read some source code. mwah. 💖✨🎀🫶

📦 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

doubleOTHERauxiliary(3)

🌸 doubleOTHERauxiliary(3) — LAPACK auxiliary, bestie edition 🌸

✨ NAME ✨

✨ SYNOPSIS ✨

🌸 Functions 🌸

Detailed Description 💖

Function Documentation ✨

subroutine clatrs3

Purpose:

Parameters

Parameters

Author

Further Details:

subroutine dlabrd

Purpose:

Parameters

👑 Author

📚 Further Details:

🌈 subroutine dlacn2 (integer N, double precision, dimension( * ) V, double precision, dimension( * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE, integer, dimension( 3 ) ISAVE)

🎯 Purpose:

🧾 Parameters

👑 Author

📚 Further Details:

🫶 Contributors:

📖 References:

🌈 subroutine dlacon (integer N, double precision, dimension( * ) V, double precision, dimension( * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE)

🎯 Purpose:

🧾 Parameters

👑 Author

🫶 Contributors:

📖 References:

🌈 subroutine dladiv (double precision A, double precision B, double precision C, double precision D, double precision P, double precision Q)

🎯 Purpose:

🧾 Parameters

👑 Author

🎯 Purpose:

🧾 Parameters

Author 💖

subroutine dlaexc (logical WANTQ, integer N, double precision, dimension( ldt, * ) T, integer LDT, double precision, dimension( ldq, * ) Q, integer LDQ, integer J1, integer N1, integer N2, double precision, dimension( * ) WORK, integer INFO) 🔥

Purpose ✨

Parameters 🎀

Author 💖

subroutine dlag2 (double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( ldb, * ) B, integer LDB, double precision SAFMIN, double precision SCALE1, double precision SCALE2, double precision WR1, double precision WR2, double precision WI) 🌈

Purpose ⭐

Parameters 🧚

Author 💖

subroutine dlag2s (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, real, dimension( ldsa, * ) SA, integer LDSA, integer INFO) 🌟

Purpose ⭐

Parameters 🧚

Author 💖

🌸 subroutine dlags2 🌸

✨ Purpose ✨

🎀 Parameters 🎀

👑 Author 👑

🌈 subroutine dlagtm 🌈

✨ Purpose ✨

🎀 Parameters 🎀

👑 Author 👑

🦋 subroutine dlagv2 🦋

✨ Purpose ✨

🎀 Parameters 🎀

👑 Author 👑

💖 Contributors 💖

🌟 subroutine dlahqr 🌟

✨ Purpose ✨

🎀 Parameters 🎀

Author

Further Details:

🦋 subroutine dlahr2

Purpose:

Parameters

Author

Further Details:

References:

💜 subroutine dlaic1

Purpose:

Parameters

💫 Author 💫

🦋 dlaln2 🦋

🌸 Purpose 🌸

`subroutine clatrs3`

`subroutine dlabrd`

🌈 `subroutine dlacn2 (integer N, double precision, dimension( * ) V, double precision, dimension( * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE, integer, dimension( 3 ) ISAVE)`

🌈 `subroutine dlacon (integer N, double precision, dimension( * ) V, double precision, dimension( * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE)`

🌈 `subroutine dladiv (double precision A, double precision B, double precision C, double precision D, double precision P, double precision Q)`

`subroutine dlaexc (logical WANTQ, integer N, double precision, dimension( ldt, * ) T, integer LDT, double precision, dimension( ldq, * ) Q, integer LDQ, integer J1, integer N1, integer N2, double precision, dimension( * ) WORK, integer INFO)` 🔥

`subroutine dlag2 (double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( ldb, * ) B, integer LDB, double precision SAFMIN, double precision SCALE1, double precision SCALE2, double precision WR1, double precision WR2, double precision WI)` 🌈

`subroutine dlag2s (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, real, dimension( ldsa, * ) SA, integer LDSA, integer INFO)` 🌟

🦋 `dlaln2` 🦋

🦋 `dlangt` 🦋

🦋 `dlanhs` 🦋

🦋 `dlansb` 🦋

`double precision function dlansp (character NORM, character UPLO, integer N, double precision, dimension( * ) AP, double precision, dimension( * ) WORK)`

`double precision function dlantb (character NORM, character UPLO, character DIAG, integer N, integer K, double precision, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) WORK)`

`double precision function dlantp (character NORM, character UPLO, character DIAG, integer N, double precision, dimension( * ) AP, double precision, dimension( * ) WORK)`

🦋 `dlantr` 💜

🔥 `dlanv2` 🔥

🌈 `dlapll` 🌈

💅 `dlapmr` 💅

🦋 `dlapmt` 🦋

🌸 `dlaqp2` 🌸

🧮 `dlaqp2` — the little pivoting QR bestie

🌈 subroutine `dlaqps` 🌈

💜 subroutine `dlaqr0` 💜