complexOTHERauxiliary(3)

📦

version 3.11.0 · 2 versions

3.11.0 · debian bookworm (latest)
3.10.0 · ubuntu jammy

, in debian bookworm

complexOTHERauxiliary(3) — LAPACK manual, bestie 💖

okay so this is the complexOTHERauxiliary(3) page from the LAPACK manual ✨. It’s part 1 of 25, so settle in 🍓.

✨ NAME ✨

complexOTHERauxiliary — complex

part 1 of 25, don’t go anywhere 💅

✨ SYNOPSIS ✨

✨ Functions ✨

okay so here's the entire function lineup from this section, bestie. every signature is copied exactly, every description is translated into our language, and absolutely nothing is skipped. 💖

subroutine clabrd (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY) — CLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form. it's giving bidiagonal glow-up, no cap. ✨
subroutine clacgv (N, X, INCX) — CLACGV conjugates a complex vector. literally just conjugation, periodt. 💅
subroutine clacn2 (N, V, X, EST, KASE, ISAVE) — CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. iykyk. 🔥
subroutine clacon (N, V, X, EST, KASE) — CLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. same energy as CLACN2, different signature. 💜
subroutine clacp2 (UPLO, M, N, A, LDA, B, LDB) — CLACP2 copies all or part of a real two-dimensional array to a complex array. real-to-complex pipeline, slay. 🦋
subroutine clacpy (UPLO, M, N, A, LDA, B, LDB) — CLACPY copies all or part of one two-dimensional array to another. it's the copy-paste icon. 🎀
subroutine clacrm (M, N, A, LDA, B, LDB, C, LDC, RWORK) — CLACRM multiplies a complex matrix by a square real matrix. complex meets real, and they ate. 🔥
subroutine clacrt (N, CX, INCX, CY, INCY, C, S) — CLACRT performs a linear transformation of a pair of complex vectors. transformation? main character energy. 👑
complex function cladiv (X, Y) — CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. she's careful, she's precise, she's not overflowing. 💅
subroutine claein (RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO) — CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. eigenvector hunt, fr fr. 🌸
subroutine claev2 (A, B, C, RT1, RT2, CS1, SN1) — CLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. tiny matrix, big output. ✨
subroutine clags2 (UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ) — CLAGS2 — okay the original description is literally just CLAGS2. that's it. no extra lore, bestie. 🙈
subroutine clagtm (TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB) — CLAGTM 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. math but make it tridiagonal. 🔥
subroutine clahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO) — CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. QR algorithm, but make it double-shift/single-shift. 💜
subroutine clahr2 (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY) — CLAHR2 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. reduction with receipts. 🧚
subroutine claic1 (JOB, J, X, SEST, W, GAMMA, SESTPR, S, C) — CLAIC1 applies one step of incremental condition estimation. one step, but make it incremental. 💫
real function clangt (NORM, N, DL, D, DU) — CLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. norm options on deck. 📏
real function clanhb (NORM, UPLO, N, K, AB, LDAB, WORK) — CLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix. Hermitian band matrix, normed. 🎶
real function clanhp (NORM, UPLO, N, AP, WORK) — CLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form. packed form, unpacked norms. 🎁
real function clanhs (NORM, N, A, LDA, WORK) — CLANHS 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 norms, no cap. 🔥
real function clanht (NORM, N, D, E) — CLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix. tridiagonal Hermitian, normed. 💜
real function clansb (NORM, UPLO, N, K, AB, LDAB, WORK) — CLANSB 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, norm vibes. 🎀
real function clansp (NORM, UPLO, N, AP, WORK) — CLANSP 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 symmetric, normed. 📦
real function clantb (NORM, UPLO, DIAG, N, K, AB, LDAB, WORK) — CLANTB 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, normed. 🔺
real function clantp (NORM, UPLO, DIAG, N, AP, WORK) — CLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. triangular packed, normed. 🎁
real function clantr (NORM, UPLO, DIAG, M, N, A, LDA, WORK) — CLANTR 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, normed. 📐
subroutine clapll (N, X, INCX, Y, INCY, SSMIN) — CLAPLL measures the linear dependence of two vectors. are they dependent? let's measure. 👀
subroutine clapmr (FORWRD, M, N, X, LDX, K) — CLAPMR rearranges rows of a matrix as specified by a permutation vector. row shuffle, as requested. 🔀
subroutine clapmt (FORWRD, M, N, X, LDX, K) — CLAPMT performs a forward or backward permutation of the columns of a matrix. column shuffle, forward or backward. 🔄
subroutine claqhb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED) — CLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ. scaling factors from cpbequ, periodt. ⚖️
subroutine claqhp (UPLO, N, AP, S, SCOND, AMAX, EQUED) — CLAQHP scales a Hermitian matrix stored in packed form. packed Hermitian, scaled. 🎀
subroutine claqp2 (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK) — CLAQP2 computes a QR factorization with column pivoting of the matrix block. QR with column pivoting, block edition. 🧱
subroutine claqps (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF) — CLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. one step, BLAS level 3 energy. ⚡
subroutine claqr0 (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO) — CLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. eigenvalues, optional Schur matrices. 👑
subroutine claqr1 (N, H, LDH, S1, S2, V) — CLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. small matrix, big shift energy. 🔥
subroutine claqr2 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK) — CLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). aggressive early deflation, she's not waiting. 💅
subroutine claqr3 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK) — CLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). same as CLAQR2, still aggressive, still early. ⚡
subroutine claqr4 (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO) — CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. eigenvalues and optional Schur matrices, again. 💜
subroutine claqr5 (WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH) — CLAQR5 performs a single small-bulge multi-shift QR sweep. small bulge, multi-shift, single sweep. 🔍
subroutine claqsb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED) — CLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. scaling factors from spbequ. ⚖️
subroutine claqsp (UPLO, N, AP, S, SCOND, AMAX, EQUED) — CLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. packed storage, sppequ scaling factors. 🎁
subroutine clar1v (N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK) — CLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. inverse column, scaled, r-th, tridiagonal. 🌈
subroutine clar2v (N, X, Y, Z, INCX, C, S, INCC) — CLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. plane rotations, both sides, 2-by-2. 🦋
subroutine clarcm (M, N, A, LDA, B, LDB, C, LDC, RWORK) — CLARCM copies all or part of a real two-dimensional array to a complex array. real to complex copy, again but different name. 🎀
subroutine clarf (SIDE, M, N, V, INCV, TAU, C, LDC, WORK) — CLARF applies an elementary reflector to a general rectangular matrix. elementary reflector, rectangular matrix, done. ✨
subroutine clarfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK) — CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix. block reflector or conjugate-transpose, choose your fighter. 💅
subroutine clarfb_gett (IDENT, M, N, K, T, LDT, A, LDA, B, LDB, WORK, LDWORK) — CLARFB_GETT — the original gives us nothing but the name here. no description, no crumbs. 💀
subroutine clarfg (N, ALPHA, X, INCX, TAU) — CLARFG generates an elementary reflector (Householder matrix). Householder matrix, generated. 🏠
subroutine clarfgp (N, ALPHA, X, INCX, TAU) — CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. same Householder energy, but beta is non-negative. 💖
subroutine clarft (DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) — CLARFT forms the triangular factor T of a block reflector H = I - vtvH. triangular factor T, block reflector, no cap. 🔺
subroutine clarfx (SIDE, M, N, V, TAU, C, LDC, WORK) — CLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10. loop unrolling when order ≤ 10, she's optimized. ⚡
subroutine clarfy (UPLO, N, V, INCV, TAU, C, LDC, WORK) — CLARFY — again, the original description is just the name. that's all, bestie. 🙈
subroutine clargv (N, X, INCX, Y, INCY, C, INCC) — CLARGV generates a vector of plane rotations with real cosines and complex sines. plane rotations, real cosines, complex sines. 🎡
subroutine clarnv (IDIST, ISEED, N, X) — CLARNV returns a vector of random numbers from a uniform or normal distribution. random numbers, uniform or normal. 🎲
subroutine clarrv (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) — CLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. eigenvectors from tridiagonal, given L, D, and eigenvalues. 🌸
subroutine clartv (N, X, INCX, Y, INCY, C, S, INCC) — CLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a pair of vectors. plane rotations to a pair of vectors. 🔄
subroutine clascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO) — CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. scalar is cto/cfrom, matrix gets multiplied. ✖️
subroutine claset (UPLO, M, N, ALPHA, BETA, A, LDA) — CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values. off-diagonal and diagonal get set. 🎯
subroutine clasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA) — CLASR applies a sequence of plane rotations to a general rectangular matrix. sequence of plane rotations, rectangular matrix. 🎡
subroutine claswp (N, A, LDA, K1, K2, IPIV, INCX) — CLASWP does a whole bunch of row swaps on a general rectangular matrix. like, a series of them. bestie it's swapping rows. that's the whole gig. 💫
subroutine clatbs (UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO) — CLATBS solves a triangular banded system of equations. triangular, banded, solved. slay. 🎀
subroutine clatdf (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV) — CLATDF takes the LU factorization of the n-by-n matrix that sgetc2 computed and works out a contribution to the reciprocal Dif-estimate. yes it's giving niche but we respect it. ✨
subroutine clatps (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO) — CLATPS solves a triangular system of equations, but the matrix is held in packed storage. compact queen. 👑
subroutine clatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW) — CLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A down to real tridiagonal form using a unitary similarity transformation. glow down fr. 🦋
subroutine clatrs (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO) — CLATRS solves a triangular system of equations with the scale factor set so it literally cannot overflow. safety first, always. 🔥
subroutine clauu2 (UPLO, N, A, LDA, INFO) — CLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices. this one's the unblocked algorithm, fyi. 🧚
subroutine clauum (UPLO, N, A, LDA, INFO) — CLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices. and this one is the blocked algorithm. same job, different vibe. 💅
subroutine crot (N, CX, INCX, CY, INCY, C, S) — CROT applies a plane rotation to a pair of complex vectors, with a real cosine and a complex sine. okay mixologist. 🍓
subroutine cspmv (UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY) — CSPMV computes a matrix-vector product for complex vectors using a complex symmetric packed matrix. 🌸
subroutine cspr (UPLO, N, ALPHA, X, INCX, AP) — CSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix. ✨
subroutine csrscl (N, SA, SX, INCX) — CSRSCL multiplies a vector by the reciprocal of a real scalar. divide-and-conquer but make it multiply. 💀
subroutine ctprfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK) — CTPRFB applies a complex 'triangular-pentagonal' block reflector to a complex matrix, and that matrix is composed of two blocks. pentagonal?? iconic honestly. 🫶
integer function icmax1 (N, CX, INCX) — ICMAX1 finds the index of the first vector element of maximum absolute value. first, not last. remember that. ⭐
integer function ilaclc (M, N, A, LDA) — ILACLC scans a matrix for its last non-zero column. where's the last one bestie. 🔍
integer function ilaclr (M, N, A, LDA) — ILACLR scans a matrix for its last non-zero row. same energy, other axis. 🌈
integer function izmax1 (N, ZX, INCX) — IZMAX1 finds the index of the first vector element of maximum absolute value. 🎀
real function scsum1 (N, CX, INCX) — SCSUM1 forms the 1-norm of the complex vector using the true absolute value. no shortcuts, only truth. 💖

✨ Detailed Description ✨

This is the group of complex other auxiliary routines. the complex-other-auxiliary-outines girlies. we contain multitudes. 💜

🎀 Function Documentation 🎀

`subroutine clabrd` `(integer M, integer N, integer NB, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( * ) TAUQ, complex, dimension( * ) TAUP, complex, dimension( ldx, * ) X, integer LDX, complex, dimension( ldy, * ) Y, integer LDY)`

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

Purpose

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

If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower bidiagonal form.

This is an auxiliary routine called by CGEBRD. (so like, it's the support act, but it ate. 🍓)

Parameters

M

M is INTEGER The number of rows in the matrix A.

N

N is INTEGER The number of columns in the matrix A.

NB

NB is INTEGER The number of leading rows and columns of A to be reduced.

A

A is COMPLEX array, dimension (LDA,N) On entry, the m by n general matrix to be reduced. On exit, the first NB rows and columns of the matrix are overwritten; the rest of the array is unchanged. If m >= n, elements on and below the diagonal in the first NB columns, with the array TAUQ, represent the unitary matrix Q as a product of elementary reflectors; and elements above the diagonal in the first NB rows, with the array TAUP, represent the unitary matrix P as a product of elementary reflectors. If m < n, elements below the diagonal in the first NB columns, with the array TAUQ, represent the unitary matrix Q as a product of elementary reflectors, and elements on and above the diagonal in the first NB rows, with the array TAUP, represent the unitary matrix P as a product of elementary reflectors. See Further Details.

LDA

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

D

D is REAL array, dimension (NB) The diagonal elements of the first NB rows and columns of the reduced matrix. D(i) = A(i,i).

E

E is REAL array, dimension (NB) The off-diagonal elements of the first NB rows and columns of the reduced matrix.

TAUQ

TAUQ is COMPLEX array, dimension (NB) The scalar factors of the elementary reflectors which represent the unitary matrix Q. See Further Details.

TAUP

TAUP is COMPLEX array, dimension (NB) The scalar factors of the elementary reflectors which represent the unitary matrix P. See Further Details.

X

X is COMPLEX array, dimension (LDX,NB) The m-by-nb matrix X required to update the unreduced part of A.

LDX

LDX is INTEGER The leading dimension of the array X. LDX >= max(1,M).

Y

Y is COMPLEX array, dimension (LDY,NB) The n-by-nb matrix Y required to update the unreduced part of A.

LDY

LDY is INTEGER The leading dimension of the array Y. LDY >= max(1,N).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details

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)

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

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

where tauq and taup are complex scalars, and v and u are complex vectors.

If m >= n, 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, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

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

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

`subroutine clacgv` `(integer N, complex, dimension( * ) X, integer INCX)`

CLACGV conjugates a complex vector. that's it, that's the function, and it devoured. 💖

Purpose

CLACGV conjugates a complex vector of length N.

Parameters

N

N is INTEGER The length of the vector X. N >= 0.

X

X is COMPLEX array, dimension (1+(N-1)*abs(INCX)) On entry, the vector of length N to be conjugated. On exit, X is overwritten with conjg(X).

INCX

INCX is INTEGER The spacing between successive elements of X.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

`subroutine clacn2` `(integer N, complex, dimension( * ) V, complex, dimension( * ) X, real EST, integer KASE, integer, dimension( 3 ) ISAVE)`

CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. reverse communication?? the way that it hands the work back to you. 💀✨

Purpose

CLACN2 estimates the 1-norm of a square, complex 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 COMPLEX array, dimension (N) On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned). rude of it honestly. 🙈

X

X is COMPLEX array, dimension (N) On an intermediate return, X should be overwritten by

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

where A**H is the conjugate transpose of A, and CLACN2 must be re-called with all the other parameters unchanged.

EST

EST is REAL On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be unchanged from the previous call to CLACN2. On exit, EST is an estimate (a lower bound) for norm(A).

KASE
KASE — KASE is INTEGER. First call? Set it to 0, bestie. On an intermediate return it'll come back as 1 or 2, and that's literally the universe telling you whether X should be overwritten by A * X or A**H * X. On the final return from CLACN2, KASE goes right back to 0. Full circle. 🦋
ISAVE — ISAVE is INTEGER array, dimension (3). It's the little memory box that saves variables between calls to SLACN2. 🎀

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details: 👑

Originally named CONEST, dated March 16, 1988.

Last modified: April, 1999

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

   CLACON     CLACN2
    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 `clacon` (integer N, complex, dimension( n ) V, complex, dimension( n ) X, real EST, integer KASE)

CLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. So iconic of it. ✨

Purpose:

CLACON estimates the 1-norm of a square, complex 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 COMPLEX array, dimension (N). On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned). Rude of W honestly but okay. 💅
X — X is COMPLEX array, dimension (N). On an intermediate return, X should be overwritten by

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

where A**H is the conjugate transpose of A, and CLACON must be re-called with all the other parameters unchanged.

EST — EST is REAL. On entry with KASE = 1 or 2 and JUMP = 3, EST should be unchanged from the previous call to CLACON. On exit, EST is an estimate (a lower bound) for norm(A).
KASE — KASE is INTEGER. On the initial call to CLACON, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether X should be overwritten by A * X or A**H * X. On the final return from CLACON, KASE will again be 0.

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

Originally named CONEST, dated March 16, 1988.

Last modified: April, 1999

Contributors: 🫶

Nick Higham, University of Manchester

References: 📚

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

subroutine `clacp2` (character UPLO, integer M, integer N, real, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)

CLACP2 copies all or part of a real two-dimensional array to a complex array. Copy-paste but make it ✨complex✨.

Purpose:

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

Parameters 🌈

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

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

M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
A — A is REAL array, dimension (LDA,N). The m by n matrix A. If UPLO = 'U', only the upper trapezium is accessed; if UPLO = 'L', only the lower trapezium is accessed. Trapezium only, no cap. 🦋
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is COMPLEX array, dimension (LDB,N). On exit, B = A in the locations specified by UPLO.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M).

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `clacpy` (character UPLO, integer M, integer N, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)

CLACPY copies all or part of one two-dimensional array to another. The blueprint behavior we all deserve. 💅

Purpose:

CLACPY copies all or part of a two-dimensional matrix A to another matrix B.

Parameters 🎀

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

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

M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
A — A is COMPLEX array, dimension (LDA,N). The m by n matrix A. If UPLO = 'U', only the upper trapezium is accessed; if UPLO = 'L', only the lower trapezium is accessed.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is COMPLEX array, dimension (LDB,N). On exit, B = A in the locations specified by UPLO.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M).

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `clacrm` (integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, complex, dimension( ldc, * ) C, integer LDC, real, dimension( * ) RWORK)

CLACRM multiplies a complex matrix by a square real matrix. Mixing worlds, main character energy only. 👑

Purpose:

CLACRM performs a very simple matrix-matrix multiplication:

        C := A * B,

where A is M by N and complex; B is N by N and real; C is M by N and complex.

Parameters ⭐

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

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `clacrt` (integer N, complex, dimension( * ) CX, integer INCX, complex, dimension( * ) CY, integer INCY, complex C, complex S)

CLACRT performs a linear transformation of a pair of complex vectors. A glow-up for two vectors at once, fr fr. ✨

Purpose:

CLACRT performs the operation

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

where c and s are complex and the vectors x and y are complex.

Parameters 🌸

N — N is INTEGER. The number of elements in the vectors CX and CY.
CX — CX is COMPLEX array, dimension (N). On input, the vector x. On output, CX is overwritten with c*x + s*y.
INCX — INCX is INTEGER. The increment between successive values of CX. INCX <> 0.
CY — CY is COMPLEX array, dimension (N). On input, the vector y. On output, CY is overwritten with -s*x + c*y.
INCY — INCY is INTEGER. The increment between successive values of CY. INCY <> 0.
C — C is COMPLEX.
S — S is COMPLEX. C and S define the matrix

      [  C   S  ].
      [ -S   C  ]

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

complex function `cladiv` (complex X, complex Y)

CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. She's careful, she's precise, she's not letting anything spill over. 🫶

Purpose:

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

Parameters 🎀

X — X is COMPLEX.
Y — Y is COMPLEX. The complex scalars X and Y.

Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `claein` (logical RIGHTV, logical NOINIT, integer N, complex, dimension( ldh, * ) H, integer LDH, complex W, complex, dimension( * ) V, complex, dimension( ldb, * ) B, integer LDB, real, dimension( * ) RWORK, real EPS3, real SMLNUM, integer INFO)

CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. Eigenvector hunting era. 🔥

Purpose:

CLAEIN uses inverse iteration to find a right or left eigenvector corresponding to the eigenvalue W of a complex upper Hessenberg matrix H.

Parameters 💜

RIGHTV
RIGHTV is LOGICAL. Bestie, this is the fork in the road 💅: = .TRUE. → we compute the right eigenvector; = .FALSE. → we compute the left eigenvector.
NOINIT is LOGICAL. Did you bring a starter or not 🫶: = .TRUE. → no initial vector supplied in V; = .FALSE. → initial vector supplied in V.
N is INTEGER. The order of the matrix H. N >= 0.
H is COMPLEX array, dimension (LDH,N). The upper Hessenberg matrix H. She's giving upper Hessenberg, no cap.
LDH is INTEGER. The leading dimension of the array H. LDH >= max(1,N).
W is COMPLEX. The eigenvalue of H whose corresponding right or left eigenvector is to be computed. This is THE girl we're chasing 💖.
V is COMPLEX array, dimension (N). On entry, if NOINIT = .FALSE., V must contain a starting vector for inverse iteration; otherwise V need not be set. On exit, V contains the computed eigenvector, normalized so that the component of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|. Iconic exit behaviour honestly.
B is COMPLEX array, dimension (LDB,N).
LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
RWORK is REAL array, dimension (N).
EPS3 is REAL. A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots. Tiny but she's doing SO much work ✨.
SMLNUM is REAL. A machine-dependent value close to the underflow threshold.
INFO is INTEGER. The vibe report at the end:
- = 0: successful exit 🌈
- = 1: inverse iteration did not converge; V is set to the last iterate.

Not converging is not a fail, bestie, it's just a "we ran out of patience" moment — V still gets the last iterate, so you at least have something to look at. 💀

Author 🎀

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

`subroutine claev2 (complex A, complex B, complex C, real RT1, real RT2, real CS1, complex SN1)` 💜

CLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. Tiny matrix, huge main character energy 👑.

Purpose:

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

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

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

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

Parameters

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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

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

`subroutine clags2 (logical UPPER, real A1, complex A2, real A3, real B1, complex B2, real B3, real CSU, complex SNU, real CSV, complex SNV, real CSQ, complex SNQ)` 🦋

CLAGS2 (that's literally all the tagline says, she's mysterious like that 🙈).

Purpose:

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

         U**H *A*Q = U**H *( A1 A2 )*Q = ( x  0  )
                           ( 0  A3 )     ( x  x  )

and

         V**H*B*Q = V**H *( B1 B2 )*Q = ( x  0  )
                          ( 0  B3 )     ( x  x  )

or if ( .NOT.UPPER ) then

         U**H *A*Q = U**H *( A1 0  )*Q = ( x  x  )
                           ( A2 A3 )     ( 0  x  )

and

         V**H *B*Q = V**H *( B1 0  )*Q = ( x  x  )
                           ( B2 B3 )     ( 0  x  )

where

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

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

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

Parameters

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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

`subroutine clagtm (character TRANS, integer N, integer NRHS, real ALPHA, complex, dimension( * ) DL, complex, dimension( * ) D, complex, dimension( * ) DU, complex, dimension( ldx, * ) X, integer LDX, real BETA, complex, dimension( ldb, * ) B, integer LDB)` 🔥

CLAGTM 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. Look at her, doing algebra with Greek letters and absolutely slaying it ✨.

Purpose:

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

Parameters

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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. The algorithm name alone is a whole aesthetic, fr fr 💅.

Purpose:

CLAHQR is an auxiliary routine called by CHSEQR to update the eigenvalues and Schur decomposition already computed by CHSEQR, by dealing with the Hessenberg submatrix in rows and columns ILO to IHI. She's the supporting actress who eats every scene, honestly.

Parameters

WANTT is LOGICAL. = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.
WANTZ is LOGICAL. = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.
N is INTEGER. The order of the matrix H. N >= 0.
ILO is INTEGER.
IHI
IHI — INTEGER, bestie. The vibe here is: H is already assumed upper triangular in rows and columns IHI+1:N, and H(ILO,ILO-1) = 0 (unless ILO = 1, bestie, no exceptions). CLAHQR mostly does its thing with the Hessenberg submatrix sitting in rows and columns ILO to IHI, but it'll reach out and apply transformations to all of H if WANTT is .TRUE.. The rules: 1 <= ILO <= max(1,IHI) and IHI <= N. 💅
H — COMPLEX array, dimension (LDH,N). On entry, this is the upper Hessenberg matrix H. On exit — and this is where it gets spicy — if INFO is zero and WANTT is .TRUE., then H is upper triangular in rows and columns ILO:IHI. If INFO is zero and WANTT is .FALSE., then the contents of H are unspecified on exit (so don't come crying to her about it, she warned you). The output state of H in case INF is positive? That's explained down below under the INFO description. ✨
LDH — INTEGER. The leading dimension of the array H. LDH >= max(1,N). That's the whole tea. 🫶
W — COMPLEX array, dimension (N). The computed eigenvalues ILO to IHI get stored in the corresponding elements of W. If WANTT is .TRUE., those eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i). Iconic, lowkey. 🌸
ILOZ — INTEGER. That's it. That's the description. Sometimes the manpage just vibes and says nothing. 😭
IHIZ — INTEGER. This one specifies the rows of Z that transformations must be applied to, if WANTZ is .TRUE.. Rules of the road: 1 <= ILOZ <= ILO and IHI <= IHIZ <= N. 🎀
Z — COMPLEX array, dimension (LDZ,N). If WANTZ is .TRUE., then on entry Z must contain the current matrix Z of transformations accumulated by CHSEQR, and on exit Z has been updated. The transformations only get applied to the submatrix Z(ILOZ:IHIZ,ILO:IHI) — not the whole thing, she's precise like that. If WANTZ is .FALSE., Z is not referenced. Not even looked at. Ghosted. 💀
LDZ — INTEGER. The leading dimension of the array Z. LDZ >= max(1,N). 💜
INFO — INTEGER. Okay bestie, here's where the plot thickens:
- = 0: successful exit. We love a happy ending. 🎉
- > 0: if INFO = i, CLAHQR failed to compute all the eigenvalues ILO to IHI in a total of 30 iterations per eigenvalue. In that case, elements i+1:ihi of W contain those eigenvalues which have been successfully computed. So it's a partial slay. 🙈
  - If INFO > 0 and WANTT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix with 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's regardless of the value of WANTT. 👑

Author 🫶

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

Contributors: 🎀

02-96 Based on modifications by David Day, Sandia National Laboratory, USA
12-04 Further modifications by Ralph Byers, University of Kansas, USA

This is a modified version of CLAHQR 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). Glow-up behavior, honestly. ✨

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

CLAHR2 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. Very "I'll do the hard part and hand you the receipts" energy. 💅

Purpose:

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

This is an auxiliary routine called by CGEHRD. iykyk. 🦋

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 — COMPLEX array, dimension (LDA,N-K+1). On entry, the n-by-(n-k+1) general matrix A. On exit, the elements on and above the k-th subdiagonal in the first NB columns are overwritten with the corresponding elements of the reduced matrix; the elements below the k-th subdiagonal, with the array TAU, represent the matrix Q as a product of elementary reflectors. The other columns of A are unchanged. See Further Details for the full lore. 📖
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N). 💅
TAU — COMPLEX array, dimension (NB). The scalar factors of the elementary reflectors. See Further Details. ✨
T — COMPLEX array, dimension (LDT,NB). The upper triangular matrix T. 👑
LDT — INTEGER. The leading dimension of the array T. LDT >= NB. 🎀
Y — COMPLEX array, dimension (LDY,NB). The n-by-nb matrix Y. 🌸
LDY — INTEGER. The leading dimension of the array Y. LDY >= N. 🫶

Author 🫶

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

Further Details: 🦋

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

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

Each H(i) has the form

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

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

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

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

   ( a   a   a   a   a )
   ( a   a   a   a   a )
   ( a   a   a   a   a )
   ( h   h   a   a   a )
   ( v1  h   a   a   a )
   ( v1  v2  a   a   a )
   ( v1  v2  a   a   a )

where a denotes an element of the original matrix A, h denotes a modified element of the upper Hessenberg matrix H, and vi denotes an element of the vector defining H(i). Chef's kiss diagram, honestly. 😭💖

This subroutine is a slight modification of LAPACK-3.0's CLAHRD 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 CLAHRD routine. (This subroutine is not backward compatible with LAPACK-3.0's CLAHRD.) So heads up, bestie — don't assume it behaves like the old one. 🚨

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 claic1 (integer JOB, integer J, complex, dimension( j ) X, real SEST, complex, dimension( j ) W, complex GAMMA, real SESTPR, complex S, complex C)`

CLAIC1 applies one step of incremental condition estimation. Just one step. It's giving "small but mighty". 💅

Purpose:

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

      [ s*x ]
xhat = [  c  ]

is an approximate singular vector of

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

in the sense that

   twonorm(Lhat*xhat) = sestpr.

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

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

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

where alpha = x**H*w. ⭐

Parameters

JOB — 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 — COMPLEX array, dimension (J). The j-vector x. 💜
SEST — REAL. Estimated singular value of j by j matrix L. ✨
W — COMPLEX array, dimension (J). The j-vector w. 🌸
GAMMA — COMPLEX. The diagonal element gamma. 🫶
SESTPR — REAL. Estimated singular value of (j+1) by (j+1) matrix Lhat. 🔥
S — COMPLEX. Sine needed in forming xhat. 🌈
C — COMPLEX. Cosine needed in forming xhat. 🎀

Author 🫶

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

🦋 `clangt`

   real function clangt (character NORM, integer N, complex, dimension( * ) DL,
       complex, dimension( * ) D, complex, dimension( * ) DU)

bestie!! CLANGT hands you back the 1-norm, the Frobenius norm, the infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. She's a norm girlie and she contains multitudes 💅

🌸 Purpose

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

🎀 Returns

   CLANGT = ( 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.

heads up bestie: max(abs(A(i,j))) is NOT a consistent matrix norm. she's built different and the math girlies know it 😭

💜 Parameters

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

🌈 `clanhb`

   real function clanhb (character NORM, character UPLO, integer N, integer K,
       complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * )
       WORK)

CLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix. banded, hermitian, unbothered ✨

🌸 Purpose

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

🎀 Returns

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

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

💜 Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in CLANHB as described above.
UPLO — UPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the band matrix A is supplied.
- = 'U': Upper triangular
- = 'L': Lower triangular
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, CLANHB is set to zero.
K — K is INTEGER. The number of super-diagonals or sub-diagonals of the band matrix A. K >= 0.
AB — AB is COMPLEX array, dimension (LDAB,N). The upper or lower triangle of the hermitian band matrix A, stored in the first K+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB as follows:
- if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
- if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+k).

Note that the imaginary parts of the diagonal elements need not be set and are assumed to be zero.

LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= K+1.
WORK — WORK is REAL 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.

🦋 `clanhp`

   real function clanhp (character NORM, character UPLO, integer N, complex,
       dimension( * ) AP, real, dimension( * ) WORK)

CLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form. packed = no wasted space, we love an efficient queen 🫶

🌸 Purpose

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

🎀 Returns

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

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

💜 Parameters

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

Note that the imaginary parts of the diagonal elements need not be set and are assumed to be zero.

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

🌈 `clanhs`

   real function clanhs (character NORM, integer N, complex, dimension( lda, *
       ) A, integer LDA, real, dimension( * ) WORK)

CLANHS 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 said "everything below the first sub-diagonal is none of my business" and honestly? iconic 💅

🌸 Purpose

CLANHS 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

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

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

💜 Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in CLANHS as described above.
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, CLANHS is set to zero.
A — A is COMPLEX array, dimension (LDA,N). The n by n upper Hessenberg matrix A; the part of A below the first sub-diagonal is not referenced.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(N,1).
WORK — WORK is REAL array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced.

not me double-checking this one 🫣 — for CLANHS it is literally only NORM = 'I' that needs WORK. not '1', not 'O'. just 'I'. read it twice bestie.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🎀 `clanht`

   real function clanht (character NORM, integer N, real, dimension( * ) D,
       complex, dimension( * ) E)

CLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix. tridiagonal but make it hermitian 💜

🌸 Purpose

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

🎀 Returns

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

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

💜 Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in CLANHT as described above.
N

😭 bestie the page literally cuts off right here, mid-parameter. the N entry for CLANHT keeps going — next part has the rest, hold on 🫶

N — N is INTEGER. It's the order of the matrix A, and N >= 0 is non-negotiable bestie. When N = 0, CLANHT gets set to zero, and honestly? iconic behaviour. 💅
D — D is REAL array, dimension (N). These are the diagonal elements of A — the main character entries of the whole matrix, no cap. ✨
E — E is COMPLEX array, dimension (N-1). The (n-1) sub-diagonal or super-diagonal elements of A. the supporting cast, but make it complex 🦋

👑 Author

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

🎀 `real function clansb`

real function clansb (character NORM, character UPLO, integer N, integer K,
    complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * )
    WORK)

CLANSB hands you 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. it's giving "tell me how big this matrix is" energy 🔥

🌟 Purpose

CLANSB 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

CLANSB

CLANSB = ( 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. 🚨 like, it said what it said and it's not sorry.

🍓 Parameters

NORM — NORM is CHARACTER*1. Specifies the value to be returned in CLANSB as described above. you're literally picking the norm, bestie.
UPLO — UPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the band matrix A is supplied.
- = 'U': Upper triangular part is supplied
- = 'L': Lower triangular part is supplied
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, CLANSB is set to zero.
K — K is INTEGER. The number of super-diagonals or sub-diagonals of the band matrix A. K >= 0.
AB — AB is COMPLEX 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 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).
LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= K+1.
WORK — WORK is REAL array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced. she's only here when she needs to be 🙈

👑 Author

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

🌈 `real function clansp`

real function clansp (character NORM, character UPLO, integer N, complex,
    dimension( * ) AP, real, dimension( * ) WORK)

CLANSP 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 = squeezed, and she ate. 💅

🌟 Purpose

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

💫 Returns

CLANSP

CLANSP = ( 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 CLANSP as described above.
UPLO — UPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrix A is supplied.
- = 'U': Upper triangular part of A is supplied
- = 'L': Lower triangular part of A is supplied
N — N is INTEGER. The order of the matrix A. N >= 0. When N = 0, CLANSP is set to zero.
AP — AP is COMPLEX array, dimension (N*(N+1)/2). The upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows:
- if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
- if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
WORK — WORK is REAL array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced. lowkey she's optional 💜

👑 Author

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

🧚 `real function clantb`

real function clantb (character NORM, character UPLO, character DIAG,
    integer N, integer K, complex, dimension( ldab, * ) AB, integer LDAB,
    real, dimension( * ) WORK)

CLANTB 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 AND banded? she's doing the most and we're obsessed. 🔥

🌟 Purpose

CLANTB 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

CLANTB

CLANTB = ( 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 CLANTB as described above.
UPLO — UPLO is CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
DIAG — DIAG is CHARACTER*1. Specifies 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, CLANTB 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 COMPLEX array, dimension (LDAB,N). The upper or lower triangular band matrix A, stored in the first k+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB as follows:
- if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
- if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+k).

Note that when DIAG = 'U', the elements of the array AB corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one. the ones are just... implied. the audacity 💀

LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= K+1.
WORK — WORK is REAL array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced. note: only 'I' gets the WORK treatment here, don't get it twisted 🫶

👑 Author

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

🌸 `real function clantp`

real function clantp (character NORM, character UPLO, character DIAG,
    integer N, complex, dimension( * ) AP, real, dimension( * ) WORK)

CLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. triangular + packed = she's compact and she's powerful. 💖

🌟 Purpose

CLANTP 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

CLANTP

CLANTP = ( 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 CLANTP as described above.
UPLO — UPLO is CHARACTER*1. Specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
DIAG — DIAG is CHARACTER*1. Specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
N — N is INTEGER

N is INTEGER 💖

The order of the matrix A. N >= 0. And when N = 0, CLANTP is set to zero — literally just hands you a zero and moves on. no notes. ✨

AP is COMPLEX array, dimension (N*(N+1)/2) 🦋

The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows:

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

Note that when DIAG = 'U', the elements of the array AP corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one. the way it just skips them and assumes bestie is a 1?? iconic behaviour 💅

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

🦋 `clantr` 🦋

real function clantr (character NORM, character UPLO, character DIAG,
    integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real,
    dimension( * ) WORK)

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

✨ Purpose ✨

CLANTR 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. it's giving "norm girlie, pick your fave" 🎀

💅 Returns 💅

CLANTR = ( 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. 😭 like she's out here doing her own thing fr fr.

🫶 Parameters 🫶

NORM is CHARACTER*1 🌸

Specifies the value to be returned in CLANTR as described above.

UPLO is CHARACTER*1 🌈

Specifies whether the matrix A is upper or lower trapezoidal.

= 'U': Upper trapezoidal
= 'L': Lower trapezoidal

Note that A is triangular instead of trapezoidal if M = N. lowkey a glow-up when they match ✨

DIAG is CHARACTER*1 🎀

Specifies whether or not the matrix A has unit diagonal.

= 'N': Non-unit diagonal
= 'U': Unit diagonal

M is INTEGER 💜

The number of rows of the matrix A. M >= 0, and if UPLO = 'U', M <= N. When M = 0, CLANTR is set to zero. 🫥

N is INTEGER 💜

The number of columns of the matrix A. N >= 0, and if UPLO = 'L', N <= M. When N = 0, CLANTR is set to zero.

A is COMPLEX 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. she's self-confident like that 💅

LDA is INTEGER 💖

The leading dimension of the array A. LDA >= max(M,1).

WORK is REAL array, dimension (MAX(1,LWORK)), where LWORK >= M when NORM = 'I'; otherwise, WORK is not referenced. 🙈

👑 Author 👑

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

🧚 `clapll` 🧚

subroutine clapll (integer N, complex, dimension( * ) X, integer INCX,
    complex, dimension( * ) Y, integer INCY, real SSMIN)

CLAPLL measures the linear dependence of two vectors. it's giving "are u two even that close or nah" 😭

✨ 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. 💫 two steps, one little number, periodt.

🫶 Parameters 🫶

N is INTEGER 🌸

The length of the vectors X and Y.

X is COMPLEX array, dimension (1+(N-1)*INCX) 🦋

On entry, X contains the N-vector X. On exit, X is overwritten. rip bestie, she's gone 💀

INCX is INTEGER ✨

The increment between successive elements of X. INCX > 0.

Y is COMPLEX array, dimension (1+(N-1)*INCY) 🦋

On entry, Y contains the N-vector Y. On exit, Y is overwritten. also gone. not me grieving two vectors at once 😭

INCY is INTEGER ✨

The increment between successive elements of Y. INCY > 0.

SSMIN is REAL 💜

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.

🌈 `clapmr` 🌈

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

CLAPMR rearranges rows of a matrix as specified by a permutation vector. rows only bestie, do not touch the columns 🫶

✨ Purpose ✨

CLAPMR 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 is LOGICAL 🎀

= .TRUE., forward permutation
= .FALSE., backward permutation

M is INTEGER 💜

The number of rows of the matrix X. M >= 0.

N is INTEGER 💜

The number of columns of the matrix X. N >= 0.

X is COMPLEX array, dimension (LDX,N) 🦋

On entry, the M by N matrix X. On exit, X contains the permuted matrix X. ✨ she's changed, but in a good way.

LDX is INTEGER 💖

The leading dimension of the array X, LDX >= MAX(1,M).

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 does the work AND puts everything back?? we love a tidy queen 👑

👑 Author 👑

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

🍓 `clapmt` 🍓

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

CLAPMT performs a forward or backward permutation of the columns of a matrix. yes queen, columns this time 💅

✨ Purpose ✨

CLAPMT 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 is LOGICAL 🎀

= .TRUE., forward permutation
= .FALSE., backward permutation

M is INTEGER 💜

The number of rows of the matrix X. M >= 0.

N is INTEGER 💜

The number of columns of the matrix X. N >= 0.

X is COMPLEX array, dimension (LDX,N) 🦋

On entry, the M by N matrix X. On exit, X contains the permuted matrix X. 🌈

LDX is INTEGER 💖

The leading dimension of the array X, LDX >= MAX(1,M).

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!! so polite, so considerate, the moment 🫶

👑 Author 👑

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

💜 `claqhb` 💜

subroutine claqhb (character UPLO, integer N, integer KD, complex,
    dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) S, real
    SCOND, real AMAX, character EQUED)

CLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ. teamwork makes the dream work ✨

✨ Purpose ✨

CLAQHB equilibrates an Hermitian band matrix A using the scaling factors in the vector S. balance restored, we love to see it 🧘‍♀️

🫶 Parameters 🫶

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 is INTEGER 💜

The order of the matrix A. N >= 0.

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. she's extra in both directions 💅

AB is COMPLEX array, dimension (LDAB,N) 🦋

On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:

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

On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**H *U or A = L*L**H of the band matrix A, in the same storage format as A. ✨ glow-up complete.

LDAB is INTEGER 💖

The leading dimension of the array AB. LDAB >= KD+1.

S

S — S is a REAL array, dimension (N). These are the scale factors for A, bestie. That's it. That's the whole vibe. ✨
SCOND — SCOND is REAL. It's the ratio of the smallest S(i) to the largest S(i). Like, the little guy vs the big guy. 📏
AMAX — AMAX is REAL. Absolute value of the biggest matrix entry. The main character of the matrix, honestly. 👑
EQUED — EQUED is CHARACTER*1. Tells you whether equilibration happened or not, no cap:
- = 'N': No equilibration. We left it alone. 🚫
- = 'Y': Equilibration WAS done, meaning A has been replaced by diag(S) * A * diag(S). Glow-up complete. 💅

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. Simple as that. 🎀
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. Them's the rules. 📐

Author

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

subroutine `claqhp` (`character UPLO`, `integer N`, `complex, dimension( * ) AP`, `real, dimension( * ) S`, `real SCOND`, `real AMAX`, `character EQUED`)

CLAQHP scales a Hermitian matrix stored in packed form. Packed form, bestie — no room for zeros, it's giving minimalist. 🦋

Purpose:

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

Parameters

UPLO — UPLO is CHARACTER*1. Says whether we're storing the upper or lower triangular part of the Hermitian matrix A:
- = 'U': Upper triangular ⬆️
- = 'L': Lower triangular ⬇️
N — N is INTEGER. The order of the matrix A. N >= 0.
AP — AP is COMPLEX array, dimension (N*(N+1)/2). On entry, it's the upper or lower triangle of the Hermitian matrix A, packed columnwise in a linear array. The j-th column of A lives in AP like this:
- if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
- if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

On exit, you get the equilibrated matrix diag(S) * A * diag(S), in the exact same storage format as A. She's transformed but she's still the same girl. 💖

S — S is REAL array, dimension (N). The scale factors for A.
SCOND — SCOND is REAL. Ratio of the smallest S(i) to the largest S(i).
AMAX — AMAX is REAL. 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 `claqp2` (`integer M`, `integer N`, `integer OFFSET`, `complex, dimension( lda, * ) A`, `integer LDA`, `integer, dimension( * ) JPVT`, `complex, dimension( * ) TAU`, `real, dimension( * ) VN1`, `real, dimension( * ) VN2`, `complex, dimension( * ) WORK`)

CLAQP2 computes a QR factorization with column pivoting of the matrix block. Pivoting but make it blocky. 🔥

Purpose:

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

Read that again: pivoted but NOT factorized. She's involved but she's not doing the work. Respect. 💅

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.
A — A is COMPLEX 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. (Say it with me: pivoted, not 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. Character development, honestly. 🌸
TAU — TAU is COMPLEX array, dimension (min(M,N)). The scalar factors of the elementary reflectors.
VN1 — VN1 is REAL array, dimension (N). The vector with the partial column norms.
VN2 — VN2 is REAL array, dimension (N). The vector with the exact column norms. Exact. No rounding slander. 💅
WORK — WORK is COMPLEX 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 `claqps` (`integer M`, `integer N`, `integer OFFSET`, `integer NB`, `integer KB`, `complex, dimension( lda, * ) A`, `integer LDA`, `integer, dimension( * ) JPVT`, `complex, dimension( * ) TAU`, `real, dimension( * ) VN1`, `real, dimension( * ) VN2`, `complex, dimension( * ) AUXV`, `complex, dimension( ldf, * ) F`, `integer LDF`)

CLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. BLAS-3 energy only, we don't do level 1 around here. 🔥

Purpose:

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

In some cases, due to catastrophic cancellations, it cannot
factorize NB columns.  Hence, the actual number of factorized
columns is returned in KB.

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

Okay so real talk: catastrophic cancellations can happen and she simply cannot factorize all NB columns. That's why you get KB back — the actual number she managed. Don't be mad at her, be mad at floating point. 😭

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. (The dream. ✨)
KB — KB is INTEGER. The number of columns actually factorized. (The reality. 🫠)
A — A is COMPLEX 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 COMPLEX array, dimension (KB). The scalar factors of the elementary reflectors.
VN1 — VN1 is REAL array, dimension (N). The vector with the partial column norms.
VN2 — VN2 is REAL array, dimension (N). The vector with the exact column norms.
AUXV — AUXV is COMPLEX array, dimension (NB). Auxiliary vector. The little helper. 🧚
F — F is COMPLEX array, dimension (LDF,NB). Matrix F**H = L * Y**H * 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 `claqr0` (`logical WANTT`, `logical WANTZ`, `integer N`, `integer ILO`, `integer IHI`, `complex, dimension( ldh, * ) H`, `integer LDH`, `complex, dimension( * ) W`, `integer ILOZ`, `integer IHIZ`, `complex, dimension( ldz, * ) Z`, `integer LDZ`, `complex, dimension( * ) WORK`, `integer LWORK`, `integer INFO`)

CLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. Eigenvalues AND maybe the Schur stuff?? She does both, she contains multitudes. 👑

Purpose:

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

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

Parameters

WANTT
WANTT is LOGICAL — bestie, set this to .TRUE. and the full Schur form T is required; set it to .FALSE. and only eigenvalues are required. ✨
WANTZ is LOGICAL — .TRUE. means the matrix of Schur vectors Z is required; .FALSE. means Schur vectors are not required. No cap. 💅
N is INTEGER — the order of the matrix H. N >= 0. 🌸
ILO is INTEGER — that's the whole vibe, literally just that. 🎀
IHI is INTEGER — okay so 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 CGEBAL, and then passed to CGEHRD when the matrix output by CGEBAL 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 is COMPLEX array, dimension (LDH,N) — on entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper triangular matrix T from the Schur decomposition (the Schur form). If INFO = 0 and WANT is .FALSE., then the contents of H are unspecified on exit. (The output value of H when INFO > 0 is given under the description of INFO below.) This subroutine may explicitly set H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. 👑
LDH is INTEGER — the leading dimension of the array H. LDH >= max(1,N). 🌈
W is COMPLEX array, dimension (N) — the computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i). ⭐
ILOZ is INTEGER — just vibes, no extra description. 🫶
IHIZ is INTEGER — specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. 💖
Z is COMPLEX 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 is INTEGER — the leading dimension of the array Z. if WANTZ is .TRUE. then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1. 💀
WORK is COMPLEX 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 LWORK typically as large as 6*N may be required for optimal performance. A workspace query to determine the optimal workspace size is recommended. If LWORK = -1, then CLAQR0 does a workspace query. In this case, CLAQR0 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, CLAQR0 failed to compute all of the eigenvalues. Elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.) If INFO > 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. If INFO > 0 and WANTT is .TRUE., then on exit

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

where U is a unitary matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. If INFO > 0 and WANTZ is .TRUE., then on exit

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

where U is the unitary matrix in (*) (regardless of the value of WANTT.) If INFO > 0 and WANTZ is .FALSE., then Z is not accessed. 😭

💖 Author

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

🫶 Contributors:

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

📚 References:

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

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

🌈 `subroutine claqr1`

subroutine claqr1 (integer N, complex, dimension( ldh, * ) H, integer LDH,
    complex S1, complex S2, complex, dimension( * ) V)

CLAQR1 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, CLAQR1 sets v to a scalar multiple of the first column of the product

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

scaling to avoid overflows and most underflows.

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

🎀 Parameters

N is INTEGER — order of the matrix H. N must be either 2 or 3. 💖
H is COMPLEX array, dimension (LDH,N) — the 2-by-2 or 3-by-3 matrix H in (*). 🌸
LDH is INTEGER — the leading dimension of H as declared in the calling procedure. LDH >= N. 🦋
S1 is COMPLEX — that's it. ✨
S2 is COMPLEX — S1 and S2 are the shifts defining K in (*) above. 💅
V is COMPLEX 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 claqr2`

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

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

💖 Purpose:

CLAQR2 is identical to CLAQR3 except that it avoids recursion by calling CLAHQR instead of CLAQR4.

Aggressive early deflation:

This subroutine accepts as input an upper Hessenberg matrix H and performs an unitary similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been overwritten by a new Hessenberg matrix that is a perturbation of an unitary similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. ✨

🦋 Parameters

WANTT is LOGICAL — if .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues. 👑
WANTZ is LOGICAL — if .TRUE., then the unitary matrix Z is updated so so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced. 💖
N is INTEGER — the order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary matrix Z. 🌸
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 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. ✨
NW is INTEGER — deflation window size. 1 <= NW <= (KBOT-KTOP+1). 💅
H is COMPLEX array, dimension (LDH,N) — on input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries. 🔥
LDH — okay so LDH is an INTEGER, the leading dimension of H just as declared in the calling subroutine. N <= LDH, periodt. 💅
ILOZ — ILOZ is an INTEGER. that's it, that's the tea. 🍓
IHIZ — IHIZ is an INTEGER. it specifies the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N. no cap. ✨
Z — Z is a COMPLEX array, dimension (LDZ,N). if WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. if WANTZ is .FALSE., then Z is unreferenced. literally just ignored, iconic. 💖
LDZ — LDZ is an INTEGER. the leading dimension of Z just as declared in the calling subroutine. 1 <= LDZ.
NS — NS is an INTEGER. the number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine. shifts, shifts, shifts. 🌸
ND — ND is an INTEGER. the number of converged eigenvalues uncovered by this subroutine. we love a converged queen. 👑
SH — SH is a COMPLEX array, dimension (KBOT). on output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT-ND). converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT). (yes, we said SR there, bestie — keeping it exactly as the original, no fixing 💀)
V — V is a COMPLEX array, dimension (LDV,NW). an NW-by-NW work array.
LDV — LDV is an INTEGER. the leading dimension of V just as declared in the calling subroutine. NW <= LDV.
NH — NH is an INTEGER. the number of columns of T. NH >= NW.
T — T is a COMPLEX array, dimension (LDT,NW).
LDT — LDT is an INTEGER. the leading dimension of T just as declared in the calling subroutine. NW <= LDT.
NV — NV is an INTEGER. the number of rows of work array WV available for workspace. NV >= NW.
WV — WV is a COMPLEX array, dimension (LDWV,NW).
LDWV — LDWV is an INTEGER. the leading dimension of W just as declared in the calling subroutine. NW <= LDV. (again, keeping the original's LDV here, no cap 🫶)
WORK — WORK is a COMPLEX 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 workspace glow-up. ✨
LWORK — LWORK is an INTEGER. the dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK.

if LWORK = -1, then a workspace query is assumed; CLAQR2 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. so don't panic, it's just a query, 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 claqr3 🌈

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

CLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). it's giving aggressive early deflation energy. 💅

Purpose:

Aggressive early deflation:

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

Parameters

WANTT — WANTT is LOGICAL. if .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). if .FALSE., then only enough of H is updated to preserve the eigenvalues.
WANTZ — WANTZ is LOGICAL. if .TRUE., then the unitary matrix Z is updated so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). if .FALSE., then Z is not referenced.
N — N is INTEGER. the order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary 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.
NW — NW is INTEGER. deflation window size. 1 <= NW <= (KBOT-KTOP+1).
H — H is a COMPLEX array, dimension (LDH,N). on input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. on output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries.
LDH — LDH is INTEGER. leading dimension of H just as declared in the calling subroutine. N <= LDH.
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 <= IHIZ <= N.
Z — Z is a COMPLEX array, dimension (LDZ,N). IF WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. if WANTZ is .FALSE., then Z is unreferenced.
LDZ — 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.
SH — SH is a COMPLEX array, dimension (KBOT). on output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT-ND). converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT).
V — V is a COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDWV,NW).
LDWV — LDWV is INTEGER. the leading dimension of W just as declared in the calling subroutine. NW <= LDV.
WORK — WORK is a COMPLEX 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.

if LWORK = -1, then a workspace query is assumed; CLAQR3 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.

✨ 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 claqr4 🌈

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

CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. it's giving main character energy. 👑

Purpose:

CLAQR4 does, like, one whole level of recursion for CLAQR0. It's a full, complete, no-notes implementation of the small bulge multi-shift QR algorithm. CLAQR0 can call it, and if the deflation window size is big enough, CLAQR3 can call it too. Here's the tea: this subroutine is identical to CLAQR0 except that it calls CLAQR2 instead of CLAQR3. That's it, that's the diff. 💅

CLAQR4 computes the eigenvalues of a Hessenberg matrix H and, if you want them, the matrices T and Z from the Schur decomposition H = Z T Z**H, where T is an upper triangular matrix (that's the Schur form, bestie) and Z is the unitary matrix of Schur vectors. ✨

Optionally Z may be postmultiplied into an input unitary matrix Q so that this routine can give you the Schur factorization of a matrix A that's already been reduced to the Hessenberg form H by the unitary matrix Q:

A = Q*H*Q**H = (QZ)*H*(QZ)**H

✨ Parameters ✨

WANTT — WANTT is LOGICAL
- = .TRUE. : you want the full Schur form T, no half measures;
- = .FALSE.: you only want eigenvalues, keep it light. 🦋
WANTZ — WANTZ is LOGICAL
- = .TRUE. : you want the matrix of Schur vectors Z;
- = .FALSE.: nope, no Schur vectors required, we're good. 🙈
N — N is INTEGER. The order of the matrix H. N >= 0.
ILO — ILO is INTEGER
IHI — IHI is INTEGER. So the vibe here is: we're assuming H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N, and if ILO > 1, then H(ILO,ILO-1) is zero. Normally ILO and IHI get set by an earlier call to CGEBAL, and then handed over to CGEHRD when the matrix that CGEBAL output gets reduced to Hessenberg form. If that's not your situation, just set ILO and IHI to 1 and N respectively — no drama. If N > 0, then 1 <= ILO <= IHI <= N. And if N = 0, then ILO = 1 and IHI = 0. 👑
H — H is COMPLEX array, dimension (LDH,N).
- On entry: the upper Hessenberg matrix H. 🌸
- On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper triangular matrix T from the Schur decomposition (the Schur form). If INFO = 0 and WANT is .FALSE., then the contents of H are unspecified on exit — literally whatever, don't look at it. (The output value of H when INFO > 0 is described down under INFO below.)
- 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 — LDH is INTEGER. The leading dimension of the array H. LDH >= max(1,N).
W — W is COMPLEX array, dimension (N). The computed eigenvalues of H(ILO:IHI,ILO:IHI) get stored in W(ILO:IHI). And if WANTT is .TRUE., then those eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i). Matching set, iconic. 💖
ILOZ — ILOZ is INTEGER
IHIZ — IHIZ is INTEGER. These two specify the rows of Z that transformations must be applied to if WANTZ is .TRUE.. Rules: 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
Z — Z is COMPLEX array, dimension (LDZ,IHI).
- If WANTZ is .FALSE., then Z is not referenced — she's off duty, not even clocked in. 🫶
- 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 described under INFO below.)
LDZ — 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 COMPLEX array, dimension LWORK. On exit, if LWORK = -1, WORK(1) returns an estimate of the optimal value for LWORK. ✨
LWORK — LWORK is INTEGER. The dimension of the array WORK. LWORK >= max(1,N) is sufficient, but honestly LWORK typically needs to be as large as 6*N for optimal performance — so a workspace query to figure out the optimal workspace size is recommended, bestie, don't skip it. 🌈
- If LWORK = -1, then CLAQR4 does a workspace query. In that case, CLAQR4 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate comes back in WORK(1). No error message related to LWORK is issued by XERBLA. And neither H nor Z are accessed. 🎀
INFO — INFO is INTEGER
- = 0: successful exit, we ate and left no crumbs. 💅
- > 0: if INFO = i, CLAQR4 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, so don't spiral. 🍓)
- If INFO > 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H.
- If INFO > 0 and WANTT is .TRUE., then on exit:

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

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

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

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

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

If INFO > 0 and WANTZ is .FALSE., then Z is not accessed. Not even a little. 🙈

🌟 Author 🌟

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

🫶 Contributors 🫶

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

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

🧚 `claqr5` — the small-bulge multi-shift QR sweep 💫

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

CLAQR5 performs a single small-bulge multi-shift QR sweep. That's the whole job description, and she does it flawlessly. 🔥

Purpose:

CLAQR5 called by CLAQR0 performs a single small-bulge multi-shift QR sweep.

Parameters

WANTT — WANTT is LOGICAL. WANTT = .true. if the triangular Schur factor is being computed. WANTT is set to .false. otherwise.
WANTZ — WANTZ is LOGICAL. WANTZ = .true. if the unitary Schur factor is being computed. WANTZ is set to .false. otherwise.
KACC22 — KACC22 is INTEGER with value 0, 1, or 2. Specifies the computation mode of far-from-diagonal orthogonal updates.
- = 0: CLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries.
- = 1: CLAQR5 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 that is no longer supported — it's giving deprecated, she retired. 💀
N — N is INTEGER. N is the order of the Hessenberg matrix H upon which this subroutine operates.
KTOP — KTOP is INTEGER
KBOT — 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.

(Yes bestie, no check. It's on you. 🫣)

NSHFTS — NSHFTS is INTEGER. NSHFTS gives the number of simultaneous shifts. NSHFTS must be positive and even — those are the rules, no cap. ⭐
S — S is COMPLEX array, dimension (NSHFTS). S contains the shifts of origin that define the multi-shift QR sweep. On output S may be reordered. 🔄
H — H is COMPLEX 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 — LDH is INTEGER. LDH is the leading dimension of H just as declared in the calling procedure. LDH >= MAX(1,N).
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 <= IHIZ <= N
Z — Z is COMPLEX array, dimension (LDZ,IHIZ).
- If WANTZ = .TRUE., then the QR Sweep unitary similarity transformation is accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
- If WANTZ = .FALSE., then Z is unreferenced. 🌸
LDZ — *LDZ is INTEGER*… and bestie, this one's story is literally continued in the next part, stay tuned 🫶
LDZ is INTEGER — bestie, LDZ is the leading dimension of Z, exactly like you declared it back in the calling procedure. And it's gotta be LDZ >= N. 💅
V is COMPLEX array, dimension (LDV,NSHFTS/2) 🌸
LDV is INTEGER — the leading dimension of V as declared in the calling procedure. LDV >= 3. ✨
U is COMPLEX array, dimension (LDU,2*NSHFTS) 🦋
LDU is INTEGER — the leading dimension of U as declared in the calling subroutine. LDU >= 2*NSHFTS. 💜
NV is INTEGER — NV is the number of rows in WV available for workspace. NV >= 1. 🎀
WV is COMPLEX array, dimension (LDWV,2*NSHFTS) 🌈
LDWV is INTEGER — the leading dimension of WV as declared in the calling subroutine. LDWV >= NV. ⭐
NH is INTEGER — NH is the number of columns in array WH available for workspace. NH >= 1. 🍓
WH is COMPLEX array, dimension (LDWH,NH) 💖
LDWH is INTEGER — leading dimension of WH just as declared in the calling procedure. LDWH >= 2*NSHFTS. 🫶

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 claqsb`

subroutine claqsb (character UPLO, integer N, integer KD, complex,
    dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) S, real
    SCOND, real AMAX, character EQUED)

CLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. like... okay so it takes the scaling factors somebody else already did the math for and applies them. teamwork, we love to see it. 🫶

Purpose:

CLAQSB equilibrates a symmetric band matrix A using the scaling factors in the vector S. ✨

Parameters

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 — is INTEGER. The order of the matrix A. N >= 0. 💅
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 — is COMPLEX array, dimension (LDAB,N).
- On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:
  - if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  - if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
- On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**H *U or A = L*L**H of the band matrix A, in the same storage format as A. 🦋
LDAB — is INTEGER. The leading dimension of the array AB. LDAB >= KD+1. ⭐
S — is REAL array, dimension (N). The scale factors for A. 🎀
SCOND — is REAL. Ratio of the smallest S(i) to the largest S(i). 💜
AMAX — is REAL. Absolute value of largest matrix entry. 🔥
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 claqsp`

subroutine claqsp (character UPLO, integer N, complex, dimension( * ) AP,
    real, dimension( * ) S, real SCOND, real AMAX, character EQUED)

CLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. packed storage besties rise up, the whole matrix is squished into one linear array and honestly? iconic. 🌸

Purpose:

CLAQSP equilibrates a symmetric matrix A using the scaling factors in the vector S. ✨

Parameters

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 — is INTEGER. The order of the matrix A. N >= 0. 💅
AP — is COMPLEX 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. 🦋
S — is REAL array, dimension (N). The scale factors for A. 🎀
SCOND — is REAL. Ratio of the smallest S(i) to the largest S(i). 💜
AMAX — is REAL. Absolute value of largest matrix entry. 🔥
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 clar1v`

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

CLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. not the whole inverse, just one column, just the one that matters. lowkey genius behaviour. 🧠✨

Purpose:

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

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

Parameters

N — is INTEGER. The order of the matrix L D L**T. 💅
B1 — is INTEGER. First index of the submatrix of L D L**T. 🌸
BN — is INTEGER. Last index of the submatrix of L D L**T. 🦋
LAMBDA — is REAL. The shift. In order to compute an accurate eigenvector, LAMBDA should be a good approximation to an eigenvalue of L D L**T. 🎯
L — is REAL array, dimension (N-1). The (n-1) subdiagonal elements of the unit bidiagonal matrix L, in elements 1 to N-1. ✨
D — is REAL array, dimension (N). The n diagonal elements of the diagonal matrix D. 💜
LD — is REAL array, dimension (N-1). The n-1 elements L(i)*D(i). 🎀
LLD — is REAL array, dimension (N-1). The n-1 elements L(i)*L(i)*D(i). 🍓
PIVMIN — is REAL. The minimum pivot in the Sturm sequence. 🔥
GAPTOL — is REAL. Tolerance that indicates when eigenvector entries are negligible w.r.t. their contribution to the residual. 👀
Z — is COMPLEX array, dimension (N).
- On input, all entries of Z must be set to 0. fr fr, all of them, zero them out first. 😭
- On output, Z contains the (scaled) r-th column of the inverse. The scaling is such that Z(R) equals 1. 🌈
WANTNC — is LOGICAL. Specifies whether NEGCNT has to be computed. 🤔
NEGCNT

🦋 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 REAL. The square of the 2-norm of Z.

🦋 MINGMA 🦋

MINGMA is REAL. 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 REAL. NRMINV = 1/SQRT( ZTZ )

🦋 RESID 🦋

RESID is REAL. The residual of the FP vector. RESID = ABS( MINGMA )/SQRT( ZTZ )

🦋 RQCORR 🦋

RQCORR is REAL. The Rayleigh Quotient correction to LAMBDA. RQCORR = MINGMA*TMP

🦋 WORK 🦋

WORK is REAL 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 clar2v` 🌈

subroutine clar2v (integer N, complex, dimension( * ) X, complex, dimension(
* ) Y, complex, dimension( * ) Z, integer INCX, real, dimension( * ) C,
complex, dimension( * ) S, integer INCC)

CLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. bestie. it does the both sides thing. iconic.

✨ Purpose ✨

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

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

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

✨ Parameters ✨

N — N is INTEGER. The number of plane rotations to be applied.
X — X is COMPLEX array, dimension (1+(N-1)*INCX). The vector x; the elements of x are assumed to be real.
Y — Y is COMPLEX array, dimension (1+(N-1)*INCX). The vector y; the elements of y are assumed to be real.
Z — Z is COMPLEX 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 REAL array, dimension (1+(N-1)*INCC). The cosines of the plane rotations.
S — S is COMPLEX 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 clarcm` 🫶

subroutine clarcm (integer M, integer N, real, dimension( lda, * ) A,
integer LDA, complex, dimension( ldb, * ) B, integer LDB, complex,
dimension( ldc, * ) C, integer LDC, real, dimension( * ) RWORK)

CLARCM copies all or part of a real two-dimensional array to a complex array. lowkey the simplest lil multiply in the whole building.

🌸 Purpose 🌸

CLARCM performs a very simple matrix-matrix multiplication:

         C := A * B,

where A is M by M and real; B is M by N and complex; C is M by N and complex.

🌸 Parameters 🌸

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

🌸 Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💅 `subroutine clarf` 💅

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

CLARF applies an elementary reflector to a general rectangular matrix. yes queen, go reflect.

⭐ Purpose ⭐

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

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

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

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

To apply H**H (the conjugate transpose of H), supply conjg(tau) instead tau.

⭐ Parameters ⭐

SIDE — SIDE is CHARACTER*1
- = 'L': form H * C
- = 'R': form C * H
M — M is INTEGER. The number of rows of the matrix C.
N — N is INTEGER. The number of columns of the matrix C.
V — V is COMPLEX 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 COMPLEX. The value tau in the representation of H.
C — C is COMPLEX array, dimension (LDC,N). On entry, the M-by-N matrix C. On exit, C is overwritten by the matrix H * C if SIDE = 'L', or C * H if SIDE = 'R'.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
WORK — WORK is COMPLEX 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 clarfb` 🎀

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

CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix. block reflector era. we love to see it 💖

🧚 Purpose 🧚

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

🧚 Parameters 🧚

SIDE — SIDE is CHARACTER*1
- = 'L': apply H or H**H from the Left
- = 'R': apply H or H**H from the Right
TRANS — TRANS is CHARACTER*1
- = 'N': apply H (No transpose)
- = 'C': apply H**H (Conjugate transpose)
DIRECT — DIRECT is CHARACTER*1. Indicates how H is formed from a product of elementary reflectors
- = 'F': H = H(1) H(2) . . . H(k) (Forward)
- = 'B': H = H(k) . . . H(2) H(1) (Backward)
STOREV — STOREV is CHARACTER*1. Indicates how the vectors which 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 (= the number of elementary reflectors whose product defines the block reflector). If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0.
V — V is COMPLEX 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.

LDV — LDV is INTEGER. The leading dimension of the array V. If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
T — T is COMPLEX 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 COMPLEX array, dimension (LDC,N). On entry, the M-by-N matrix C. On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
LDC — and that's literally where the page cuts off, bestie 💀 mid-parameter, no cap. next part picks it up.
LDC is INTEGER — the leading dimension of the array C. LDC >= max(1,M).
WORK is COMPLEX array, dimension (LDWORK,K).
LDWORK is INTEGER — the leading dimension of the array WORK. 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.

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

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

CLARFB_GETT — bestie, say it with me 🫶

✨ Purpose: ✨

CLARFB_GETT applies a complex Householder block reflector H from the left to a complex (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 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 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 order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. 0 <= K <= N.
T — T is COMPLEX array, dimension (LDT,K). The upper-triangular K-by-K matrix T in the representation of the block reflector.
LDT — LDT is INTEGER. The leading dimension of the array T. LDT >= K.
A — A is COMPLEX array, dimension (LDA,N).

On entry:

In the K-by-N upper-trapezoidal part A: input matrix A.
In the columns below the diagonal: columns of V1 (ones are not stored on the diagonal).

On exit: A is overwritten by rectangular K-by-N product H*A.

See Further Details section.

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

On entry:

In the M-by-(N-K) right block: input matrix B.
In the M-by-N left block: columns of V2.

On exit: B is overwritten by rectangular M-by-N product H*B.

See Further Details section.

LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M).
WORK — WORK is COMPLEX array, dimension (LDWORK,max(K,N-K)).
LDWORK — LDWORK is INTEGER. The leading dimension of the array WORK. LDWORK>=max(1,K).

🎀 Author 🎀

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 Contributors: 🌈

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

💜 Further Details: 💜

(1) Description of the Algebraic Operation.

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

A = ( A1, A2 ).

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

B = ( B1, B2 ).

Perform the operation:

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

On input:

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

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

where the column blocks are:

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

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

V = ( V1 )
    ( V2 )

where:

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

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

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

On output:

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

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

where the column blocks are:

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

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

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

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

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

If IDENT != 'I':

The computation for column block 1:

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

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

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

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

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

If IDENT == 'I':

The operation for column block 1:

A1_out: = A1_in - V1*T*A1_in

B1_out: = - V2*T*A1_in

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

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

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

(2) Description of the Algorithmic Computation.

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

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

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

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

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

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

okay so when V1 is literally just the identity matrix, you get to skip a whole bunch of steps, and the glow-up is real. here's the receipts:

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

🌈 the combined algorithm (cases (A) + (B), merged, iconic) 🌈

Combine these cases (A) and (B) together, this is the resulting algorithm — and bestie, this is the one that actually runs:

               if ( N > K ) then

                 (First Step - column block 2)

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

               else

               (Second Step - column block 1)

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

               end if

✨ `subroutine clarfg` ✨

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

CLARFG generates an elementary reflector (Householder matrix).

🎀 Purpose 🎀

okay so CLARFG generates a complex elementary reflector H of order n, such that

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

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

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

where tau is a complex scalar and v is a complex (n-1)-element vector. Note that H is not hermitian — like, do not come for me assuming it is, it isn't 💅

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

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

🫶 Parameters 🫶

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

👑 Author 👑

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

💜 `subroutine clarfgp` 💜

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

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

🌸 Purpose 🌸

CLARFGP generates a complex elementary reflector H of order n, such that

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

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

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

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

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

🍓 Parameters 🍓

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

⭐ Author ⭐

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

🔥 `subroutine clarft` 🔥

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

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

✨ Purpose ✨

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

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

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

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

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

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

🧚 Parameters 🧚

DIRECT — DIRECT is CHARACTER*1. Specifies the order in which the elementary reflectors are multiplied to form the block reflector:
- = 'F': H = H(1) H(2) . . . H(k) (Forward)
- = 'B': H = H(k) . . . H(2) H(1) (Backward)
STOREV — STOREV is CHARACTER*1. Specifies how the vectors which define the elementary reflectors are stored (see also Further Details):
- = 'C': columnwise
- = 'R': rowwise
N — N is INTEGER. The order of the block reflector H. N >= 0.
K — K is INTEGER. The order of the triangular factor T (= the number of elementary reflectors). K >= 1.
V — V is COMPLEX array, dimension
- (LDV,K) if STOREV = 'C'
- (LDV,N) if STOREV = 'R'

The matrix V. See further details.

LDV — LDV is INTEGER. The leading dimension of the array V. If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
TAU — TAU is COMPLEX array, dimension (K). TAU(i) must contain the scalar factor of the elementary reflector H(i).
T — T is COMPLEX array, dimension (LDT,K). The k by k triangular factor T of the block reflector. If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is lower triangular. The rest of the array is not used.
LDT — LDT is INTEGER. The leading dimension of the array T. LDT >= K.

🌈 Author 🌈

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

🎀 Further Details: 🎀

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

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

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

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

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

🦋 `subroutine clarfx` 🦋

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

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

💫 Purpose 💫

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

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

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

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

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

🌸 Parameters 🌸

SIDE — bestie hold that thought, this one finishes in the next part 💖
SIDE is CHARACTER*1
- = 'L': form H * C
- = 'R': form C * H
M

M is INTEGER

The number of rows of the matrix C.

N

N is INTEGER

The number of columns of the matrix C.

V

V is COMPLEX array, dimension (M) if SIDE = 'L' or (N) if SIDE = 'R'

The vector v in the representation of H.

TAU

TAU is COMPLEX

The value tau in the representation of H.

C

C is COMPLEX array, dimension (LDC,N)

On entry, the m by n matrix C.

On exit, C is overwritten by the matrix H * C if SIDE = 'L', or C * H if SIDE = 'R'.

LDC

LDC is INTEGER

The leading dimension of the array C. LDC >= max(1,M).

WORK

WORK is COMPLEX array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R'

WORK is not referenced if H has order < 11. (it literally just sits there, unbothered, when H is smol 💅)

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLARFY ✨

🌸 Purpose:

CLARFY applies an elementary reflector, or Householder matrix, H, to an n x n Hermitian 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. (free pass, no work today bestie)

💖 Parameters

UPLO

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the Hermitian 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 COMPLEX 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. (do NOT be out here passing zero, that's not the vibe 😭)

TAU

TAU is COMPLEX

The value tau as described above.

C

C is COMPLEX 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 COMPLEX array, dimension (N)

🎀 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 subroutine clargv (integer N, complex, dimension( * ) X, integer INCX, complex, dimension( * ) Y, integer INCY, real, dimension( * ) C, integer INCC)

CLARGV generates a vector of plane rotations with real cosines and complex sines. (cosines stay real, sines get to be complex, it's giving main character energy)

🌈 Purpose:

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

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

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

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

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

✨ Parameters

N

N is INTEGER

The number of plane rotations to be generated.

X

X is COMPLEX array, dimension (1+(N-1)*INCX)

On entry, the vector x.

On exit, x(i) is overwritten by r(i), for i = 1,...,n.

INCX

INCX is INTEGER

The increment between elements of X. INCX > 0.

Y

Y is COMPLEX 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 REAL 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.

🌸 Further Details:

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

This version has a few statements commented out for thread safety (machine parameters are computed on each entry). 10 feb 03, SJH. (thread safety said "I got this", and honestly? it did)

⭐ subroutine clarnv (integer IDIST, integer, dimension( 4 ) ISEED, integer N, complex, dimension( * ) X)

CLARNV returns a vector of random numbers from a uniform or normal distribution. ✨ rng era ✨

🎀 Purpose:

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

🦋 Parameters

IDIST

IDIST is INTEGER

Specifies the distribution of the random numbers:

= 1: real and imaginary parts each uniform (0,1)
= 2: real and imaginary parts each uniform (-1,1)
= 3: real and imaginary parts each normal (0,1)
= 4: uniformly distributed on the disc abs(z) < 1
= 5: uniformly distributed on the circle abs(z) = 1

ISEED

ISEED is INTEGER array, dimension (4)

On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. (odd!! not even, bestie, the whole thing breaks otherwise 💀)

On exit, the seed is updated.

N

N is INTEGER

The number of random numbers to be generated.

X

X is COMPLEX array, dimension (N)

The generated random numbers.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💖 Further Details:

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

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

CLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. (big one. huge even. iconic behavior)

💅 Purpose:

CLARRV 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 SLARRE. (should have been — as in, that's not optional, that's the assignment)

🌸 Parameters

N

N is INTEGER

The order of the matrix. N >= 0.

VL

VL is REAL

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 REAL

Upper bound of the interval that contains the desired eigenvalues. VL < VU. Needed to compute gaps on the left or right end of the extremal eigenvalues in the desired RANGE.

D

D is REAL array, dimension (N)

On entry, the N diagonal elements of the diagonal matrix D.

On exit, D may be overwritten.

L

L is REAL 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 SLARRE.

On exit, L is overwritten.

PIVMIN

PIVMIN is REAL

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 REAL

RTOL1

RTOL1 is REAL

RTOL2

RTOL2 is REAL

Parameters for bisection. An interval [LEFT,RIGHT] has converged if

RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )

W
W — bestie, this is a REAL array, dimension (N). The first M elements are holding the APPROXIMATE eigenvalues for which eigenvectors are to be computed. You gotta group those eigenvalues by split-off block and order them from smallest to largest within the block (yes, the output array W from SLARRE is what's expected here, no cap). 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. Periodt.
WERR — WERR is REAL array, dimension (N). The first M elements contain the semiwidth of the uncertainty interval of the corresponding eigenvalue in W. 💅
WGAP — WGAP is REAL array, dimension (N). The separation from the right neighbor eigenvalue in W.
IBLOCK — IBLOCK is 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 — INDEXW is 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 — GERS is REAL 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 — Z is COMPLEX 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 — LDZ is INTEGER. The leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', LDZ >= max(1,N).
ISUPPZ — ISUPPZ is 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 — WORK is REAL array, dimension (12*N).
IWORK — IWORK is INTEGER array, dimension (7*N).
INFO — INFO is INTEGER. Here's the breakdown, bestie:
- = 0: successful exit
- > 0: A problem occurred in CLARRV.
- < 0: One of the called subroutines signaled an internal problem. Needs inspection of the corresponding parameter IINFO for further information.
- =-1: Problem in SLARRB when refining a child's eigenvalues.
- =-2: Problem in SLARRF when computing the RRR of a child. When a child is inside a tight cluster, it can be difficult to find an RRR. A partial remedy from the user's point of view is to make the parameter MINRGP smaller and recompile. However, as the orthogonality of the computed vectors is proportional to 1/MINRGP, the user should be aware that he might be trading in precision when he decreases MINRGP.
- =-3: Problem in SLARRB 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 clartv ✨

   subroutine clartv (integer N, complex, dimension( * ) X, integer INCX,
       complex, dimension( * ) Y, integer INCY, real, dimension( * ) C,
       complex, dimension( * ) S, integer INCC)

CLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a pair of vectors. It's giving glow-up for your vectors. ✨

Purpose:

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

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

Parameters:

N — N is INTEGER. The number of plane rotations to be applied.
X — X is COMPLEX array, dimension (1+(N-1)*INCX). The vector x.
INCX — INCX is INTEGER. The increment between elements of X. INCX > 0.
Y — Y is COMPLEX array, dimension (1+(N-1)*INCY). The vector y.
INCY — INCY is INTEGER. The increment between elements of Y. INCY > 0.
C — C is REAL array, dimension (1+(N-1)*INCC). The cosines of the plane rotations.
S — S is COMPLEX 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 clascl ✨

   subroutine clascl (character TYPE, integer KL, integer KU, real CFROM, real
       CTO, integer M, integer N, complex, dimension( lda, * ) A, integer LDA,
       integer INFO)

CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. No cap.

Purpose:

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

Parameters:

TYPE — TYPE is CHARACTER*1. TYPE indices the storage type of the input matrix.
- = 'G': A is a full matrix.
- = 'L': A is a lower triangular matrix.
- = 'U': A is an upper triangular matrix.
- = 'H': A is an upper Hessenberg matrix.
- = 'B': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the lower half stored.
- = 'Q': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the upper half stored.
- = 'Z': A is a band matrix with lower bandwidth KL and upper bandwidth KU. See CGBTRF for storage details.
KL — KL is INTEGER. The lower bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'.
KU — KU is INTEGER. The upper bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'.
CFROM — CFROM is REAL.
CTO — CTO is REAL. The matrix A is multiplied by CTO/CFROM. A(I,J) is computed without over/underflow if the final result CTO*A(I,J)/CFROM can be represented without over/underflow. CFROM must be nonzero.
M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
A — A is COMPLEX array, dimension (LDA,N). The matrix to be multiplied by CTO/CFROM. See TYPE for the storage type.
LDA — LDA is INTEGER. The leading dimension of the array A. Here are the rules:

If TYPE = 'G', 'L', 'U', 'H', LDA >= max(1,M);
   TYPE = 'B', LDA >= KL+1;
   TYPE = 'Q', LDA >= KU+1;
   TYPE = 'Z', LDA >= 2*KL+KU+1.

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

✨ Author ✨

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

✨ subroutine claset ✨

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

CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values. It's like setting the vibe of your matrix. 💅

Purpose:

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

Parameters:

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

✨ Author ✨

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

🎀 `clasr` — the plane-rotation glow-up 🎀

subroutine clasr (character SIDE, character PIVOT, character DIRECT, integer
    M, integer N, real, dimension( * ) C, real, dimension( * ) S, complex,
    dimension( lda, * ) A, integer LDA)

CLASR applies a sequence of plane rotations to a general rectangular matrix. bestie. 💅

✨ Purpose ✨

okay so CLASR takes a whole sequence of real plane rotations and applies them to a complex matrix A — from either the left or the right, whatever you're feeling. 🦋

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

   A := P*A

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

   A := A*P**T

where P is an orthogonal matrix consisting of a sequence of z plane rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', and P**T is the transpose of P. no cap. 💖

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

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

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

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

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

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

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

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

where R(k) appears as a rank-2 modification to the identity matrix in rows and columns k and k+1. the way it just slides in there 😭

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

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

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

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

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

where R(k) appears in rows and columns k and z. and here's the iconic part: the rotations are performed without ever forming P(k) explicitly. not me being efficient 💀

🌸 Parameters 🌸

SIDE — SIDE is CHARACTER*1. specifies whether the plane rotation matrix P is applied to A on the left or the right.
- = 'L': Left, compute A := P*A
- = 'R': Right, compute A:= A*P**T
PIVOT — PIVOT is CHARACTER*1. specifies the plane for which P(k) is a plane rotation matrix.
- = 'V': Variable pivot, the plane (k,k+1)
- = 'T': Top pivot, the plane (1,k+1)
- = 'B': Bottom pivot, the plane (k,z)
DIRECT — DIRECT is CHARACTER*1. specifies whether P is a forward or backward sequence of plane rotations.
- = 'F': Forward, P = P(z-1)*...*P(2)*P(1)
- = 'B': Backward, P = P(1)*P(2)*...*P(z-1)
M — M is INTEGER. the number of rows of the matrix A. if m <= 1, an immediate return is effected. (we love an early exit 👑)
N — N is INTEGER. the number of columns of the matrix A. if n <= 1, an immediate return is effected.
C — C is REAL array, dimension
- (M-1) if SIDE = 'L'
- (N-1) if SIDE = 'R'

the cosines c(k) of the plane rotations. 💖

S — S is REAL array, dimension
- (M-1) if SIDE = 'L'
- (N-1) if SIDE = 'R'

the sines s(k) of the plane rotations. the 2-by-2 plane rotation part of the matrix P(k), R(k), has the form

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

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

👑 Author 👑

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

🫶 `claswp` — row swap central 🫶

subroutine claswp (integer N, complex, dimension( lda, * ) A, integer LDA,
    integer K1, integer K2, integer, dimension( * ) IPIV, integer INCX)

CLASWP performs a series of row interchanges on a general rectangular matrix. 🔥

✨ Purpose ✨

CLASWP performs a series of row interchanges on the matrix A. one row interchange is initiated for each of rows K1 through K2 of A. literally just swapping rows, bestie, one at a time. 💅

🌸 Parameters 🌸

N — N is INTEGER. the number of columns of the matrix A.
A — A is COMPLEX 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 — LDA is INTEGER. the leading dimension of the array A.
K1 — K1 is INTEGER. the first element of IPIV for which a row interchange will be done.
K2 — K2 is INTEGER. (K2-K1+1) is the number of elements of IPIV for which a row interchange will be done.
IPIV — IPIV is INTEGER array, dimension (K1+(K2-K1)*abs(INCX)). the vector of pivot indices. only the elements in positions K1 through K1+(K2-K1)*abs(INCX) of IPIV are accessed. IPIV(K1+(K-K1)*abs(INCX)) = L implies rows K and L are to be interchanged. 🎀
INCX — INCX is INTEGER. the increment between successive values of IPIV. if INCX is negative, the pivots are applied in reverse order. (yes, walking backwards, main character energy 🙈)

👑 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

💜 `clatbs` — triangular banded solver, but make it safe 💜

subroutine clatbs (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, integer KD, complex, dimension( ldab, * )
    AB, integer LDAB, complex, dimension( * ) X, real SCALE, real,
    dimension( * ) CNORM, integer INFO)

CLATBS solves a triangular banded system of equations. no cap. 🌈

✨ Purpose ✨

CLATBS solves one of the triangular systems

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

with scaling to prevent overflow, where A is an upper or lower triangular band matrix. here A**T denotes the transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. 🌸

if the unscaled problem will not cause overflow, the Level 2 BLAS routine CTBSV 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. 💀 (yes bestie, that's a real branch, read it again)

🌸 Parameters 🌸

UPLO — UPLO is CHARACTER*1. specifies whether the matrix A is upper or lower triangular.
- = 'U': Upper triangular
- = 'L': Lower triangular
TRANS — TRANS is CHARACTER*1. specifies the operation applied to A.
- = 'N': Solve A * x = s*b (No transpose)
- = 'T': Solve A**T * x = s*b (Transpose)
- = 'C': Solve A**H * x = s*b (Conjugate transpose)
DIAG — DIAG is CHARACTER*1. specifies whether or not the matrix A is unit triangular.
- = 'N': Non-unit triangular
- = 'U': Unit triangular
NORMIN — NORMIN is CHARACTER*1. specifies whether CNORM has been set or not.
- = 'Y': CNORM contains the column norms on entry
- = 'N': CNORM is not set on entry. on exit, the norms will be computed and stored in CNORM.
N — N is INTEGER. the order of the matrix A. N >= 0.
KD — KD is INTEGER. the number of subdiagonals or superdiagonals in the triangular matrix A. KD >= 0.
AB — AB is COMPLEX 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 COMPLEX 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 REAL. the scaling factor s for the triangular system
- A * x = s*b, A**T * x = s*b, or A**H * x = s*b.

if SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0. 🚨

CNORM

CNORM is a REAL array, dimension (N) 🫶

if NORMIN = 'Y', then CNORM is an input argument, and CNORM(j) holds the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) has to be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) has to be greater than or equal to the 1-norm. like, no negotiating on this one bestie 💅
if NORMIN = 'N', then 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, CTBSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every single operation. she is thorough, okay?? 👑

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 CTBSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow). math girlies stay winning 🧚

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. not the singular-matrix jump scare 😭

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

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

We simultaneously compute two bounds

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

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

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

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

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

subroutine clatdf 👑

subroutine clatdf (integer IJOB, integer N, complex, dimension( ldz, * ) Z,
    integer LDZ, complex, dimension( * ) RHS, real RDSUM, real RDSCAL,
    integer, dimension( * ) IPIV, integer, dimension( * ) JPIV)

CLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. dif-estimate era 💅

Purpose: ✨

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

The factorization of Z returned by CGETC2 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 CGECON, 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 COMPLEX array, dimension (LDZ, N)

On entry, the LU part of the factorization of the n-by-n matrix Z computed by CGETC2: Z = P * L * U * Q

LDZ

LDZ is INTEGER

The leading dimension of the array Z. LDA >= max(1, N).

RHS

RHS is COMPLEX 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 REAL

On entry, the sum of squares of computed contributions to the Dif-estimate under computation by CTGSYL, 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 CTGSY2 is called by CTGSYL.

RDSCAL

RDSCAL is REAL

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 CTGSY2 is called by CTGSYL.

IPIV

IPIV is INTEGER array, dimension (N).

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

JPIV

JPIV is INTEGER array, dimension (N).

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

Author 🦋

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details: 💜

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

Contributors: 🌸

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

References: ⭐

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

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

subroutine clatps 💖

subroutine clatps (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, complex, dimension( * ) AP, complex,
    dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)

CLATPS solves a triangular system of equations with the matrix held in packed storage. packed storage girlies rise 🎀

Purpose: ✨

CLATPS solves one of the triangular systems

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

with scaling to prevent overflow, where A is an upper or lower triangular matrix stored in packed form. Here A**T denotes the transpose of A, A**H denotes the conjugate transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine CTPSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned. she said "no overflow on my watch" 💅

Parameters 🎀

UPLO

UPLO is CHARACTER*1

Specifies whether the matrix A is upper or lower triangular.

= 'U': Upper triangular
= 'L': Lower triangular

TRANS

TRANS is CHARACTER*1

Specifies the operation applied to A.

= 'N': Solve A * x = s*b (No transpose)
= 'T': Solve A**T * x = s*b (Transpose)
= 'C': Solve A**H * x = s*b (Conjugate transpose)

DIAG

DIAG is CHARACTER*1

Specifies whether or not the matrix A is unit triangular.

= 'N': Non-unit triangular
= 'U': Unit triangular

NORMIN

NORMIN is CHARACTER*1

Specifies whether CNORM has been set or not.

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

N

N is INTEGER

The order of the matrix A. N >= 0.

AP

AP is COMPLEX array, dimension (N*(N+1)/2)

The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows:

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

X

X is COMPLEX 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 — REAL, and it is literally the scaling factor s for the triangular system

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

If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0. 💀 so if that zero shows up, bestie, it's not a bug, it's a vibe: you're holding a non-trivial answer.

CNORM — REAL array, dimension (N)

If NORMIN = 'Y', then CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. Two flavors, okay?

If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm ✨
If TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm ✨

If NORMIN = 'N', then 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 🎉 slay
- < 0: if INFO = -k, the k-th argument had an illegal value 😭 somebody passed something cursed

Author 🌸

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

Further Details 💅

Okay so here's the tea: a rough bound on x gets computed. If that bound is less than overflow, CTPSV is called. Otherwise? Specific code is used that checks for possible overflow or divide-by-zero at every single operation. No cap. 🫶

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 CTPSV 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, then x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. 👀 plot twist.

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

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

We simultaneously compute two bounds

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

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

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

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

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

subroutine `clatrd`

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

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

Purpose 🎀

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

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

This is an auxiliary routine called by CHETRD. iykyk 💅

Parameters ✨

UPLO — CHARACTER*1. Specifies whether the upper or lower triangular part of the Hermitian 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 — COMPLEX array, dimension (LDA,N). On entry, the Hermitian matrix A. If UPLO = 'U', the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. 🙈

On exit:

if UPLO = 'U', the last NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements above the diagonal, with the array TAU, represent the unitary matrix Q as a product of elementary reflectors;
if UPLO = 'L', the first NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements below the diagonal, with the array TAU, represent the unitary matrix Q as a product of elementary reflectors.

See Further Details. 👀

LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N).
E — 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 — COMPLEX 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 — COMPLEX array, dimension (LDW,NB). The n-by-nb matrix W required to update the unreduced part of A.
LDW — INTEGER. The leading dimension of the array W. LDW >= max(1,N).

Author 🌸

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

Further Details 💅

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

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

Each H(i) has the form

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

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

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

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

Each H(i) has the form

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

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

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

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

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

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

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

subroutine `clatrs`

subroutine clatrs (character UPLO, character TRANS, character DIAG,
    character NORMIN, integer N, complex, dimension( lda, * ) A, integer
    LDA, complex, dimension( * ) X, real SCALE, real, dimension( * ) CNORM,
    integer INFO)

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

Purpose 🎀

CLATRS solves one of the triangular systems

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

with scaling to prevent overflow. Here A is an upper or lower triangular matrix, A**T denotes the transpose of A, A**H denotes the conjugate transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine CTRSV 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 — 👀 hold that thought bestie, this one's spilling into the next part…
TRANS — CHARACTER*1. Specifies the operation applied to A. bestie, pick your fighter ✨
- 'N': Solve A * x = s*b (No transpose)
- 'T': Solve A**T * x = s*b (Transpose)
- 'C': Solve A**H * x = s*b (Conjugate transpose)
DIAG — CHARACTER*1. Specifies whether or not the matrix A is unit triangular. the tea is: is the diagonal just vibes or actual numbers?
- 'N': Non-unit triangular
- 'U': Unit triangular
NORMIN — CHARACTER*1. Specifies whether CNORM has been set or not. did you bring the norms or are we computing them?
- '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. no negative energy allowed 💅
A — 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. so iconic of her to skip the diagonal 🎀
LDA — INTEGER. The leading dimension of the array A. LDA >= max (1,N). respect the leading dimension or else 😭
X — COMPLEX array, dimension (N). On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x. glow up: from b to x 🦋
SCALE — REAL. The scaling factor s for the triangular system

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

If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0. if this is zero, we have a problem bestie 💀

CNORM — 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. norms on norms on norms 👑
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, CTRSV 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 CTRSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).

The bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. If the computed bound is greater than a large constant, x is scaled to prevent overflow, but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.

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

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

We simultaneously compute two bounds

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

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

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

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

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

`subroutine clauu2 (character UPLO, integer N, complex, dimension( lda, * ) A, integer LDA, integer INFO)`

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

Purpose

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

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

This is the unblocked form of the algorithm, calling Level 2 BLAS. level 2 BLAS? we love a reliable bestie ✨

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. no negative order, periodt
A — COMPLEX array, dimension (LDA,N). On entry, the triangular factor U or L. On exit, if UPLO = 'U', the upper triangle of A is overwritten with the upper triangle of the product U * U**H; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**H * L. she's changing in public, and that's fine 🦋
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N). mind the gap (the leading dimension) 😭
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 clauum (character UPLO, integer N, complex, dimension( lda, * ) A, integer LDA, integer INFO)`

CLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). blocked algorithm? highkey iconic 👑

Purpose

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

This is the blocked form of the algorithm, calling Level 3 BLAS. level 3 BLAS, the main character energy we deserve 🔥

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. same energy as before, still no negatives
A — COMPLEX array, dimension (LDA,N). On entry, the triangular factor U or L. On exit, if UPLO = 'U', the upper triangle of A is overwritten with the upper triangle of the product U * U**H; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**H * L. overwriting? she's transforming, not copying 🦋
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,N). leading dimension stays winning 💅
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 crot (integer N, complex, dimension( * ) CX, integer INCX, complex, dimension( * ) CY, integer INCY, real C, complex S)`

CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors. rotation? bestie, we're spinning 💫

Purpose

CROT applies a plane rotation, where the cos (C) is real and the sin (S) is complex, and the vectors CX and CY are complex. real cos, complex sin, a duo for the ages ✨

Parameters

N — INTEGER. The number of elements in the vectors CX and CY. count 'em
CX — COMPLEX array, dimension (N). On input, the vector X. On output, CX is overwritten with C*X + S*Y. glow up: X becomes C*X + S*Y 🦋
INCX — INTEGER. The increment between successive values of CX. INCX <> 0. zero increment? not on my watch 💀
CY — COMPLEX array, dimension (N). On input, the vector Y. On output, CY is overwritten with -CONJG(S)*X + C*Y. conjugate? we love a little twist 🌀
INCY — INTEGER. The increment between successive values of CY. INCX <> 0. yes, it says INCX here, we're keeping it exactly 🙈
C — REAL. the real one
S —

                     S is COMPLEX
                     C and S define a rotation
                        [  C          S  ]
                        [ -conjg(S)   C  ]
                     where C*C + S*CONJG(S) = 1.0.

so yeah — S is COMPLEX, and C plus S together define a rotation, and they gotta satisfy C*C + S*CONJG(S) = 1.0. that's the vibe, no notes 💅

✨ Author ✨

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

subroutine cspmv (character UPLO, integer N, complex ALPHA, complex,
    dimension( * ) AP, complex, dimension( * ) X, integer INCX, complex
    BETA, complex, dimension( * ) Y, integer INCY)

CSPMV computes a matrix-vector product for complex vectors using a complex symmetric packed matrix. okay so this one's a little workhorse 🐴

🎯 Purpose:

CSPMV performs the matrix-vector operation

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

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

📋 Parameters

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

Unchanged on exit. 🫶

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

✨ Author ✨

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

subroutine cspr (character UPLO, integer N, complex ALPHA, complex,
    dimension( * ) X, integer INCX, complex, dimension( * ) AP)

CSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix. she's giving ✨glow-up✨ energy, fr fr

🎯 Purpose:

CSPR performs the symmetric rank 1 operation

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

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

📋 Parameters

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

Unchanged on exit.

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

⚠️ bestie listen: note that the imaginary parts of the diagonal elements need not be set, they are assumed to be zero, and on exit they are set to zero. 😭 don't come for me, that's just how it is

✨ Author ✨

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

subroutine csrscl (integer N, real SA, complex, dimension( * ) SX, integer
    INCX)

CSRSCL multiplies a vector by the reciprocal of a real scalar. tiny function, huge heart 🎀

🎯 Purpose:

CSRSCL multiplies an n-element complex 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. so like, she's careful. she's mindful. we love her 💜

📋 Parameters

N — N is INTEGER. The number of components of the vector x.
SA — SA is REAL. The scalar a which is used to divide each component of x.

🚨 SA must be >= 0, or the subroutine will divide by zero. that's not a vibe check you wanna fail 💀

SX — SX is COMPLEX array, dimension (1+(N-1)*abs(INCX)). The n-element vector x.
INCX — INCX is INTEGER. The increment between successive values of the vector SX.
- > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n

✨ Author ✨

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

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

CTPRFB applies a complex 'triangular-pentagonal' block reflector to a complex matrix, which is composed of two blocks. she's the main character of block reflectors, no cap 👑

🎯 Purpose:

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

📋 Parameters

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

The pentagonal matrix V, which contains the elementary reflectors H(1), H(2), ..., H(K). See Further Details. 🌸

LDV
T — a COMPLEX array, dimension (LDT,K). It's the triangular K-by-K matrix T chillin' inside the block-reflector representation. That's the whole vibe of it.
LDT — INTEGER, bestie. The leading dimension of T. You need LDT >= K, no exceptions, periodt.
A — a COMPLEX 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 gets overwritten by the corresponding block of H*C or H**H*C or C*H or C*H**H. See Further Details, babes.
LDA — INTEGER. Leading dimension of A. If SIDE = 'L', then LDA >= max(1,K); if SIDE = 'R', then LDA >= max(1,M). Simple as.
B — a COMPLEX array, dimension (LDB,N). On entry, the M-by-N matrix B. On exit, B is overwritten by the corresponding block of H*C or H**H*C or C*H or C*H**H. Same deal, see Further Details 💖
LDB — INTEGER. Leading dimension of B. LDB >= max(1,M).
WORK — a COMPLEX array, dimension (LDWORK,N) if SIDE = 'L', and (LDWORK,K) if SIDE = 'R'. She's the workspace, she's doing her best.
LDWORK — INTEGER. Leading dimension of 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.

Iconic lineup ngl 💅

🌈 Further Details:

Okay so this is the part where they explain what C even is, and it's actually kind of gorgeous once you see it:

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.

Like?? The blocks just assemble. It's giving LEGO but make it linear algebra 🧚

✨ `integer function icmax1 (integer N, complex, dimension(*) CX, integer INCX)`

ICMAX1 finds the index of the first vector element of maximum absolute value. She's a finder, she's a looker, she's not playing games.

🎯 Purpose:

ICMAX1 finds the index of the first vector element of maximum absolute value.

Based on ICAMAX from Level 1 BLAS.
The change is to use the 'genuine' absolute value.

Note the glow-up: it's based on ICAMAX from Level 1 BLAS, but the twist is it uses the 'genuine' absolute value. Different flavour, better result 💅

🎀 Parameters

N — INTEGER. How many elements are in the vector CX.
CX — a COMPLEX array, dimension (N). This is the vector CX. The ICMAX1 function returns the index of its first element of maximum absolute value.
INCX — INTEGER. The spacing between successive values of CX. INCX >= 1. Bestie said no zero strides today.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 Contributors:

Nick Higham for use with CLACON.

✨ `integer function ilaclc (integer M, integer N, complex, dimension( lda, * ) A, integer LDA)`

ILACLC scans a matrix for its last non-zero column. She's the one who checks the back of the fridge and tells you what's still there.

🎯 Purpose:

ILACLC scans A for its last non-zero column.

🎀 Parameters

M — INTEGER. The number of rows of the matrix A.
N — INTEGER. The number of columns of the matrix A.
A — a COMPLEX array, dimension (LDA,N). The m by n matrix A. That's it, that's the girl.
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M).

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ `integer function ilaclr (integer M, integer N, complex, dimension( lda, * ) A, integer LDA)`

ILACLR scans a matrix for its last non-zero row. Same energy, other axis. The way that she just... goes down the rows and finds the last one with something in it 🦋

🎯 Purpose:

ILACLR scans A for its last non-zero row.

🎀 Parameters

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

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ `integer function izmax1 (integer N, complex16, dimension() ZX, integer INCX)`

IZMAX1 finds the index of the first vector element of maximum absolute value. Yes it's basically ICMAX1 but the double-precision twin — sisters, not twins, but honestly twins.

🎯 Purpose:

IZMAX1 finds the index of the first vector element of maximum absolute value.

Based on IZAMAX from Level 1 BLAS.
The change is to use the 'genuine' absolute value.

Again with the glow-up: based on IZAMAX from Level 1 BLAS, and the change is that it uses the 'genuine' absolute value. Consistency is her brand ✨

🎀 Parameters

N — INTEGER. The number of elements in the vector ZX.
ZX — a COMPLEX*16 array, dimension (N). The vector ZX. The IZMAX1 function returns the index of its first element of maximum absolute value.
INCX — INTEGER. The spacing between successive values of ZX. INCX >= 1.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 Contributors:

Nick Higham for use with ZLACON.

✨ `real function scsum1 (integer N, complex, dimension( * ) CX, integer INCX)`

SCSUM1 forms the 1-norm of the complex vector using the true absolute value. She's out here summing, she's out here norming, she's thriving.

🎯 Purpose:

SCSUM1 takes the sum of the absolute values of a complex
vector and returns a single precision result.

Based on SCASUM from the Level 1 BLAS.
The change is to use the 'genuine' absolute value.

Based on SCASUM from the Level 1 BLAS, and yeah, the change is the 'genuine' absolute value. We love a consistent queen 🙈

🎀 Parameters

N — INTEGER. The number of elements in the vector CX.
CX — a COMPLEX array, dimension (N). The vector whose elements will be summed.
INCX — INTEGER. The spacing between successive values of CX. INCX > 0. Strictly greater than zero this time, not >= 1 — read it carefully bestie.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 Contributors:

Nick Higham for use with CLACON.

👑 Author

Generated automatically by Doxygen for LAPACK from the source code.

And that's the whole page, babes 💖 This one's LAPACK, version 3.11.0, dated Sun Nov 27 2022, and the man page is complexOTHERauxiliary(3) — you made it to the end of all 25 parts, which is genuinely main character energy 🌈 Go drink some water, 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

complexOTHERauxiliary(3)

complexOTHERauxiliary(3) — LAPACK manual, bestie 💖

✨ NAME ✨

✨ SYNOPSIS ✨

✨ Functions ✨

✨ Detailed Description ✨

🎀 Function Documentation 🎀

subroutine clabrd (integer M, integer N, integer NB, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( * ) TAUQ, complex, dimension( * ) TAUP, complex, dimension( ldx, * ) X, integer LDX, complex, dimension( ldy, * ) Y, integer LDY)

Purpose

Parameters

Author

Further Details

subroutine clacgv (integer N, complex, dimension( * ) X, integer INCX)

Purpose

Parameters

Author

subroutine clacn2 (integer N, complex, dimension( * ) V, complex, dimension( * ) X, real EST, integer KASE, integer, dimension( 3 ) ISAVE)

Purpose

Parameters

Author 💖

Further Details: 👑

Contributors: 🫶

References: 📚

subroutine clacon (integer N, complex, dimension( n ) V, complex, dimension( n ) X, real EST, integer KASE)

subroutine clacp2 (character UPLO, integer M, integer N, real, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)

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

subroutine clacrm (integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, complex, dimension( ldc, * ) C, integer LDC, real, dimension( * ) RWORK)

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

complex function cladiv (complex X, complex Y)

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

Author 🎀

subroutine claev2 (complex A, complex B, complex C, real RT1, real RT2, real CS1, complex SN1) 💜

subroutine clags2 (logical UPPER, real A1, complex A2, real A3, real B1, complex B2, real B3, real CSU, complex SNU, real CSV, complex SNV, real CSQ, complex SNQ) 🦋

subroutine clagtm (character TRANS, integer N, integer NRHS, real ALPHA, complex, dimension( * ) DL, complex, dimension( * ) D, complex, dimension( * ) DU, complex, dimension( ldx, * ) X, integer LDX, real BETA, complex, dimension( ldb, * ) B, integer LDB) 🔥

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

Author 🫶

Contributors: 🎀

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

Purpose:

Parameters

Author 🫶

Further Details: 🦋

References: 📚

subroutine claic1 (integer JOB, integer J, complex, dimension( j ) X, real SEST, complex, dimension( j ) W, complex GAMMA, real SESTPR, complex S, complex C)

Purpose:

Parameters

Author 🫶

🦋 clangt

🌸 Purpose

🎀 Returns

💜 Parameters

👑 Author

🌈 clanhb

🌸 Purpose

🎀 Returns

💜 Parameters

👑 Author

🦋 clanhp

🌸 Purpose

🎀 Returns

💜 Parameters

👑 Author

🌈 clanhs

🌸 Purpose

🎀 Returns

💜 Parameters

👑 Author

🎀 clanht

🌸 Purpose

🎀 Returns

💜 Parameters

👑 Author

🎀 real function clansb

🌟 Purpose

💫 Returns

🍓 Parameters

👑 Author

🌈 real function clansp

🌟 Purpose

💫 Returns

`subroutine clabrd` `(integer M, integer N, integer NB, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( * ) TAUQ, complex, dimension( * ) TAUP, complex, dimension( ldx, * ) X, integer LDX, complex, dimension( ldy, * ) Y, integer LDY)`

`subroutine clacgv` `(integer N, complex, dimension( * ) X, integer INCX)`

`subroutine clacn2` `(integer N, complex, dimension( * ) V, complex, dimension( * ) X, real EST, integer KASE, integer, dimension( 3 ) ISAVE)`

subroutine `clacon` (integer N, complex, dimension( n ) V, complex, dimension( n ) X, real EST, integer KASE)

subroutine `clacp2` (character UPLO, integer M, integer N, real, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)

subroutine `clacpy` (character UPLO, integer M, integer N, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)

subroutine `clacrm` (integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, complex, dimension( ldc, * ) C, integer LDC, real, dimension( * ) RWORK)

subroutine `clacrt` (integer N, complex, dimension( * ) CX, integer INCX, complex, dimension( * ) CY, integer INCY, complex C, complex S)

complex function `cladiv` (complex X, complex Y)

subroutine `claein` (logical RIGHTV, logical NOINIT, integer N, complex, dimension( ldh, * ) H, integer LDH, complex W, complex, dimension( * ) V, complex, dimension( ldb, * ) B, integer LDB, real, dimension( * ) RWORK, real EPS3, real SMLNUM, integer INFO)

`subroutine claev2 (complex A, complex B, complex C, real RT1, real RT2, real CS1, complex SN1)` 💜

`subroutine clags2 (logical UPPER, real A1, complex A2, real A3, real B1, complex B2, real B3, real CSU, complex SNU, real CSV, complex SNV, real CSQ, complex SNQ)` 🦋

`subroutine clagtm (character TRANS, integer N, integer NRHS, real ALPHA, complex, dimension( * ) DL, complex, dimension( * ) D, complex, dimension( * ) DU, complex, dimension( ldx, * ) X, integer LDX, real BETA, complex, dimension( ldb, * ) B, integer LDB)` 🔥

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

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

`subroutine claic1 (integer JOB, integer J, complex, dimension( j ) X, real SEST, complex, dimension( j ) W, complex GAMMA, real SESTPR, complex S, complex C)`

🦋 `clangt`

🌈 `clanhb`

🦋 `clanhp`

🌈 `clanhs`

🎀 `clanht`

🎀 `real function clansb`

🌈 `real function clansp`

🧚 `real function clantb`

🌸 `real function clantp`

🦋 `clantr` 🦋

🧚 `clapll` 🧚

🌈 `clapmr` 🌈

🍓 `clapmt` 🍓

💜 `claqhb` 💜

subroutine `claqhp` (`character UPLO`, `integer N`, `complex, dimension( * ) AP`, `real, dimension( * ) S`, `real SCOND`, `real AMAX`, `character EQUED`)

subroutine `claqp2` (`integer M`, `integer N`, `integer OFFSET`, `complex, dimension( lda, * ) A`, `integer LDA`, `integer, dimension( * ) JPVT`, `complex, dimension( * ) TAU`, `real, dimension( * ) VN1`, `real, dimension( * ) VN2`, `complex, dimension( * ) WORK`)

🌈 `subroutine claqr1`

🔥 `subroutine claqr2`