auxOTHERcomputational(3)

📦

version 3.11.0 · 2 versions

3.11.0 · debian bookworm (latest)
3.10.0 · ubuntu jammy

💅✨ auxOTHERcomputational(3) ✨💅

LAPACK — and yes bestie, this is the page for auxOTHERcomputational, aka the auxiliary Computational routines. It's a whole friend group of them. This page is long, so it's been chopped into 21 parts, and you're holding part 1. Buckle up. 🎀

🌸 NAME 🌸

auxOTHERcomputational — the auxiliary Computational routines. That's the name, that's the vibe, we move. 💫

✨ SYNOPSIS ✨

🦋 Functions 🦋

okay so here's the full guest list for this party. Every single one of them, no skipping, no "oh we'll get to it later." 👑

character *1 function chla_transtype (TRANS) — CHLA_TRANSTYPE
subroutine dbdsdc (UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO) — DBDSDC
subroutine dbdsqr (UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO) — DBDSQR
subroutine ddisna (JOB, M, N, D, SEP, INFO) — DDISNA
subroutine dlaed0 (ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO) — DLAED0 used by DSTEDC. Computes all the eigenvalues and the matching eigenvectors of an unreduced symmetric tridiagonal matrix, using the divide and conquer method. slay. 🪄
subroutine dlaed1 (N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO) — DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after it gets modified by a rank-one symmetric matrix. This one shows up when the original matrix is tridiagonal. 💅
subroutine dlaed2 (K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO) — DLAED2 used by DSTEDC. Merges eigenvalues and deflates the secular equation. Used when the original matrix is tridiagonal. ✨
subroutine dlaed3 (K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO) — DLAED3 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal. 🌸
subroutine dlaed4 (N, I, D, Z, DELTA, RHO, DLAM, INFO) — DLAED4 used by DSTEDC. Finds a single root of the secular equation. one (1) root, bestie. 🎯
subroutine dlaed5 (I, D, Z, DELTA, RHO, DLAM) — DLAED5 used by DSTEDC. Solves the 2-by-2 secular equation. smol but mighty. 🐣
subroutine dlaed6 (KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO) — DLAED6 used by DSTEDC. Computes one Newton step in the solution of the secular equation. ✏️
subroutine dlaed7 (ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO) — DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This time it's for when the original matrix is dense. 🧠
subroutine dlaed8 (ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO) — DLAED8 used by DSTEDC. Merges eigenvalues and deflates the secular equation. Used when the original matrix is dense. 💜
subroutine dlaed9 (K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO) — DLAED9 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense. 🌈
subroutine dlaeda (N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO) — DLAEDA used by DSTEDC. Computes the Z vector that determines the rank-one modification of the diagonal matrix. Used when the original matrix is dense. 🦋
subroutine dlagtf (N, A, LAMBDA, B, C, TOL, D, IN, INFO) — DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges. 🔥
subroutine dlamrg (N1, N2, A, DTRD1, DTRD2, INDEX) — DLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order. the ultimate group project. 🤝
subroutine dlartgs (X, Y, SIGMA, CS, SN) — DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem. (yes it's called a bulge, we don't make the rules 💀)
subroutine dlasq1 (N, D, E, WORK, INFO) — DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr. ⭐
subroutine dlasq2 (N, Z, INFO) — DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z, to high relative accuracy. Used by sbdsqr and sstegr. 👑
subroutine dlasq3 (I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU) — DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr. 🌀
subroutine dlasq4 (I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G) — DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. 🔍
subroutine dlasq5 (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS) — DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. 🏓
subroutine dlasq6 (I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2) — DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. 🏓
subroutine dlasrt (ID, N, D, INFO) — DLASRT sorts numbers in increasing or decreasing order. 🍓
subroutine dstebz (RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO) — DSTEBZ
subroutine dstedc (COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO) — DSTEDC
subroutine dsteqr (COMPZ, N, D, E, Z, LDZ, WORK, INFO) — DSTEQR
subroutine dsterf (N, D, E, INFO) — DSTERF
integer function iladiag (DIAG) — ILADIAG
integer function ilaprec (PREC) — ILAPREC
integer function ilatrans (TRANS) — ILATRANS
integer function ilauplo (UPLO) — ILAUPLO
subroutine sbdsdc (UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO) — SBDSDC
subroutine sbdsqr (UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO) — SBDSQR
subroutine sdisna (JOB, M, N, D, SEP, INFO) — SDISNA
subroutine slaed0 (ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO) — SLAED0 used by SSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method. 🪄
subroutine slaed1 (N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO) — SLAED1 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal. 💅
subroutine slaed2 (K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO) — SLAED2 used by SSTEDC. Merges eigenvalues and deflates the secular equation. Used when the original matrix is tridiagonal. ✨
subroutine slaed3 (K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO) — SLAED3 used by SSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal. 🌸
subroutine slaed4 (N, I, D, Z, DELTA, RHO, DLAM, INFO) — SLAED4 used by SSTEDC. Finds a single root of the secular equation. 🎯
subroutine slaed5 (I, D, Z, DELTA, RHO, DLAM) — SLAED5 used by SSTEDC. Solves the 2-by-2 secular equation. 🐣
subroutine slaed6 (KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO) — SLAED6 used by SSTEDC. Computes one Newton step in solution of the secular equation. ✏️
subroutine slaed7 (ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO) — SLAED7 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense. 🧠
subroutine slaed8 (ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO) — SLAED8 used by SSTEDC. Merges eigenvalues and deflates the secular equation. Used when the original matrix is dense. 💜
subroutine slaed9 (K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO) — SLAED9 used by SSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense. 🌈
subroutine slaeda (N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO) — SLAEDA used by SSTEDC. Computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense. 🦋
subroutine slagtf (N, A, LAMBDA, B, C, TOL, D, IN, INFO) — SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges. 🔥
subroutine slamrg (N1, N2, A, STRD1, STRD2, INDEX) — SLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order. 🤝
subroutine slartgs (X, Y, SIGMA, CS, SN) — SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem. 💀
subroutine slasq1 (N, D, E, WORK, INFO) — SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr. ⭐
subroutine slasq2 (N, Z, INFO) — SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr. 👑
subroutine slasq3 (I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU) — SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr. 🌀
subroutine slasq4 (I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G) — SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. 🔍
subroutine slasq5 (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS) — SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. 🏓
subroutine slasq6 (I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2) — SLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. 🏓
subroutine slasrt (ID, N, D, INFO) — SLASRT sorts numbers in increasing or decreasing order. 🍓
subroutine spttrf (N, D, E, INFO) — SPTTRF
subroutine sstebz (RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO) — SSTEBZ
subroutine sstedc (COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO) — SSTEDC
subroutine ssteqr (COMPZ, N, D, E, Z, LDZ, WORK, INFO) — SSTEQR
subroutine ssterf (N, D, E, INFO) — SSTERF

🫶 Detailed Description 🫶

This is the group of auxiliary Computational routines. That's the whole description. Short, sweet, iconic, no notes. 💖

📚 Function Documentation 📚

🌟 character*1 function chla_transtype (integer TRANS) 🌟

CHLA_TRANSTYPE

🎯 Purpose:

okay so bestie, when the integer that got passed in isn't one of the transposition-operator constants, CHLA_TRANSTYPE just hits you with the letter 'X'. That's the "nah, that's not a transposition operator" flag, lowkey iconic of it. ✨

Otherwise — if it is a legit transposition op — CHLA_TRANSTYPE returns the constant value that matches TRANS. That's the whole deal. No cap. 🎀

🌸 Author 🌸

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

🦋 subroutine `dbdsdc`

subroutine dbdsdc (character UPLO, character COMPQ, integer N, double
    precision, dimension( * ) D, double precision, dimension( * ) E, double
    precision, dimension( ldu, * ) U, integer LDU, double precision,
    dimension( ldvt, * ) VT, integer LDVT, double precision, dimension( * )
    Q, integer, dimension( * ) IQ, double precision, dimension( * ) WORK,
    integer, dimension( * ) IWORK, integer INFO)

DBDSDC

💖 Purpose 💖

DBDSDC computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B:

B = U * S * VT

using a divide and conquer method, where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices of left and right singular vectors, respectively. DBDSDC can be used to compute all singular values, and optionally, singular vectors or singular vectors in compact form.

This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See DLASD3 for details. 🙈

The code currently calls DLASDQ if singular values only are desired. However, it can be slightly modified to compute singular values using the divide and conquer method.

✨ Parameters ✨

UPLO — UPLO is CHARACTER*1
- = 'U': B is upper bidiagonal.
- = 'L': B is lower bidiagonal.
COMPQ — COMPQ is CHARACTER*1. Specifies whether singular vectors are to be computed as follows:
- = 'N': Compute singular values only;
- = 'P': Compute singular values and compute singular vectors in compact form;
- = 'I': Compute singular values and singular vectors.
N — N is INTEGER. The order of the matrix B. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B.
E — E is DOUBLE PRECISION array, dimension (N-1). On entry, the elements of E contain the offdiagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E has been destroyed. 💀
U — U is DOUBLE PRECISION array, dimension (LDU,N). If COMPQ = 'I', then: On exit, if INFO = 0, U contains the left singular vectors of the bidiagonal matrix. For other values of COMPQ, U is not referenced.
LDU — LDU is INTEGER. The leading dimension of the array U. LDU >= 1. If singular vectors are desired, then LDU >= max( 1, N ).
VT — VT is DOUBLE PRECISION array, dimension (LDVT,N). If COMPQ = 'I', then: On exit, if INFO = 0, VT**T contains the right singular vectors of the bidiagonal matrix. For other values of COMPQ, VT is not referenced.
LDVT — LDVT is INTEGER. The leading dimension of the array VT. LDVT >= 1. If singular vectors are desired, then LDVT >= max( 1, N ).
Q — Q is DOUBLE PRECISION array, dimension (LDQ). If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, Q contains all the DOUBLE PRECISION data in

LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))

words of memory, where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about 25). For other values of COMPQ, Q is not referenced.

IQ — IQ is INTEGER array, dimension (LDIQ). If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, IQ contains all INTEGER data in

LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))

WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
- If COMPQ = 'N' then LWORK >= (4 * N).
- If COMPQ = 'P' then LWORK >= (6 * N).
- If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
IWORK — IWORK is INTEGER array, dimension (8*N)
INFO — INFO is INTEGER
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: The algorithm failed to compute a singular value. The update process of divide and conquer failed. 😭

🌈 Author 🌈

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

🫶 Contributors: 🫶

Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

⭐ subroutine `dbdsqr` ⭐

subroutine dbdsqr (character UPLO, integer N, integer NCVT, integer NRU,
    integer NCC, double precision, dimension( * ) D, double precision,
    dimension( * ) E, double precision, dimension( ldvt, * ) VT, integer
    LDVT, double precision, dimension( ldu, * ) U, integer LDU, double
    precision, dimension( ldc, * ) C, integer LDC, double precision,
    dimension( * ) WORK, integer INFO)

DBDSQR

💜 Purpose 💜

DBDSQR computes the singular values and, optionally, the right and/or left singular vectors from the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B using the implicit zero-shift QR algorithm. The SVD of B has the form

B = Q * S * P**T

where S is the diagonal matrix of singular values, Q is an orthogonal matrix of left singular vectors, and P is an orthogonal matrix of right singular vectors. If left singular vectors are requested, this subroutine actually returns U*Q instead of Q, and, if right singular vectors are requested, this subroutine returns P**T*VT instead of P**T, for given real input matrices U and VT. When U and VT are the orthogonal matrices that reduce a general matrix A to bidiagonal form:

A = U*B*VT

as computed by DGEBRD, then

A = (U*Q) * S * (P**T*VT)

is the SVD of A. Optionally, the subroutine may also compute Q**T*C for a given real input matrix C.

See 'Computing Small Singular Values of Bidiagonal Matrices With Guaranteed High Relative Accuracy,' by J. Demmel and W. Kahan, LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, no. 5, pp. 873-912, Sept 1990) and 'Accurate singular values and differential qd algorithms,' by B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics Department, University of California at Berkeley, July 1992 for a detailed description of the algorithm. 📚

✨ Parameters ✨

UPLO — UPLO is CHARACTER*1
- = 'U': B is upper bidiagonal;
- = 'L': B is lower bidiagonal.
N — N is INTEGER. The order of the matrix B. N >= 0.
NCVT — NCVT is INTEGER. The number of columns of the matrix VT. NCVT >= 0.
NRU — NRU is INTEGER. The number of rows of the matrix U. NRU >= 0.
NCC — NCC is INTEGER. The number of columns of the matrix C. NCC >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B in decreasing order.
E — E is DOUBLE PRECISION array, dimension (N-1). On entry, the N-1 offdiagonal elements of the bidiagonal matrix B. On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E will contain the diagonal and superdiagonal elements of a bidiagonal matrix orthogonally equivalent to the one given as input.
VT — VT is DOUBLE PRECISION array, dimension (LDVT, NCVT). On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P**T * VT. Not referenced if NCVT = 0.
LDVT — LDVT is INTEGER. The leading dimension of the array VT. LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
U — U is DOUBLE PRECISION array, dimension (LDU, N). On entry, an NRU-by-N matrix U. On exit, U is overwritten by U * Q. Not referenced if NRU = 0.
LDU — LDU is INTEGER. The leading dimension of the array U. LDU >= max(1,NRU).
C — C is DOUBLE PRECISION array, dimension (LDC, NCC). On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q**T * C. Not referenced if NCC = 0.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
WORK — WORK is DOUBLE PRECISION array, dimension (4*(N-1))
INFO — and that's where this part of the page cuts off, bestie — the rest of INFO rolls in with the next bit 🍓
> 0:
- if NCVT = NRU = NCC = 0,
  - = 1, a split was marked by a positive value in E
  - = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop)
  - = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)
- else NCVT = NRU = NCC = 0, the algorithm did not converge; D and E contain the elements of a bidiagonal matrix which is orthogonally similar to the input matrix B; if INFO = i, i elements of E have not converged to zero.

Like, bestie, if you see a positive INFO here it is NOT a cute lil victory lap. It means the algorithm ✨gave up✨. Read the number, it's telling you exactly which vibe check it failed 💀

🔧 Internal Parameters:

TOLMUL — DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
- TOLMUL controls the convergence criterion of the QR loop.
- If it is positive, TOLMUL*EPS is the desired relative precision in the computed singular values.
- If it is negative, abs(TOLMUL*EPS*sigma_max) is the desired absolute accuracy in the computed singular values (this corresponds to relative accuracy abs(TOLMUL*EPS) in the largest singular value).
- abs(TOLMUL) should be between 1 and 1/EPS, and preferably between 10 (for fast convergence) and .1/EPS (so that there's actually some accuracy in the results).
- Default is to lose at either one eighth or 2 of the available decimal digits in each computed singular value (whichever is smaller).
MAXITR — INTEGER, default = 6
- MAXITR controls the maximum number of passes of the algorithm through its inner loop. The algorithms stops (and so fails to converge) if the number of passes through the inner loop exceeds MAXITR*N**2.

So basically TOLMUL is the "how picky are we" knob and MAXITR is the "how long are we willing to wait" knob. Turn them wrong and it's giving never-converging energy 🫠

📝 Note:

Bug report from Cezary Dendek. On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is removed since it can overflow pretty easily (for N larger or equal than 18,919). We instead use MAXITDIVN = MAXITR*N.

Shoutout Cezary for catching that 🙏 integer overflow is the jumpscare of numerical computing, no cap.

✍️ Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine ddisna (character JOB, integer M, integer N, double precision, dimension( * ) D, double precision, dimension( * ) SEP, integer INFO)

DDISNA

🎯 Purpose:

DDISNA computes the reciprocal condition numbers for the eigenvectors of a real symmetric or complex Hermitian matrix or for the left or right singular vectors of a general m-by-n matrix. The reciprocal condition number is the 'gap' between the corresponding eigenvalue or singular value and the nearest other one.

The bound on the error, measured by angle in radians, in the I-th computed vector is given by

DLAMCH( 'E' ) * ( ANORM / SEP( I ) )

where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of the error bound.

DDISNA may also be used to compute error bounds for eigenvectors of the generalized symmetric definite eigenproblem.

So like, it's the "how much do I trust this vector" calculator 💅 the reciprocal condition number is literally the distance to the nearest neighbour, and if two eigenvalues are basically sitting on top of each other, that gap is tiny and the error bound blows up. Iykyk.

📋 Parameters

JOB — JOB is CHARACTER*1. Specifies for which problem the reciprocal condition numbers should be computed:
- = 'E': the eigenvectors of a symmetric/Hermitian matrix;
- = 'L': the left singular vectors of a general matrix;
- = 'R': the right singular vectors of a general matrix.
M — M is INTEGER. The number of rows of the matrix. M >= 0.
N — N is INTEGER. If JOB = 'L' or 'R', the number of columns of the matrix, in which case N >= 0. Ignored if JOB = 'E'.
D — D is DOUBLE PRECISION array, dimension (M) if JOB = 'E', dimension (min(M,N)) if JOB = 'L' or 'R'. The eigenvalues (if JOB = 'E') or singular values (if JOB = 'L' or 'R') of the matrix, in either increasing or decreasing order. If singular values, they must be non-negative.
SEP — SEP is DOUBLE PRECISION array, dimension (M) if JOB = 'E', dimension (min(M,N)) if JOB = 'L' or 'R'. The reciprocal condition numbers of the vectors.
INFO — INFO is INTEGER:
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.

Note that D must be sorted — either increasing or decreasing, pick a lane 🚦 — and if you're feeding it singular values they cannot be negative. Fr fr, it will not be cute about it.

✍️ Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine dlaed0 (integer ICOMPQ, integer QSIZ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldq, * ) Q, integer LDQ, double precision, dimension( ldqs, * ) QSTORE, integer LDQS, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding
eigenvectors of an unreduced symmetric tridiagonal matrix using the
divide and conquer method.

🎯 Purpose:

DLAED0 computes all eigenvalues and corresponding eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method.

Divide and conquer, bestie — she splits the problem up like it's group project drama and then reassembles it. Main character energy 🎀

📋 Parameters

ICOMPQ — ICOMPQ is INTEGER:
- = 0: Compute eigenvalues only.
- = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
- = 2: Compute eigenvalues and eigenvectors of tridiagonal matrix.
QSIZ — QSIZ is INTEGER. The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the main diagonal of the tridiagonal matrix. On exit, its eigenvalues.
E — E is DOUBLE PRECISION array, dimension (N-1). The off-diagonal elements of the tridiagonal matrix. On exit, E has been destroyed.

R.I.P. E 🪦 she does not survive the function. If you needed those off-diagonals later, bestie, save a copy first.

Q — Q is DOUBLE PRECISION array, dimension (LDQ, N). On entry, Q must contain an N-by-N orthogonal matrix.
- If ICOMPQ = 0, Q is not referenced.
- If ICOMPQ = 1, on entry, Q is a subset of the columns of the orthogonal matrix used to reduce the full matrix to tridiagonal form corresponding to the subset of the full matrix which is being decomposed at this time.
- If ICOMPQ = 2, on entry, Q will be the identity matrix. On exit, Q contains the eigenvectors of the tridiagonal matrix.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. If eigenvectors are desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
QSTORE — QSTORE is DOUBLE PRECISION array, dimension (LDQS, N). Referenced only when ICOMPQ = 1. Used to store parts of the eigenvector matrix when the updating matrix multiplies take place.
LDQS — LDQS is INTEGER. The leading dimension of the array QSTORE. If ICOMPQ = 1, then LDQS >= max(1,N). In any case, LDQS >= 1.
WORK — WORK is DOUBLE PRECISION array.
- If ICOMPQ = 0 or 1, the dimension of WORK must be at least 1 + 3*N + 2*N*lg N + 3*N**2 ( lg( N ) = smallest integer k such that 2^k >= N ).
- If ICOMPQ = 2, the dimension of WORK must be at least 4*N + N**2.
IWORK — IWORK is INTEGER array.
- If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 6 + 6*N + 5*N*lg N ( lg( N ) = smallest integer k such that 2^k >= N ).
- If ICOMPQ = 2, the dimension of IWORK must be at least 3 + 5*N.

Those workspace sizes are not a suggestion 🌸 they are the whole entire law. Undersize WORK or IWORK and it's giving undefined behaviour, and nobody wants that on their conscience.

INFO — INFO is INTEGER:
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1).

✍️ Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

subroutine dlaed1 (integer N, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal
matrix after modification by a rank-one symmetric matrix. Used when the
original matrix is tridiagonal.

🎯 Purpose:

DLAED1 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles the case in which eigenvalues only or eigenvalues and eigenvectors of a full symmetric matrix (which was reduced to tridiagonal form) are desired.

  T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

   where Z = Q**T*u, u is a vector of length N with ones in the
   CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine DLAED2. 💅

The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine DLAED4 (as called by DLAED3). This routine also calculates the eigenvectors of the current problem. ✨

The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem. 👑

Parameters 💖

N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix.
Q — Q is DOUBLE PRECISION array, dimension (LDQ,N). On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INDXQ is INTEGER array, dimension (N). On entry, the permutation which separately sorts the two subproblems in D into ascending order. On exit, the permutation which will reintegrate the subproblems back into sorted order, i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.
RHO — RHO is DOUBLE PRECISION. The subdiagonal entry used to create the rank-1 modification.
CUTPNT — CUTPNT is INTEGER. The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N/2.
WORK — WORK is DOUBLE PRECISION array, dimension (4*N + N**2).
IWORK — IWORK is INTEGER array, dimension (4*N).
INFO — INFO is INTEGER. = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not converge.

Author 🌟

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

Contributors: 🫶

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

✨ subroutine dlaed2 ✨

subroutine dlaed2 (integer K, integer N, integer N1, double precision,
    dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ,
    integer, dimension( * ) INDXQ, double precision RHO, double precision,
    dimension( * ) Z, double precision, dimension( * ) DLAMDA, double
    precision, dimension( * ) W, double precision, dimension( * ) Q2,
    integer, dimension( * ) INDX, integer, dimension( * ) INDXC, integer,
    dimension( * ) INDXP, integer, dimension( * ) COLTYP, integer INFO)

DLAED2 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal. 🔥

Purpose: 🎯

DLAED2 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny entry in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one. It's giving efficiency, bestie. 💅

Parameters 📋

K — K is INTEGER. The number of non-deflated eigenvalues, and the order of the related secular equation. 0 <= K <=N.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
N1 — N1 is INTEGER. The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= N1 <= N/2.
D — D is DOUBLE PRECISION array, dimension (N). On entry, D contains the eigenvalues of the two submatrices to be combined. On exit, D contains the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.
Q — Q is DOUBLE PRECISION array, dimension (LDQ, N). On entry, Q contains the eigenvectors of two submatrices in the two square blocks with corners at (1,1), (N1,N1) and (N1+1, N1+1), (N,N). On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INDXQ is INTEGER array, dimension (N). The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have N1 added to their values. Destroyed on exit.
RHO — RHO is DOUBLE PRECISION. On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by DLAED3.
Z — Z is DOUBLE PRECISION array, dimension (N). On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z have been destroyed by the updating process.
DLAMDA — DLAMDA is DOUBLE PRECISION array, dimension (N). A copy of the first K eigenvalues which will be used by DLAED3 to form the secular equation.
W — W is DOUBLE PRECISION array, dimension (N). The first k values of the final deflation-altered z-vector which will be passed to DLAED3.
Q2 — Q2 is DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2). A copy of the first K eigenvectors which will be used by DLAED3 in a matrix multiply (DGEMM) to solve for the new eigenvectors.
INDX — INDX is INTEGER array, dimension (N). The permutation used to sort the contents of DLAMDA into ascending order.
INDXC — INDXC is INTEGER array, dimension (N). The permutation used to arrange the columns of the deflated Q matrix into three groups: the first group contains non-zero elements only at and above N1, the second contains non-zero elements only below N1, and the third is dense.
INDXP — INDXP is INTEGER array, dimension (N). The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.
COLTYP — COLTYP is INTEGER array, dimension (N). During execution, a label which will indicate which of the following types a column in the Q2 matrix is: 1 : non-zero in the upper half only; 2 : dense; 3 : non-zero in the lower half only; 4 : deflated. On exit, COLTYP(i) is the number of columns of type i, for i=1 to 4 only.
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.

Contributors: 🫶

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

✨ subroutine dlaed3 ✨

subroutine dlaed3 (integer K, integer N, integer N1, double precision,
    dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ,
    double precision RHO, double precision, dimension( * ) DLAMDA, double
    precision, dimension( * ) Q2, integer, dimension( * ) INDX, integer,
    dimension( * ) CTOT, double precision, dimension( * ) W, double
    precision, dimension( * ) S, integer INFO)

DLAED3 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal. Iconic. 💅

Purpose: 🎯

DLAED3 finds the roots of the secular equation, as defined by the values in D, W, and RHO, between 1 and K. It makes the appropriate calls to DLAED4 and then updates the eigenvectors by multiplying the matrix of eigenvectors of the pair of eigensystems being combined by the matrix of eigenvectors of the K-by-K system which is solved here. ✨

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

Parameters 📋

K — K is INTEGER. The number of terms in the rational function to be solved by DLAED4. K >= 0.
N — N is INTEGER. The number of rows and columns in the Q matrix. N >= K (deflation may result in N>K).
N1 — N1 is INTEGER. The location of the last eigenvalue in the leading submatrix. min(1,N) <= N1 <= N/2.
D — D is DOUBLE PRECISION array, dimension (N). D(I) contains the updated eigenvalues for 1 <= I <= K.
Q — Q is DOUBLE PRECISION array, dimension (LDQ,N). Initially the first K columns are used as workspace. On output the columns 1 to K contain the updated eigenvectors.
LDQ
LDQ — INTEGER. The leading dimension of the array Q. Gotta be LDQ >= max(1,N), no exceptions bestie, the array will simply not cope otherwise. 💅
RHO — DOUBLE PRECISION. It's the value of the parameter in the rank one update equation. RHO >= 0 is required. Like, required required. 🫶
DLAMDA — DOUBLE PRECISION array, dimension (K). The first K elements hold the old roots of the deflated updating problem — these are the poles of the secular equation, iconic honestly. Heads up though: it may be changed on output by having its lowest order bit set to zero on Cray X-MP, Cray Y-MP, Cray-2, or Cray C-90, as described above. Yes that's real, yes that's cursed, we move. 😭
Q2 — DOUBLE PRECISION array, dimension (LDQ2*N). The first K columns of this matrix contain the non-deflated eigenvectors for the split problem. Slay. ✨
INDX — INTEGER array, dimension (N). The permutation used to arrange the columns of the deflated Q matrix into three groups (see DLAED2). Important detail: the rows of the eigenvectors found by DLAED4 must be likewise permuted before the matrix multiply can take place. Don't skip that step, it's giving chaos otherwise. 🙈
CTOT — INTEGER array, dimension (4). A count of the total number of the various types of columns in Q, as described in INDX. And the fourth column type? Any column which has been deflated. Gone. Deflated. 💀
W — DOUBLE PRECISION array, dimension (K). The first K elements contain the components of the deflation-adjusted updating vector. Destroyed on output, so don't get attached, she's not coming back. 😭
S — DOUBLE PRECISION array, dimension (N1 + 1)*K. Will contain the eigenvectors of the repaired matrix, which will be multiplied by the previously accumulated eigenvectors to update the system. Repaired matrix, we love a glow-up. 💖
INFO — INTEGER.
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value. Rude of you tbh.
- > 0: if INFO = 1, an eigenvalue did not converge. Not the vibe. 💀

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Modified by Francoise Tisseur, University of Tennessee

🌟 subroutine `dlaed4` — the single-root finder ✨

subroutine dlaed4 (integer N, integer I, double precision, dimension( * ) D,
    double precision, dimension( * ) Z, double precision, dimension( * )
    DELTA, double precision RHO, double precision DLAM, integer INFO)

DLAED4 used by DSTEDC. Finds a single root of the secular equation. One root. Singular. Main character energy, honestly. 🎀

🎯 Purpose:

This subroutine computes the I-th updated eigenvalue of a symmetric rank-one modification to a diagonal matrix whose elements are given in the array d, and that

           D(i) < D(j)  for  i < j

and that RHO > 0. This is arranged by the calling routine, and is no loss in generality. The rank-one modified system is thus

           diag( D )  +  RHO * Z * Z_transpose.

where we assume the Euclidean norm of Z is 1. ✨

The method consists of approximating the rational functions in the secular equation by simpler interpolating rational functions. Fancy? Yes. Necessary? Also yes bestie. 💅

📋 Parameters

N — INTEGER. The length of all arrays. Every single one. 🫶
I — INTEGER. The index of the eigenvalue to be computed. 1 <= I <= N.
D — DOUBLE PRECISION array, dimension (N). The original eigenvalues. It is assumed that they are in order, D(I) < D(J) for I < J. Sorted and unbothered. 👑
Z — DOUBLE PRECISION array, dimension (N). The components of the updating vector.
DELTA — DOUBLE PRECISION array, dimension (N). If N > 2, DELTA contains (D(j) - lambda_I) in its j-th component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5 for detail. The vector DELTA contains the information necessary to construct the eigenvectors by DLAED3 and DLAED9. It's giving lore. 📜
RHO — DOUBLE PRECISION. The scalar in the symmetric updating formula.
DLAM — DOUBLE PRECISION. The computed lambda_I, the I-th updated eigenvalue. ✨
INFO — INTEGER.
- = 0: successful exit.
- > 0: if INFO = 1, the updating process failed. 😭

🔍 Internal Parameters:

Logical variable ORGATI (origin-at-i?) is used for distinguishing whether D(i) or D(i+1) is treated as the origin.

          ORGATI = .true.    origin at i
          ORGATI = .false.   origin at i+1

Logical variable SWTCH3 (switch-for-3-poles?) is for noting if we are working with THREE poles! Three! A whole trio. 🎀

MAXIT is the maximum number of iterations allowed for each eigenvalue.

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

🌈 subroutine `dlaed5` — the 2-by-2 girlie 💖

subroutine dlaed5 (integer I, double precision, dimension( 2 ) D, double
    precision, dimension( 2 ) Z, double precision, dimension( 2 ) DELTA,
    double precision RHO, double precision DLAM)

DLAED5 used by DSTEDC. Solves the 2-by-2 secular equation. Tiny. Adorable. Fully solved. ✨

🎯 Purpose:

This subroutine computes the I-th eigenvalue of a symmetric rank-one modification of a 2-by-2 diagonal matrix

           diag( D )  +  RHO * Z * transpose(Z) .

The diagonal elements in the array D are assumed to satisfy

           D(i) < D(j)  for  i < j .

We also assume RHO > 0 and that the Euclidean norm of the vector Z is one. No cap. 🫶

📋 Parameters

I — INTEGER. The index of the eigenvalue to be computed. I = 1 or I = 2. That's it. That's the whole range. 🎀
D — DOUBLE PRECISION array, dimension (2). The original eigenvalues. We assume D(1) < D(2).
Z — DOUBLE PRECISION array, dimension (2). The components of the updating vector.
DELTA — DOUBLE PRECISION array, dimension (2). The vector DELTA contains the information necessary to construct the eigenvectors.
RHO — DOUBLE PRECISION. The scalar in the symmetric updating formula.
DLAM — DOUBLE PRECISION. The computed lambda_I, the I-th updated eigenvalue. 💅

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

🔥 subroutine `dlaed6` — one Newton step, that's the whole job 💀

subroutine dlaed6 (integer KNITER, logical ORGATI, double precision RHO,
    double precision, dimension( 3 ) D, double precision, dimension( 3 ) Z,
    double precision FINIT, double precision TAU, integer INFO)

DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation. One step. Iconic minimalism. ✨

🎯 Purpose:

DLAED6 computes the positive or negative root (closest to the origin) of

                         z(1)        z(2)        z(3)
        f(x) =   rho + --------- + ---------- + ---------
                        d(1)-x      d(2)-x      d(3)-x

It is assumed that

              if ORGATI = .true. the root is between d(2) and d(3);
              otherwise it is between d(1) and d(2)

This routine will be called by DLAED4 when necessary. In most cases, the root sought is the smallest in magnitude, though it might not be in some extremely rare situations. So like, usually the smallest. Usually. 👀

📋 Parameters

KNITER — INTEGER. Refer to DLAED4 for its significance.
ORGATI — LOGICAL. If ORGATI is true, the needed root is between d(2) and d(3); otherwise it is between d(1) and d(2). See DLAED4 for further details.
RHO — DOUBLE PRECISION. Refer to the equation f(x) above. ⬆️
D — DOUBLE PRECISION array, dimension (3). D satisfies d(1) < d(2) < d(3). Perfectly ordered, we stan. 💅
Z — DOUBLE PRECISION array, dimension (3). Each of the elements in z must be positive. All of them. Positivity check passed bestie. 🫶
FINIT — DOUBLE PRECISION. The value of f at 0. It is more accurate than the one evaluated inside this routine (if someone wants to do so). So, like, use it, it's the good one. ✨
TAU — DOUBLE PRECISION. The root of the equation f(x).
INFO — INTEGER.
- = 0: successful exit
- > 0: if INFO = 1, failure to converge 😭

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

📜 Further Details:

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

05/10/06: Modified from a new version of Ren-Cang Li, use Gragg-Thornton-Warner cubic convergent scheme for better stability. Stability upgrade, we love to see it. 🌈

🫶 Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

👑 subroutine `dlaed7` — the updated eigensystem girlie 🦋

subroutine dlaed7 (integer ICOMPQ, integer N, integer QSIZ, integer TLVLS,
    integer CURLVL, integer CURPBM, double precision, dimension( * ) D,
    double precision, dimension( ldq, * ) Q, integer LDQ, integer,
    dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double
    precision, dimension( * ) QSTORE, integer, dimension( * ) QPTR, integer,
    dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension(
    * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, double precision,
    dimension( 2, * ) GIVNUM, double precision, dimension( * ) WORK,
    integer, dimension( * ) IWORK, integer INFO)

DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense. Dense matrix, dense lore, we're seated. 💅

🎯 Purpose:

DLAED7 is literally the updated-eigensystem glow-up for a diagonal matrix after a rank-one symmetric matrix has done its little makeover. ✨ It's only for the eigenproblem that needs all eigenvalues and optionally eigenvectors of a dense symmetric matrix that's already been reduced to tridiagonal form. If you want all eigenvalues and eigenvectors of a symmetric tridiagonal matrix, that's DLAED1's job, not this one. 💅

T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**Tu, u is a vector of length N with ones in the
CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm is a three-stage saga:

The first stage deflates the size of the problem when there are multiple eigenvalues or if there's a zero in the Z vector. For each such occurrence, the dimension of the secular equation problem is reduced by one. This stage is performed by the routine DLAED8.
The second stage calculates the updated eigenvalues. This is done by finding the roots of the secular equation via the routine DLAED4 (as called by DLAED9). This routine also calculates the eigenvectors of the current problem.
The final stage computes the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.

🎀 Parameters 🎀

ICOMPQ — ICOMPQ is INTEGER. = 0: compute eigenvalues only. = 1: compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
QSIZ — QSIZ is INTEGER. The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
TLVLS — TLVLS is INTEGER. The total number of merging levels in the overall divide and conquer tree.
CURLVL — CURLVL is INTEGER. The current level in the overall merge routine, 0 <= CURLVL <= TLVLS.
CURPBM — CURPBM is INTEGER. The current problem in the current level in the overall merge routine (counting from upper left to lower right).
D — D is DOUBLE PRECISION array, dimension (N). On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix.
Q — Q is DOUBLE PRECISION array, dimension (LDQ, N). On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INDXQ is INTEGER array, dimension (N). The permutation which will reintegrate the subproblem just solved back into sorted order, i.e., D( INDXQ( I = 1, N ) ) will be in ascending order.
RHO — RHO is DOUBLE PRECISION. The subdiagonal element used to create the rank-1 modification.
CUTPNT — CUTPNT is INTEGER. Contains the location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.
QSTORE — QSTORE is DOUBLE PRECISION array, dimension (N**2+1). Stores eigenvectors of submatrices encountered during divide and conquer, packed together. QPTR points to beginning of the submatrices.
QPTR — QPTR is INTEGER array, dimension (N+2). List of indices pointing to beginning of submatrices stored in QSTORE. The submatrices are numbered starting at the bottom left of the divide and conquer tree, from left to right and bottom to top.
PRMPTR — PRMPTR is INTEGER array, dimension (N lg N). Contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and also the size of the full, non-deflated problem.
PERM — PERM is INTEGER array, dimension (N lg N). Contains the permutations (from deflation and sorting) to be applied to each eigenblock.
GIVPTR — GIVPTR is INTEGER array, dimension (N lg N). Contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations.
GIVCOL — GIVCOL is INTEGER array, dimension (2, N lg N). Each pair of numbers indicates a pair of columns to take place in a Givens rotation.
GIVNUM — GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N). Each number indicates the S value to be used in the corresponding Givens rotation.
WORK — WORK is DOUBLE PRECISION array, dimension (3*N+2*QSIZ*N).
IWORK — IWORK is INTEGER array, dimension (4*N).
INFO — INFO is INTEGER. = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not converge.

👑 Author 👑

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

💖 Contributors 💖

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

🧚 subroutine dlaed8 🧚

   subroutine dlaed8 (integer ICOMPQ, integer K, integer N, integer QSIZ,
       double precision, dimension( * ) D, double precision, dimension( ldq, *
       ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO,
       integer CUTPNT, double precision, dimension( * ) Z, double precision,
       dimension( * ) DLAMDA, double precision, dimension( ldq2, * ) Q2,
       integer LDQ2, double precision, dimension( * ) W, integer, dimension( *
       ) PERM, integer GIVPTR, integer, dimension( 2, * ) GIVCOL, double
       precision, dimension( 2, * ) GIVNUM, integer, dimension( * ) INDXP,
       integer, dimension( * ) INDX, integer INFO)

DLAED8 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.

🎀 Purpose 🎀

DLAED8 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny element in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one.

💅 Parameters 💅

ICOMPQ — ICOMPQ is INTEGER. = 0: compute eigenvalues only. = 1: compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
K — K is INTEGER. The number of non-deflated eigenvalues, and the order of the related secular equation.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
QSIZ — QSIZ is INTEGER. The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the eigenvalues of the two submatrices to be combined. On exit, the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.
Q — Q is DOUBLE PRECISION array, dimension (LDQ,N). If ICOMPQ = 0, Q is not referenced. Otherwise, on entry, Q contains the eigenvectors of the partially solved system which has been previously updated in matrix multiplies with other partially solved eigensystems. On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INDXQ is INTEGER array, dimension (N). The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have CUTPNT added to their values in order to be accurate.
RHO — RHO is DOUBLE PRECISION. On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by DLAED3.
CUTPNT — CUTPNT is INTEGER. The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.
Z — Z is DOUBLE PRECISION array, dimension (N). On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z are destroyed by the updating process.
DLAMDA — DLAMDA is DOUBLE PRECISION array, dimension (N). A copy of the first K eigenvalues which will be used by DLAED3 to form the secular equation.
Q2 —

`Q2` 🫶

Q2 — a DOUBLE PRECISION array, dimension (LDQ2,N). Okay so if ICOMPQ = 0, bestie, Q2 is not referenced at all, it's just vibing. Otherwise it holds a copy of the first K eigenvectors, which DLAED7 is gonna use in a matrix multiply (DGEMM) to update the new eigenvectors. ✨
LDQ2 — LDQ2 is INTEGER. It's the leading dimension of the array Q2. And LDQ2 >= max(1,N), no exceptions. 💅
W — W is a DOUBLE PRECISION array, dimension (N). The first k values of the final deflation-altered z-vector live here, and they get passed along to DLAED3. 🎀
PERM — PERM is an INTEGER array, dimension (N). It's the permutations (from deflation and sorting, the whole glow-up) that get applied to each eigenblock.
GIVPTR — GIVPTR is INTEGER. It's literally just how many Givens rotations took place in this subproblem. Counting queen. 👑
GIVCOL — GIVCOL is an INTEGER array, dimension (2, N). Each pair of numbers indicates a pair of columns to take place in a Givens rotation.
GIVNUM — GIVNUM is a DOUBLE PRECISION array, dimension (2, N). Each number here indicates the S value used in the corresponding Givens rotation.
INDXP — INDXP is an INTEGER array, dimension (N). It's the permutation used to place deflated values of D at the end of the array. So INDXP(1:K) points to the nondeflated D-values, and INDXP(K+1:N) points to the deflated eigenvalues. Two little friend groups. 🫂
INDX — INDX is an INTEGER array, dimension (N). The permutation used to sort the contents of D into ascending order. Organization is self-care, fr fr. 🌸
INFO — INFO is INTEGER, and it's giving report card:
- = 0: successful exit, we ate. ✅
- < 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.

Contributors 🌈

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

🦋 subroutine `dlaed9` 🦋

subroutine dlaed9 (integer K, integer KSTART, integer KSTOP, integer N,
    double precision, dimension( * ) D, double precision, dimension( ldq, *
    ) Q, integer LDQ, double precision RHO, double precision, dimension( * )
    DLAMDA, double precision, dimension( * ) W, double precision, dimension(
    lds, * ) S, integer LDS, integer INFO)

DLAED9 is used by DSTEDC. It finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense. Periodt. 💅

Purpose 🎯

DLAED9 finds the roots of the secular equation, as defined by the values in D, Z, and RHO, between KSTART and KSTOP. It makes the appropriate calls to DLAED4 and then stores the new matrix of eigenvectors for use in calculating the next level of Z vectors. Very main character of it. ✨

Parameters 🌸

K — K is INTEGER. The number of terms in the rational function to be solved by DLAED4. And K >= 0, bestie.
KSTART — KSTART is INTEGER. That's it, that's the whole description, no notes. 🎀
KSTOP — KSTOP is INTEGER. The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP, are to be computed. And 1 <= KSTART <= KSTOP <= K. The range has to actually make sense, okay?
N — N is INTEGER. The number of rows and columns in the Q matrix. N >= K (delation may result in N > K). Not me side-eyeing that spelling, moving on. 💀
D — D is a DOUBLE PRECISION array, dimension (N). D(I) contains the updated eigenvalues for KSTART <= I <= KSTOP.
Q — Q is a DOUBLE PRECISION array, dimension (LDQ,N).
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= max( 1, N ). 🫶
RHO — RHO is DOUBLE PRECISION. The value of the parameter in the rank one update equation. RHO >= 0 required — like, actually required, don't play. 🚨
DLAMDA — DLAMDA is a DOUBLE PRECISION array, dimension (K). The first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation. Iconic poles. 🧚
W — W is a DOUBLE PRECISION array, dimension (K). The first K elements of this array contain the components of the deflation-adjusted updating vector.
S — S is a DOUBLE PRECISION array, dimension (LDS, K). It will contain the eigenvectors of the repaired matrix, which will be stored for subsequent Z vector calculation and multiplied by the previously accumulated eigenvectors to update the system. Busy little array, honestly. 💖
LDS — LDS is INTEGER. The leading dimension of S. LDS >= max( 1, K ).
INFO — INFO is INTEGER:
- = 0: successful exit. ✅
- < 0: if INFO = -i, the i-th argument had an illegal value. 🙈
- > 0: if INFO = 1, an eigenvalue did not converge. 😭

Author 💜

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors 🌈

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

🧚 subroutine `dlaeda` 🧚

subroutine dlaeda (integer N, integer TLVLS, integer CURLVL, integer CURPBM,
    integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer,
    dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, double
    precision, dimension( 2, * ) GIVNUM, double precision, dimension( * ) Q,
    integer, dimension( * ) QPTR, double precision, dimension( * ) Z, double
    precision, dimension( * ) ZTEMP, integer INFO)

DLAEDA is used by DSTEDC. It computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense. We love a dense-matrix moment. 🔥

Purpose 🎯

DLAEDA computes the Z vector corresponding to the merge step in the CURLVLth step of the merge process with TLVLS steps for the CURPBMth problem. It's giving very specific. And correct. ✨

Parameters 🌸

N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
TLVLS — TLVLS is INTEGER. The total number of merging levels in the overall divide and conquer tree. 🌳
CURLVL — CURLVL is INTEGER. The current level in the overall merge routine, 0 <= curlvl <= tlvls.
CURPBM — CURPBM is INTEGER. The current problem in the current level in the overall merge routine (counting from upper left to lower right). Very organized, we stan. 💅
PRMPTR — PRMPTR is an INTEGER array, dimension (N lg N). It contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and incidentally the size of the full, non-deflated problem.
PERM — PERM is an INTEGER array, dimension (N lg N). It contains the permutations (from deflation and sorting) to be applied to each eigenblock.
GIVPTR — GIVPTR is an INTEGER array, dimension (N lg N). It contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations. Counting again, love that for us. 👑
GIVCOL — GIVCOL is an INTEGER array, dimension (2, N lg N). Each pair of numbers indicates a pair of columns to take place in a Givens rotation.
GIVNUM — GIVNUM is a DOUBLE PRECISION array, dimension (2, N lg N). Each number indicates the S value to be used in the corresponding Givens rotation.
Q — Q is a DOUBLE PRECISION array, dimension (N**2). It contains the square eigenblocks from previous levels, and the starting positions for blocks are given by QPTR.
QPTR — QPTR is an INTEGER array, dimension (N+2). It contains a list of pointers which indicate where in Q an eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates the size of the block. Math girlies rise. 🌈
Z — Z is a DOUBLE PRECISION array, dimension (N). On output this vector contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix).
ZTEMP — ZTEMP is a DOUBLE PRECISION array, dimension (N). Just chilling there. 🎀
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.

Contributors 🌈

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

💖 subroutine `dlagtf` 💖

subroutine dlagtf (integer N, double precision, dimension( * ) A, double
    precision LAMBDA, double precision, dimension( * ) B, double precision,
    dimension( * ) C, double precision TOL, double precision, dimension( * )
    D, integer, dimension( * ) IN, integer INFO)

DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges. Glow-up via Gaussian elimination, basically. ✨

Purpose 🎯

DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n tridiagonal matrix and lambda is a scalar, as

   T - lambda*I = PLU,

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The factorization is obtained by Gaussian elimination with partial pivoting and implicit row scaling. Clean, efficient, no crumbs. 💅

The parameter LAMBDA is included in the routine so that DLAGTF may be used, in conjunction with DLAGTS, to obtain eigenvectors of T by inverse iteration. Bestie duo behavior. 🫶

Parameters 🌸

N — N is INTEGER. The order of the matrix T.
A — A is a DOUBLE PRECISION array, dimension (N). On entry, A must contain the diagonal elements of T.

On exit, A is overwritten by the n diagonal elements of the upper triangular matrix U of the factorization of T. She gets replaced, that's just her arc. 💖

LAMBDA — LAMBDA is DOUBLE PRECISION. On entry, the scalar lambda.
B —

still in the parameter list, bestie — we're deep in the exit values now 💖

B — B is DOUBLE PRECISION array, dimension (N-1). On entry, B must contain the (n-1) super-diagonal elements of T. On exit, B is overwritten by the (n-1) super-diagonal elements of the matrix U of the factorization of T. in goes T, out comes U bestie 💅
C — C is DOUBLE PRECISION array, dimension (N-1). On entry, C must contain the (n-1) sub-diagonal elements of T. On exit, C is overwritten by the (n-1) sub-diagonal elements of the matrix L of the factorization of T. same deal, different diagonal 🦋
TOL — TOL is DOUBLE PRECISION. On entry, a relative tolerance used to indicate whether or not the matrix (T - lambda*I) is nearly singular. TOL should normally be chosen as approximately the largest relative error in the elements of T. For example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5*10**(-4). If TOL is supplied as less than eps, where eps is the relative machine precision, then the value eps is used in place of TOL. so like, if you pass something smol, it gets bumped up to eps, no cap ⭐
D — D is DOUBLE PRECISION array, dimension (N-2). On exit, D is overwritten by the (n-2) second super-diagonal elements of the matrix U of the factorization of T. the second super-diagonal, the understudy of diagonals ✨
IN — IN is INTEGER array, dimension (N). On exit, IN contains details of the permutation matrix P. If an interchange occurred at the kth step of the elimination, then IN(k) = 1, otherwise IN(k) = 0. The element IN(n) returns the smallest positive integer j such that

   abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,

where norm( A(j) ) denotes the sum of the absolute values of the jth row of the matrix A. If no such j exists then IN(n) is returned as zero. If IN(n) is returned as positive, then a diagonal element of U is small, indicating that (T - lambda*I) is singular or nearly singular, main character IN(n) literally snitching on your singularity 🎀

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

Author

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

🌸 subroutine dlamrg 🌸

subroutine dlamrg (integer N1, integer N2, double precision, dimension( * )
    A, integer DTRD1, integer DTRD2, integer, dimension( * ) INDEX)

DLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order. a little merge-matchmaker, she's setting up sorted sets fr fr 💜

Purpose:

DLAMRG will create a permutation list which will merge the elements of A (which is composed of two independently sorted sets) into a single set which is sorted in ascending order.

Parameters

N1 — N1 is INTEGER
N2 — N2 is INTEGER. These arguments contain the respective lengths of the two sorted lists to be merged. two lists, two lengths, one destiny 🫶
A — A is DOUBLE PRECISION array, dimension (N1+N2). The first N1 elements of A contain a list of numbers which are sorted in either ascending or descending order. Likewise for the final N2 elements.
DTRD1 — DTRD1 is INTEGER
DTRD2 — DTRD2 is INTEGER. These are the strides to be taken through the array A. Allowable strides are 1 and -1. They indicate whether a subset of A is sorted in ascending (DTRDx = 1) or descending (DTRDx = -1) order. up or down, pick your vibe check 💅
INDEX — INDEX is INTEGER array, dimension (N1+N2). On exit this array will contain a permutation such that if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be sorted in ascending order. the glow-up pipeline: A → INDEX → B, sorted and iconic ✨

Author

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

🦋 subroutine dlartgs 🦋

subroutine dlartgs (double precision X, double precision Y, double precision
    SIGMA, double precision CS, double precision SN)

DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem. yes bestie, we're creating a bulge on purpose, it's giving chaos but make it numerical 🔥

Purpose:

DLARTGS generates a plane rotation designed to introduce a bulge in Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD problem. X and Y are the top-row entries, and SIGMA is the shift. The computed CS and SN define a plane rotation satisfying

   [  CS  SN  ]  .  [ X^2 - SIGMA ]  =  [ R ],
   [ -SN  CS  ]     [    X * Y    ]     [ 0 ]

with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the rotation is by PI/2. PI/2 rotation, she said "I'll just turn all the way, whatever" 💀

Parameters

X — X is DOUBLE PRECISION. The (1,1) entry of an upper bidiagonal matrix.
Y — Y is DOUBLE PRECISION. The (1,2) entry of an upper bidiagonal matrix.
SIGMA — SIGMA is DOUBLE PRECISION. The shift. the shift, the moment, the plot twist ⭐
CS — CS is DOUBLE PRECISION. The cosine of the rotation.
SN — SN is DOUBLE PRECISION. The sine of the rotation. cosine and sine, the duo that ate 🎀

Author

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

💖 subroutine dlasq1 💖

subroutine dlasq1 (integer N, double precision, dimension( * ) D, double
    precision, dimension( * ) E, double precision, dimension( * ) WORK,
    integer INFO)

DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr. she's the backstage icon behind sbdsqr, no cap 👑

Purpose:

DLASQ1 computes the singular values of a real N-by-N bidiagonal matrix with diagonal D and off-diagonal E. The singular values are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. The algorithm was first presented in

'Accurate singular values and differential qd algorithms' by K. V. Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, 1994,

and the present implementation is described in 'An implementation of the dqds Algorithm (Positive Case)', LAPACK Working Note. read the lore, bestie, it's literally cited 📚

Parameters

N — N is INTEGER. The number of rows and columns in the matrix. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, D contains the diagonal elements of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in decreasing order. big to smol, descending, as she deserves 💅
E — E is DOUBLE PRECISION array, dimension (N). On entry, elements E(1:N-1) contain the off-diagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E is overwritten. poof, gone, replaced, iconic ✨
WORK — WORK is DOUBLE PRECISION array, dimension (4*N) scratch space, she needs her room 🌸
INFO — INFO is INTEGER
- = 0: successful exit 👑
- < 0: if INFO = -i, the i-th argument had an illegal value 💀
- > 0: the algorithm failed
  - = 1, a split was marked by a positive value in E
  - = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit D and E represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into DLASQ1 ngl that's kind of a serve, you can just run it back 🔥
  - = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)

Author

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

⭐ subroutine dlasq2 ⭐

subroutine dlasq2 (integer N, double precision, dimension( * ) Z, integer
    INFO)

DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr. two callers, one legend, we stan a reused queen 🎀

Purpose:

DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd array Z to high relative accuracy are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. yes it says it twice, we're keeping it, that's the original's whole vibe fr fr 💜

To see the relation of Z to the tridiagonal matrix, let L be a unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and let U be an upper bidiagonal matrix with 1's above and diagonal Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the symmetric tridiagonal to which it is similar. even indices go down, odd indices go up, that's the split, iykyk 🦋

Note : DLASQ2 defines a logical variable, IEEE, which is true on machines which follow ieee-754 floating-point standard in their handling of infinities and NaNs, and false otherwise. This variable is passed to DLASQ3. a whole bool just for the IEEE girlies ⭐

Parameters

N — N is INTEGER. The number of rows and columns in the matrix. N >= 0.
Z — Z is DOUBLE PRECISION array, dimension ( 4*N ). On entry Z holds the qd array. On exit, entries 1 to N hold the eigenvalues in decreasing order, Z( 2*N+1 ) holds the trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of shifts that failed. she's a packed array bestie, everything has an address, the tea is all in there 🍓
INFO — INFO is INTEGER
- = 0: successful exit 👑
- < 0: if the i-th argument is a scalar and had an illegal value, then INFO = -i, if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j) yes the array one is a whole arithmetic situation, that's just how it is 💅
- > 0: the algorithm failed
  - = 1, a split was marked by a positive value in E
  - = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit Z holds a qd array with the same eigenvalues as the given Z. not a loss, a handoff, the eigenvalues are still in there 🔥
  - = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)

Author

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

Further Details:

🌸 okay so before we even get into dlasq3 proper, the source keeps a lil' housekeeping note in the header — the internal receipts, if you will. Here's what they say, verbatim-ish but in my voice:

I0:N0 — defines a current unreduced segment of Z. 💅
The shifts are accumulated in SIGMA.
Iteration count lives in ITER.
Ping-pong is controlled by PP (it alternates between 0 and 1).

🦋 `subroutine dlasq3` 🦋

subroutine dlasq3 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, double precision DMIN, double precision SIGMA, double precision DESIG, double precision QMAX, integer NFAIL, integer ITER, integer NDIV, logical IEEE, integer TTYPE, double precision DMIN1, double precision DMIN2, double precision DN, double precision DN1, double precision DN2, double precision G, double precision TAU)

DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr. 🎀

Purpose:

DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. In case of failure it changes shifts, and tries again until output is positive.

Parameters

I0 — I0 is INTEGER. First index. ✨
N0 — N0 is INTEGER. Last index. ✨
Z — Z is DOUBLE PRECISION array, dimension ( 4*N0 ). Z holds the qd array. 🫶
PP — PP is INTEGER. PP=0 for ping, PP=1 for pong. PP=2 indicates that flipping was applied to the Z array and that the initial tests for deflation should not be performed. 👀
DMIN — DMIN is DOUBLE PRECISION. Minimum value of d.
SIGMA — SIGMA is DOUBLE PRECISION. Sum of shifts used in current segment.
DESIG — DESIG is DOUBLE PRECISION. Lower order part of SIGMA.
QMAX — QMAX is DOUBLE PRECISION. Maximum value of q.
NFAIL — NFAIL is INTEGER. Increment NFAIL by 1 each time the shift was too big. 😭
ITER — ITER is INTEGER. Increment ITER by 1 for each iteration.
NDIV — NDIV is INTEGER. Increment NDIV by 1 for each division.
IEEE — IEEE is LOGICAL. Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
TTYPE — TTYPE is INTEGER. Shift type.
DMIN1 — DMIN1 is DOUBLE PRECISION.
DMIN2 — DMIN2 is DOUBLE PRECISION.
DN — DN is DOUBLE PRECISION.
DN1 — DN1 is DOUBLE PRECISION.
DN2 — DN2 is DOUBLE PRECISION.
G — G is DOUBLE PRECISION.
TAU — TAU is DOUBLE PRECISION.

💜 and the tea on that last lil' group (DMIN1, DMIN2, DN, DN1, DN2, G, TAU): these are passed as arguments in order to save their values between calls to DLASQ3. It's giving ✨memory✨ but make it arguments.

Author

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

🧚 `subroutine dlasq4` 🧚

subroutine dlasq4 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, integer N0IN, double precision DMIN, double precision DMIN1, double precision DMIN2, double precision DN, double precision DN1, double precision DN2, double precision TAU, integer TTYPE, double precision G)

DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. 🎀

Purpose:

DLASQ4 computes an approximation TAU to the smallest eigenvalue using values of d from the previous transform.

Parameters

I0 — I0 is INTEGER. First index. ✨
N0 — N0 is INTEGER. Last index. ✨
Z — Z is DOUBLE PRECISION array, dimension ( 4*N0 ). Z holds the qd array. 🫶
PP — PP is INTEGER. PP=0 for ping, PP=1 for pong.
N0IN — N0IN is INTEGER. The value of N0 at start of EIGTEST. 👑
DMIN — DMIN is DOUBLE PRECISION. Minimum value of d.
DMIN1 — DMIN1 is DOUBLE PRECISION. Minimum value of d, excluding D( N0 ).
DMIN2 — DMIN2 is DOUBLE PRECISION. Minimum value of d, excluding D( N0 ) and D( N0-1 ).
DN — DN is DOUBLE PRECISION. d(N)
DN1 — DN1 is DOUBLE PRECISION. d(N-1)
DN2 — DN2 is DOUBLE PRECISION. d(N-2)
TAU — TAU is DOUBLE PRECISION. This is the shift. 💅
TTYPE — TTYPE is INTEGER. Shift type.
G — G is DOUBLE PRECISION. G is passed as an argument in order to save its value between calls to DLASQ4.

Author

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

Further Details:

CNST1 = 9/16

⭐ that's it, that's the whole further details section. no notes. iconic.

🌈 `subroutine dlasq5` 🌈

subroutine dlasq5 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, double precision TAU, double precision SIGMA, double precision DMIN, double precision DMIN1, double precision DMIN2, double precision DN, double precision DNM1, double precision DNM2, logical IEEE, double precision EPS)

DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. 🎀

Purpose:

DLASQ5 computes one dqds transform in ping-pong form, one version for IEEE machines another for non IEEE machines. ✨

Parameters

I0 — I0 is INTEGER. First index. ✨
N0 — N0 is INTEGER. Last index. ✨
Z — Z is DOUBLE PRECISION array, dimension ( 4*N ). Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument. 👀
PP — PP is INTEGER. PP=0 for ping, PP=1 for pong.
TAU — TAU is DOUBLE PRECISION. This is the shift. 💅
SIGMA — SIGMA is DOUBLE PRECISION. This is the accumulated shift up to this step.
DMIN — DMIN is DOUBLE PRECISION. Minimum value of d.
DMIN1 — DMIN1 is DOUBLE PRECISION. Minimum value of d, excluding D( N0 ).
DMIN2 — DMIN2 is DOUBLE PRECISION. Minimum value of d, excluding D( N0 ) and D( N0-1 ).
DN — DN is DOUBLE PRECISION. d(N0), the last value of d.
DNM1 — DNM1 is DOUBLE PRECISION. d(N0-1).
DNM2 — DNM2 is DOUBLE PRECISION. d(N0-2).
IEEE — IEEE is LOGICAL. Flag for IEEE or non IEEE arithmetic.
EPS — EPS is DOUBLE PRECISION. This is the value of epsilon used. 🫶

Author

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

🍓 `subroutine dlasq6` 🍓

subroutine dlasq6 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, double precision DMIN, double precision DMIN1, double precision DMIN2, double precision DN, double precision DNM1, double precision DNM2)

DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. 🎀

Purpose:

DLASQ6 computes one dqd (shift equal to zero) transform in ping-pong form, with protection against underflow and overflow. 🛡️

Parameters

I0 — I0 is INTEGER. First index. ✨
N0 — N0 is INTEGER. Last index. ✨
Z — Z is DOUBLE PRECISION array, dimension ( 4*N ). Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument. 👀
PP — PP is INTEGER. PP=0 for ping, PP=1 for pong.
DMIN — DMIN is DOUBLE PRECISION. Minimum value of d.
DMIN1 — DMIN1 is DOUBLE PRECISION. Minimum value of d, excluding D( N0 ).
DMIN2 — DMIN2 is DOUBLE PRECISION. Minimum value of d, excluding D( N0 ) and D( N0-1 ).
DN — DN is DOUBLE PRECISION. d(N0), the last value of d.
DNM1 — DNM1 is DOUBLE PRECISION. d(N0-1).
DNM2 — DNM2 is DOUBLE PRECISION. d(N0-2).

Author

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

👑 `subroutine dlasrt` 👑

subroutine dlasrt (character ID, integer N, double precision, dimension( * ) D, integer INFO)

DLASRT sorts numbers in increasing or decreasing order. 🎀

Purpose:

Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ).

Use Quick Sort, reverting to Insertion sort on arrays of size <= 20. Dimension of STACK limits N to about 2**32. 🔥

Parameters

ID — ID is CHARACTER*1.
- = 'I': sort D in increasing order;
- = 'D': sort D in decreasing order. 💅
N — N is INTEGER. The length of the array D.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the array to be sorted. On exit, D has been sorted into increasing order ( D(1) <= ... <= D(N) ) or into decreasing order ( D(1) >= ... >= D(N) ), depending on ID. ✨
INFO — 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 dstebz` 💜

subroutine dstebz (character RANGE, character ORDER, integer N, double precision VL, double precision VU, integer IL, integer IU, double precision ABSTOL, double precision, dimension( * ) D, double precision, dimension( * ) E, integer M, integer NSPLIT, double precision, dimension( * ) W, integer, dimension( * ) IBLOCK, integer, dimension( * ) ISPLIT, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

DSTEBZ 🎀

Purpose:

DSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix T. The user may ask for all eigenvalues, all eigenvalues in the half-open interval (VL, VU], or the IL-th through IU-th eigenvalues. ✨

okay so before you even touch the args, there's a whole vibe-check on scaling you cannot skip: to dodge overflow, the matrix has to be scaled so that its largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value — and for the absolute bestest accuracy, it shouldn't be much smaller than that either. ✨

And if you wanna read the receipts yourself, see W. Kahan, 'Accurate Eigenvalues of a Symmetric Tridiagonal Matrix', Report CS41, Computer Science Dept., Stanford University, July 21, 1966. 💅

🌸 Parameters 🌸

RANGE — RANGE is CHARACTER*1
- = 'A': ('All') — all eigenvalues will be found.
- = 'V': ('Value') — all eigenvalues in the half-open interval (VL, VU] will be found.
- = 'I': ('Index') — the IL-th through IU-th eigenvalues (of the entire matrix) will be found.
ORDER — ORDER is CHARACTER*1
- = 'B': ('By Block') — the eigenvalues will be grouped by split-off block (see IBLOCK, ISPLIT) and ordered from smallest to largest within the block.
- = 'E': ('Entire matrix') — the eigenvalues for the entire matrix will be ordered from smallest to largest.
N — N is INTEGER. The order of the tridiagonal matrix T. N >= 0. 👑
VL — VL is DOUBLE PRECISION. If RANGE='V', this is the lower bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.
VU — VU is DOUBLE PRECISION. If RANGE='V', this is the upper bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.
IL — IL is INTEGER. If RANGE='I', this is the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.
IU — IU is INTEGER. If RANGE='I', this is the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.
ABSTOL — ABSTOL is DOUBLE PRECISION. The absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered located if it's been determined to lie in an interval whose width is ABSTOL or less. If ABSTOL is less than or equal to zero, then ULP*|T| will be used, where |T| means the 1-norm of T.

bestie listen — eigenvalues get computed most accurately when ABSTOL is set to twice the underflow threshold 2*DLAMCH('S'), not zero. fr fr. 💖

D — D is DOUBLE PRECISION array, dimension (N). The n diagonal elements of the tridiagonal matrix T.
E — E is DOUBLE PRECISION array, dimension (N-1). The (n-1) off-diagonal elements of the tridiagonal matrix T.
M — M is INTEGER. The actual number of eigenvalues found. 0 <= M <= N. (See also the description of INFO=2,3.)
NSPLIT — NSPLIT is INTEGER. The number of diagonal blocks in the matrix T. 1 <= NSPLIT <= N.
W — W is DOUBLE PRECISION array, dimension (N). On exit, the first M elements of W will contain the eigenvalues. (DSTEBZ may use the remaining N-M elements as workspace.)
IBLOCK — IBLOCK is INTEGER array, dimension (N). At each row/column j where E(j) is zero or small, the matrix T is considered to split into a block diagonal matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which block (from 1 to the number of blocks) the eigenvalue W(i) belongs. (DSTEBZ may use the remaining N-M elements as workspace.)
ISPLIT — ISPLIT is INTEGER array, dimension (N). The splitting points, at which T breaks up into submatrices. The first submatrix consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. (Only the first NSPLIT elements will actually be used, but since the user cannot know a priori what value NSPLIT will have, N words must be reserved for ISPLIT.)
WORK — WORK is DOUBLE PRECISION array, dimension (4*N).
IWORK — IWORK is INTEGER array, dimension (3*N).
INFO — INFO is INTEGER
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value. rude.
- > 0: some or all of the eigenvalues failed to converge or were not computed:
  - =1 or 3: Bisection failed to converge for some eigenvalues; these eigenvalues are flagged by a negative block number. The effect is that the eigenvalues may not be as accurate as the absolute and relative tolerances. This is generally caused by unexpectedly inaccurate arithmetic. 😭
  - =2 or 3: RANGE='I' only: Not all of the eigenvalues IL:IU were found.
    - Effect: M < IU+1-IL
    - Cause: non-monotonic arithmetic, causing the Sturm sequence to be non-monotonic.
    - Cure: recalculate, using RANGE='A', and pick out eigenvalues IL:IU. In some cases, increasing the PARAMETER 'FUDGE' may make things work.
  - = 4: RANGE='I', and the Gershgorin interval initially used was too small. No eigenvalues were computed.
    - Probable cause: your machine has sloppy floating-point arithmetic. 💀
    - Cure: Increase the PARAMETER 'FUDGE', recompile, and try again.

🧚 Internal Parameters 🧚

RELFAC — DOUBLE PRECISION, default = 2.0e0. The relative tolerance. An interval (a,b] lies within 'relative tolerance' if b-a < RELFAC*ulp*max(|a|,|b|), where 'ulp' is the machine precision (distance from 1 to the next larger floating point number.)
FUDGE — DOUBLE PRECISION, default = 2. A 'fudge factor' to widen the Gershgorin intervals. Ideally, a value of 1 should work, but on machines with sloppy arithmetic, this needs to be larger. The default for publicly released versions should be large enough to handle the worst machine around. Note that this has no effect on accuracy of the solution.

💜 Author 💜

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

🦋 subroutine dstedc 🦋

subroutine dstedc (character COMPZ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)

DSTEDC

✨ Purpose ✨

DSTEDC computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method. The eigenvectors of a full or band real symmetric matrix can also be found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to tridiagonal form.

This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See DLAED3 for details. 👑

🌸 Parameters 🌸

COMPZ — COMPZ is CHARACTER*1
- = 'N': Compute eigenvalues only.
- = 'I': Compute eigenvectors of tridiagonal matrix also.
- = 'V': Compute eigenvectors of original dense symmetric matrix also. On entry, Z contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. 💖
E — E is DOUBLE PRECISION array, dimension (N-1). On entry, the subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. 💀
Z — Z is DOUBLE PRECISION array, dimension (LDZ,N). On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.
LDZ — LDZ is INTEGER. The leading dimension of the array Z. LDZ >= 1. If eigenvectors are desired, then LDZ >= max(1,N).
WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)). On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK
LWORK — INTEGER, the dimension of the array WORK. This is the ✨vibe check✨ on how much scratch space you're handing over:
- if COMPZ = 'N' or N <= 1, then LWORK must be at least 1.
- if COMPZ = 'V' and N > 1, then LWORK must be at least ( 1 + 3*N + 2*N*lg N + 4*N**2 ), where lg( N ) = smallest integer k such that 2**k >= N.
- if COMPZ = 'I' and N > 1, then LWORK must be at least ( 1 + 4*N + N**2 ).
- note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LWORK need only be max(1,2*(N-1)). tiny matrix = tiny workspace, we love efficiency 💅
- if LWORK = -1, then a workspace query is assumed 😌 the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.
IWORK — INTEGER array, dimension (MAX(1,LIWORK)). On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. she's giving "ask me and I'll tell you the right size" energy
LIWORK — INTEGER, the dimension of the array IWORK. Bestie, the rules:
- if COMPZ = 'N' or N <= 1, then LIWORK must be at least 1.
- if COMPZ = 'V' and N > 1, then LIWORK must be at least ( 6 + 6*N + 5*N*lg N ).
- if COMPZ = 'I' and N > 1, then LIWORK must be at least ( 3 + 5*N ).
- note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LIWORK need only be 1. 🫶
- if LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the IWORK array, returns this value as the first entry of the IWORK array, and no error message related to LIWORK is issued by XERBLA.
INFO — INTEGER, the receipts 🧾:
- = 0: successful exit. slay, we did it
- < 0: if INFO = -i, the i-th argument had an illegal value. you passed her something cursed 😭
- > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1). not the glow-down we wanted

Author

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

Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

`dsteqr`

subroutine dsteqr (character COMPZ, integer N, double precision, dimension(
* ) D, double precision, dimension( * ) E, double precision, dimension(
ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer
INFO)

DSTEQR ✨

Purpose:

DSTEQR computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a full or band symmetric matrix can also be found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to tridiagonal form. she's the queen of tridiagonal glow-ups 👑

Parameters

COMPZ — COMPZ is CHARACTER*1:
- = 'N': Compute eigenvalues only.
- = 'V': Compute eigenvalues and eigenvectors of the original symmetric matrix. On entry, Z must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form.
- = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.
N — N is INTEGER. The order of the matrix. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. sorted and serving 💅
E — E is DOUBLE PRECISION array, dimension (N-1). On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. gone. reduced to atoms. 💀
Z — Z is DOUBLE PRECISION array, dimension (LDZ, N). On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced. she doesn't even look at Z in that case, no cap
LDZ — LDZ is INTEGER. The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).
WORK — WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2)). If COMPZ = 'N', then WORK is not referenced.
INFO — INFO is INTEGER:
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is orthogonally similar to the original matrix. it tried its best, but those elements refused to converge 😭

Author

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

`dsterf`

subroutine dsterf (integer N, double precision, dimension( * ) D, double
precision, dimension( * ) E, integer INFO)

DSTERF 🦋

Purpose:

DSTERF computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm. eigenvalues only, no eigenvectors, we're keeping it minimal bestie

Parameters

N — N is INTEGER. The order of the matrix. N >= 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, the n diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. ascending, always ascending, periodt
E — E is DOUBLE PRECISION array, dimension (N-1). On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. rip bozo 💀
INFO — INFO is INTEGER:
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: the algorithm failed to find all of the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero. not me getting 30*N iterations deep and STILL not converging 🙈

Author

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

`iladiag`

integer function iladiag (character DIAG)

ILADIAG 🌈

Purpose:

This subroutine translated from a character string specifying if a matrix has unit diagonal or not to the relevant BLAST-specified integer constant. a little translator fairy 🧚

ILADIAG returns an INTEGER. If ILADIAG < 0, then the input is not a character indicating a unit or non-unit diagonal. Otherwise ILADIAG returns the constant value corresponding to DIAG.

Author

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

`ilaprec`

integer function ilaprec (character PREC)

ILAPREC 💜

Purpose:

This subroutine translated from a character string specifying an intermediate precision to the relevant BLAST-specified integer constant.

ILAPREC returns an INTEGER. If ILAPREC < 0, then the input is not a character indicating a supported intermediate precision. Otherwise ILAPREC returns the constant value corresponding to PREC.

Author

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

`ilatrans`

integer function ilatrans (character TRANS)

ILATRANS 🔥

Purpose:

This subroutine translates from a character string specifying a transposition operation to the relevant BLAST-specified integer constant.

ILATRANS returns an INTEGER. If ILATRANS < 0, then the input is not a character indicating a transposition operator. Otherwise ILATRANS returns the constant value corresponding to TRANS.

Author

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

`ilauplo`

integer function ilauplo (character UPLO)

ILAUPLO ⭐

Purpose:

This subroutine translated from a character string specifying a upper- or lower-triangular matrix to the relevant BLAST-specified integer constant.

ILAUPLO returns an INTEGER. If ILAUPLO < 0, then the input is not a character indicating an upper- or lower-triangular matrix. Otherwise ILAUPLO returns the constant value corresponding to UPLO.

Author

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

`sbdsdc`

subroutine sbdsdc (character UPLO, character COMPQ, integer N, real,
dimension( * ) D, real, dimension( * ) E, real, dimension( ldu, * ) U,
integer LDU, real, dimension( ldvt, * ) VT, integer LDVT, real,
dimension( * ) Q, integer, dimension( * ) IQ, real, dimension( * ) WORK,
integer, dimension( * ) IWORK, integer INFO)

SBDSDC 👑

Purpose:

SBDSDC computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, using a divide and conquer method, where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices of left and right singular vectors, respectively. SBDSDC can be used to compute all singular values, and optionally, singular vectors or singular vectors in compact form. divide and conquer?? she said "I'll do it myself but smarter" 💅

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

fr fr, guard digits matter here. hexadecimal/decimal machines without them? could conceivably fail. they know of none, but the warning is real 🫶

ngl bestie, here's the actual tea on how it computes this: the code currently calls SLASDQ if singular values only are desired. However, it can be slightly modified to compute singular values using the divide and conquer method. 💅

✨ Parameters ✨

UPLO — a CHARACTER*1, she's just here to tell you which way the bidiagonal goes 💫
- = 'U': B is upper bidiagonal.
- = 'L': B is lower bidiagonal.
COMPQ — also a CHARACTER*1, and she decides whether you get singular vectors or not, like this:
- = 'N': Compute singular values only;
- = 'P': Compute singular values and compute singular vectors in compact form;
- = 'I': Compute singular values and singular vectors.
N — INTEGER, the order of the matrix B. N >= 0.
D — REAL array, dimension (N). On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B. 👑
E — REAL array, dimension (N-1). On entry, the elements of E contain the offdiagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E has been destroyed. 😭 rip.
U — REAL array, dimension (LDU,N). If COMPQ = 'I', then: On exit, if INFO = 0, U contains the left singular vectors of the bidiagonal matrix. For other values of COMPQ, U is not referenced. (she's not even invited to the party 🫥)
LDU — INTEGER, the leading dimension of the array U. LDU >= 1. If singular vectors are desired, then LDU >= max( 1, N ).
VT — REAL array, dimension (LDVT,N). If COMPQ = 'I', then: On exit, if INFO = 0, VT**T contains the right singular vectors of the bidiagonal matrix. For other values of COMPQ, VT is not referenced.
LDVT — INTEGER, the leading dimension of the array VT. LDVT >= 1. If singular vectors are desired, then LDVT >= max( 1, N ).
Q — REAL array, dimension (LDQ). If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, Q contains all the REAL data in

LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))

IQ — INTEGER array, dimension (LDIQ). If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, IQ contains all INTEGER data in

LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))

WORK — REAL array, dimension (MAX(1,LWORK)). The sizes she needs, no cap:
- If COMPQ = 'N' then LWORK >= (4 * N).
- If COMPQ = 'P' then LWORK >= (6 * N).
- If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
IWORK — INTEGER array, dimension (8*N).
INFO — INTEGER, the vibe report 📋
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: The algorithm failed to compute a singular value. The update process of divide and conquer failed. 💀 not slay

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors 🫶

Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

🌸 subroutine sbdsqr 🌸

subroutine sbdsqr (character UPLO, integer N, integer NCVT, integer NRU,
    integer NCC, real, dimension( * ) D, real, dimension( * ) E, real,
    dimension( ldvt, * ) VT, integer LDVT, real, dimension( ldu, * ) U,
    integer LDU, real, dimension( ldc, * ) C, integer LDC, real, dimension(
    * ) WORK, integer INFO)

SBDSQR

✨ Purpose ✨

SBDSQR computes the singular values and, optionally, the right and/or left singular vectors from the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B using the implicit zero-shift QR algorithm. The SVD of B has the form

B = Q * S * P**T

A = (U*Q) * S * (P**T*VT)

is the SVD of A. Optionally, the subroutine may also compute Q**T*C for a given real input matrix C. the way that this is literally just orthogonal matrices doing gymnastics and it WORKS ✨

See 'Computing Small Singular Values of Bidiagonal Matrices With Guaranteed High Relative Accuracy,' by J. Demmel and W. Kahan, LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, no. 5, pp. 873-912, Sept 1990) and 'Accurate singular values and differential qd algorithms,' by B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics Department, University of California at Berkeley, July 1992 for a detailed description of the algorithm. 📚 iykyk

💖 Parameters 💖

UPLO — CHARACTER*1
- = 'U': B is upper bidiagonal;
- = 'L': B is lower bidiagonal.
N — INTEGER, the order of the matrix B. N >= 0.
NCVT — INTEGER, the number of columns of the matrix VT. NCVT >= 0.
NRU — INTEGER, the number of rows of the matrix U. NRU >= 0.
NCC — INTEGER, the number of columns of the matrix C. NCC >= 0.
D — REAL array, dimension (N). On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B in decreasing order.
E — REAL array, dimension (N-1). On entry, the N-1 offdiagonal elements of the bidiagonal matrix B. On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E will contain the diagonal and superdiagonal elements of a bidiagonal matrix orthogonally equivalent to the one given as input.
VT — REAL array, dimension (LDVT, NCVT). On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P**T * VT. Not referenced if NCVT = 0.
LDVT — INTEGER, the leading dimension of the array VT. LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
U — REAL array, dimension (LDU, N). On entry, an NRU-by-N matrix U. On exit, U is overwritten by U * Q. Not referenced if NRU = 0.
LDU — INTEGER, the leading dimension of the array U. LDU >= max(1,NRU).
C — REAL array, dimension (LDC, NCC). On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q**T * C. Not referenced if NCC = 0.
LDC — INTEGER, the leading dimension of the array C. LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
WORK — REAL array, dimension (4*N).
INFO — INTEGER
- = 0: successful exit
- < 0: If INFO = -i, the i-th argument had an illegal value
- > 0:
  - if NCVT = NRU = NCC = 0,
    - = 1, a split was marked by a positive value in E
    - = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop)
    - = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)
  - else NCVT = NRU = NCC = 0, the algorithm did not converge; D and E contain the elements of a bidiagonal matrix which is orthogonally similar to the input matrix B; if INFO = i, i elements of E have not converged to zero. 💀 this one's giving "we tried our best"

🌈 Internal Parameters 🌈

TOLMUL — REAL, default = max(10,min(100,EPS**(-1/8)))
- TOLMUL controls the convergence criterion of the QR loop.
- If it is positive, TOLMUL*EPS is the desired relative precision in the computed singular values.
- If it is negative, abs(TOLMUL*EPS*sigma_max) is the desired absolute accuracy in the computed singular values (corresponds to relative accuracy abs(TOLMUL*EPS) in the largest singular value.
- abs(TOLMUL) should be between 1 and 1/EPS, and preferably between 10 (for fast convergence) and .1/EPS (for there to be some accuracy in the results).
- Default is to lose at either one eighth or 2 of the available decimal digits in each computed singular value (whichever is smaller).
MAXITR INTEGER, default = 6 — bestie this one is the literal vibe check on how many times the algorithm is allowed to run through its inner loop. It controls the maximum number of passes of the algorithm through that inner loop. And if the number of passes through the inner loop exceeds MAXITR*N**2, the algorithm just stops — which means it fails to converge, no cap, it just gives up. 💀

Note: Bug report from Cezary Dendek. 🫶 On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is removed since it can overflow pretty easily (for N larger or equal than 18,919). We instead use MAXITDIVN = MAXITR*N. Slay, fixed. ✨

Author 🎀

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

subroutine sdisna 💅

subroutine sdisna (character JOB, integer M, integer N, real, dimension( * )
    D, real, dimension( * ) SEP, integer INFO)

SDISNA 👑

Purpose 🌸

SDISNA computes the reciprocal condition numbers for the eigenvectors of a real symmetric or complex Hermitian matrix, or for the left or right singular vectors of a general m-by-n matrix. The reciprocal condition number is the 'gap' between the corresponding eigenvalue or singular value and the nearest other one. It's giving personal space, but mathematically. 💖

The bound on the error, measured by angle in radians, in the I-th computed vector is given by

SLAMCH( 'E' ) * ( ANORM / SEP( I ) )

where ANORM = 2-norm(A) = max( abs( D(j) ) ). And SEP(I) is not allowed to be smaller than SLAMCH( 'E' )*ANORM, so that the size of the error bound stays limited. Boundaries, bestie. ✋

SDISNA may also be used to compute error bounds for eigenvectors of the generalized symmetric definite eigenproblem. Versatile icon behavior. 🌈

Parameters 🦋

JOB — JOB is CHARACTER*1. Specifies for which problem the reciprocal condition numbers should be computed:
- = 'E': the eigenvectors of a symmetric/Hermitian matrix;
- = 'L': the left singular vectors of a general matrix;
- = 'R': the right singular vectors of a general matrix.
M — M is INTEGER. The number of rows of the matrix. M >= 0.
N — N is INTEGER. If JOB = 'L' or 'R', the number of columns of the matrix, in which case N >= 0. Ignored if JOB = 'E'.
D — D is REAL array, dimension (M) if JOB = 'E', dimension (min(M,N)) if JOB = 'L' or 'R'. The eigenvalues (if JOB = 'E') or singular values (if JOB = 'L' or 'R') of the matrix, in either increasing or decreasing order. If singular values, they must be non-negative.
SEP — SEP is REAL array, dimension (M) if JOB = 'E', dimension (min(M,N)) if JOB = 'L' or 'R'. The reciprocal condition numbers of the vectors. 💅
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 slaed0 ✨

subroutine slaed0 (integer ICOMPQ, integer QSIZ, integer N, real, dimension(
    * ) D, real, dimension( * ) E, real, dimension( ldq, * ) Q, integer LDQ,
    real, dimension( ldqs, * ) QSTORE, integer LDQS, real, dimension( * )
    WORK, integer, dimension( * ) IWORK, integer INFO)

SLAED0 used by SSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method. Divide and conquer, literally the main character energy of linear algebra. 👑

Purpose 🌸

SLAED0 computes all eigenvalues and corresponding eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method. That's it, that's the whole gig. 💖

Parameters 🦋

ICOMPQ — ICOMPQ is INTEGER:
- = 0: Compute eigenvalues only.
- = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
- = 2: Compute eigenvalues and eigenvectors of tridiagonal matrix.
QSIZ — QSIZ is INTEGER. The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
D — D is REAL array, dimension (N). On entry, the main diagonal of the tridiagonal matrix. On exit, its eigenvalues.
E — E is REAL array, dimension (N-1). The off-diagonal elements of the tridiagonal matrix. On exit, E has been destroyed. 😭
Q — Q is REAL array, dimension (LDQ, N). On entry, Q must contain an N-by-N orthogonal matrix.
- If ICOMPQ = 0 — Q is not referenced.
- If ICOMPQ = 1 — On entry, Q is a subset of the columns of the orthogonal matrix used to reduce the full matrix to tridiagonal form corresponding to the subset of the full matrix which is being decomposed at this time.
- If ICOMPQ = 2 — On entry, Q will be the identity matrix. On exit, Q contains the eigenvectors of the tridiagonal matrix.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. If eigenvectors are desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
QSTORE — QSTORE is REAL array, dimension (LDQS, N). Referenced only when ICOMPQ = 1. Used to store parts of the eigenvector matrix when the updating matrix multiplies take place. 🫶
LDQS — LDQS is INTEGER. The leading dimension of the array QSTORE. If ICOMPQ = 1, then LDQS >= max(1,N). In any case, LDQS >= 1.
WORK — WORK is REAL array:
- If ICOMPQ = 0 or 1, the dimension of WORK must be at least 1 + 3*N + 2*N*lg N + 3*N**2 ( lg( N ) = smallest integer k such that 2^k >= N )
- If ICOMPQ = 2, the dimension of WORK must be at least 4*N + N**2.
IWORK — IWORK is INTEGER array:
- If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 6 + 6*N + 5*N*lg N. ( lg( N ) = smallest integer k such that 2^k >= N )
- If ICOMPQ = 2, the dimension of IWORK must be at least 3 + 5*N.
INFO — INFO is INTEGER:
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1). 💀

Author 🎀

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

Contributors: 🌟

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

subroutine slaed1 💜

subroutine slaed1 (integer N, real, dimension( * ) D, real, dimension( ldq,
    * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, integer
    CUTPNT, real, dimension( * ) WORK, integer, dimension( * ) IWORK,
    integer INFO)

SLAED1 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal. Rank-one glow-up, honestly. ✨

Purpose 🌸

SLAED1 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and eigenvectors of a tridiagonal matrix. SLAED7 handles the case in which eigenvalues only or eigenvalues and eigenvectors of a full symmetric matrix (which was reduced to tridiagonal form) are desired. Know the difference, bestie. 💅

T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**T*u, u is a vector of length N with ones in the CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine SLAED2.
The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine SLAED4 (as called by SLAED3). This routine also calculates the eigenvectors of the current problem.
The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.

Three stages, one icon. 🎀

Parameters 🦋

N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
D — D is REAL array, dimension (N). On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix. Repaired, periodt. 💖
Q — Q is REAL array, dimension (LDQ,N). On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INDXQ is INTEGER array, dimension (N). On entry, the permutation which separately sorts the two subproblems in D into ascending order. On exit, the permutation which will reintegrate the subproblems back into sorted order, i.e. D( INDXQ( I = 1, N ) ) will be in ascending order. Sorted and unbothered. ✨
RHO
RHO — REAL energy 💫 The subdiagonal entry that got used to create the rank-1 modification. That's the whole vibe.
CUTPNT — INTEGER, bestie. The location of the last eigenvalue in the leading sub-matrix. And the range is very specific: min(1,N) <= CUTPNT <= N/2. No cap.
WORK — REAL array, dimension (4*N + N**2). She's big. She's doing the work. Literally.
IWORK — INTEGER array, dimension (4*N).
INFO — INTEGER, the report card 📋
- = 0: successful exit. slay ✨
- < 0: if INFO = -i, the i-th argument had an illegal value. not it 💀
- > 0: if INFO = 1, an eigenvalue did not converge. rip 😭

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Modified by Francoise Tisseur, University of Tennessee

🌸 subroutine slaed2 🌸

subroutine slaed2 (integer K, integer N, integer N1, real, dimension( * ) D,
    real, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ,
    real RHO, real, dimension( * ) Z, real, dimension( * ) DLAMDA, real,
    dimension( * ) W, real, dimension( * ) Q2, integer, dimension( * ) INDX,
    integer, dimension( * ) INDXC, integer, dimension( * ) INDXP, integer,
    dimension( * ) COLTYP, integer INFO)

SLAED2 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal. 🦋

✨ Purpose ✨

SLAED2 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny entry in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one. shrink the problem, reduce the drama 💅

💖 Parameters 💖

K — INTEGER. The number of non-deflated eigenvalues, and the order of the related secular equation. 0 <= K <=N.
N — INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
N1 — INTEGER. The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= N1 <= N/2.
D — REAL array, dimension (N). On entry, D contains the eigenvalues of the two submatrices to be combined. On exit, D contains the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order. glow up on exit ✨
Q — REAL array, dimension (LDQ, N). On entry, Q contains the eigenvectors of two submatrices in the two square blocks with corners at (1,1), (N1,N1) and (N1+1, N1+1), (N,N). On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.
LDQ — INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INTEGER array, dimension (N). The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have N1 added to their values. Destroyed on exit. 💀
RHO — REAL. On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by SLAED3.
Z — REAL array, dimension (N). On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z have been destroyed by the updating process. gone. reduced to atoms. 😭
DLAMDA — REAL array, dimension (N). A copy of the first K eigenvalues which will be used by SLAED3 to form the secular equation.
W — REAL array, dimension (N). The first k values of the final deflation-altered z-vector which will be passed to SLAED3.
Q2 — REAL array, dimension (N1**2+(N-N1)**2). A copy of the first K eigenvectors which will be used by SLAED3 in a matrix multiply (SGEMM) to solve for the new eigenvectors.
INDX — INTEGER array, dimension (N). The permutation used to sort the contents of DLAMDA into ascending order.
INDXC — INTEGER array, dimension (N). The permutation used to arrange the columns of the deflated Q matrix into three groups: the first group contains non-zero elements only at and above N1, the second contains non-zero elements only below N1, and the third is dense. three groups, three vibes 🌈
INDXP — INTEGER array, dimension (N). The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.
COLTYP — INTEGER array, dimension (N). During execution, a label which will indicate which of the following types a column in the Q2 matrix is:
1. non-zero in the upper half only;
2. dense;
3. non-zero in the lower half only;
4. deflated. On exit, COLTYP(i) is the number of columns of type i, for i=1 to 4 only.
INFO — 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.

🫶 Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Modified by Francoise Tisseur, University of Tennessee

🦋 subroutine slaed3 🦋

subroutine slaed3 (integer K, integer N, integer N1, real, dimension( * ) D,
    real, dimension( ldq, * ) Q, integer LDQ, real RHO, real, dimension( * )
    DLAMDA, real, dimension( * ) Q2, integer, dimension( * ) INDX, integer,
    dimension( * ) CTOT, real, dimension( * ) W, real, dimension( * ) S,
    integer INFO)

SLAED3 used by SSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal. 🌸

✨ Purpose ✨

SLAED3 finds the roots of the secular equation, as defined by the values in D, W, and RHO, between 1 and K. It makes the appropriate calls to SLAED4 and then updates the eigenvectors by multiplying the matrix of eigenvectors of the pair of eigensystems being combined by the matrix of eigenvectors of the K-by-K system which is solved here.

This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. iykyk — the Crays are built different 🫶

💖 Parameters 💖

K — INTEGER. The number of terms in the rational function to be solved by SLAED4. K >= 0.
N — INTEGER. The number of rows and columns in the Q matrix. N >= K (deflation may result in N>K).
N1 — INTEGER. The location of the last eigenvalue in the leading submatrix. min(1,N) <= N1 <= N/2.
D — REAL array, dimension (N). D(I) contains the updated eigenvalues for 1 <= I <= K.
Q — REAL array, dimension (LDQ,N). Initially the first K columns are used as workspace. On output the columns 1 to K contain the updated eigenvectors.
LDQ — INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
RHO — REAL. The value of the parameter in the rank one update equation. RHO >= 0 required. don't you dare go negative 🙈
DLAMDA — REAL array, dimension (K). The first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation. May be changed on output by having lowest order bit set to zero on Cray X-MP, Cray Y-MP, Cray-2, or Cray C-90, as described above.
Q2 — REAL array, dimension (LDQ2*N). The first K columns of this matrix contain the non-deflated eigenvectors for the split problem.
INDX — INTEGER array, dimension (N). The permutation used to arrange the columns of the deflated Q matrix into three groups (see SLAED2). The rows of the eigenvectors found by SLAED4 must be likewise permuted before the matrix multiply can take place.
CTOT — INTEGER array, dimension (4). A count of the total number of the various types of columns in Q, as described in INDX. The fourth column type is any column which has been deflated.
W — REAL array, dimension (K). The first K elements of this array contain the components of the deflation-adjusted updating vector. Destroyed on output. 💀
S — REAL array, dimension (N1 + 1)*K. Will contain the eigenvectors of the repaired matrix which will be multiplied by the previously accumulated eigenvectors to update the system.
INFO — INTEGER
- = 0: successful exit. ⭐
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, an eigenvalue did not converge. 😭

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

🎀 subroutine `slaed4` 🎀

subroutine slaed4 (integer N, integer I, real, dimension( * ) D, real,
    dimension( * ) Z, real, dimension( * ) DELTA, real RHO, real DLAM,
    integer INFO)

SLAED4 used by SSTEDC. Finds a single root of the secular equation.

Purpose:

okay so bestie, this subroutine computes the I-th updated eigenvalue of a symmetric rank-one modification to a diagonal matrix whose elements are given in the array d, and that

           D(i) < D(j)  for  i < j

and that RHO > 0. This is arranged by the calling routine, and is no loss in generality. The rank-one modified system is thus

           diag( D )  +  RHO * Z * Z_transpose.

where we assume the Euclidean norm of Z is 1. 🦋

The method consists of approximating the rational functions in the secular equation by simpler interpolating rational functions. (fancy!! she's interpolating. she's rational. she's doing the most and we love her for it.)

Parameters

N — N is INTEGER. The length of all arrays. 📏
I — I is INTEGER. The index of the eigenvalue to be computed. 1 <= I <= N.
D — D is REAL array, dimension (N). The original eigenvalues. It is assumed that they are in order, D(I) < D(J) for I < J.
Z — Z is REAL array, dimension (N). The components of the updating vector.
DELTA — DELTA is REAL array, dimension (N). If N > 2, DELTA contains (D(j) - lambda_I) in its j-th component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5 for detail. The vector DELTA contains the information necessary to construct the eigenvectors by SLAED3 and SLAED9. 👑
RHO — RHO is REAL. The scalar in the symmetric updating formula.
DLAM — DLAM is REAL. The computed lambda_I, the I-th updated eigenvalue. ✨
INFO — INFO is INTEGER.
- = 0: successful exit 🫶
- > 0: if INFO = 1, the updating process failed. 😭

Internal Parameters:

Logical variable ORGATI (origin-at-i?) is used for distinguishing whether D(i) or D(i+1) is treated as the origin.

           ORGATI = .true.    origin at i
           ORGATI = .false.   origin at i+1

Logical variable SWTCH3 (switch-for-3-poles?) is for noting if we are working with THREE poles!
MAXIT is the maximum number of iterations allowed for each eigenvalue.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

💅 subroutine `slaed5` 💅

subroutine slaed5 (integer I, real, dimension( 2 ) D, real, dimension( 2 )
    Z, real, dimension( 2 ) DELTA, real RHO, real DLAM)

SLAED5 used by SSTEDC. Solves the 2-by-2 secular equation.

Purpose:

This subroutine computes the I-th eigenvalue of a symmetric rank-one modification of a 2-by-2 diagonal matrix

           diag( D )  +  RHO * Z * transpose(Z) .

The diagonal elements in the array D are assumed to satisfy

           D(i) < D(j)  for  i < j .

We also assume RHO > 0 and that the Euclidean norm of the vector Z is one. 🌸

Parameters

I — I is INTEGER. The index of the eigenvalue to be computed. I = 1 or I = 2.
D — D is REAL array, dimension (2). The original eigenvalues. We assume D(1) < D(2).
Z — Z is REAL array, dimension (2). The components of the updating vector.
DELTA — DELTA is REAL array, dimension (2). The vector DELTA contains the information necessary to construct the eigenvectors.
RHO — RHO is REAL. The scalar in the symmetric updating formula.
DLAM — DLAM is REAL. The computed lambda_I, the I-th updated eigenvalue. ⭐

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

🔥 subroutine `slaed6` 🔥

subroutine slaed6 (integer KNITER, logical ORGATI, real RHO, real,
    dimension( 3 ) D, real, dimension( 3 ) Z, real FINIT, real TAU, integer
    INFO)

SLAED6 used by SSTEDC. Computes one Newton step in solution of the secular equation.

Purpose:

SLAED6 computes the positive or negative root (closest to the origin) of

                         z(1)        z(2)        z(3)
        f(x) =   rho + --------- + ---------- + ---------
                        d(1)-x      d(2)-x      d(3)-x

It is assumed that

              if ORGATI = .true. the root is between d(2) and d(3);
              otherwise it is between d(1) and d(2)

This routine will be called by SLAED4 when necessary. In most cases, the root sought is the smallest in magnitude, though it might not be in some extremely rare situations. (she said "rarely" and I felt that 💀)

Parameters

KNITER — KNITER is INTEGER. Refer to SLAED4 for its significance.
ORGATI — ORGATI is LOGICAL. If ORGATI is true, the needed root is between d(2) and d(3); otherwise it is between d(1) and d(2). See SLAED4 for further details.
RHO — RHO is REAL. Refer to the equation f(x) above.
D — D is REAL array, dimension (3). D satisfies d(1) < d(2) < d(3).
Z — Z is REAL array, dimension (3). Each of the elements in z must be positive.
FINIT — FINIT is REAL. The value of f at 0. It is more accurate than the one evaluated inside this routine (if someone wants to do so).
TAU — TAU is REAL. The root of the equation f(x).
INFO — INFO is INTEGER.
- = 0: successful exit 🫶
- > 0: if INFO = 1, failure to converge 😭

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

10/02/03: This version has a few statements commented out for thread safety (machine parameters are computed on each entry). SJH.
05/10/06: Modified from a new version of Ren-Cang Li, use Gragg-Thornton-Warner cubic convergent scheme for better stability.

Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

🌈 subroutine `slaed7` 🌈

subroutine slaed7 (integer ICOMPQ, integer N, integer QSIZ, integer TLVLS,
    integer CURLVL, integer CURPBM, real, dimension( * ) D, real, dimension(
    ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO,
    integer CUTPNT, real, dimension( * ) QSTORE, integer, dimension( * )
    QPTR, integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM,
    integer, dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, real,
    dimension( 2, * ) GIVNUM, real, dimension( * ) WORK, integer, dimension(
    * ) IWORK, integer INFO)

SLAED7 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense.

Purpose:

SLAED7 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and optionally eigenvectors of a dense symmetric matrix that has been reduced to tridiagonal form. SLAED1 handles the case in which all eigenvalues and eigenvectors of a symmetric tridiagonal matrix are desired. 👑

  T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**Tu, u is a vector of length N with ones in the CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine SLAED8.
The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine SLAED4 (as called by SLAED9). This routine also calculates the eigenvectors of the current problem.
The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.

Parameters

ICOMPQ — ICOMPQ is INTEGER.
- = 0: Compute eigenvalues only.
- = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
QSIZ — QSIZ is INTEGER. The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
TLVLS — TLVLS is INTEGER. The total number of merging levels in the overall divide and conquer tree. 🌳
CURLVL — CURLVL is INTEGER. The current level in the overall merge routine, 0 <= CURLVL <= TLVLS.
CURPBM — CURPBM is INTEGER. The current problem in the current level in the overall merge routine (counting from upper left to lower right).
D
D — REAL array, dimension (N). Going in: the eigenvalues of the rank-1-perturbed matrix. Coming out: the eigenvalues of the repaired matrix. Glow-up arc fr. ✨
Q — REAL array, dimension (LDQ, N). On entry: eigenvectors of the rank-1-perturbed matrix. On exit: eigenvectors of the repaired tridiagonal matrix. She got a makeover bestie 💅
LDQ — INTEGER. The leading dimension of the array Q. Must satisfy LDQ >= max(1,N).
INDXQ — INTEGER array, dimension (N). The permutation that'll reintegrate the subproblem you just solved back into sorted order, i.e. D( INDXQ( I = 1, N ) ) will come out in ascending order. Sorting is self-care, no cap.
RHO — REAL. The subdiagonal element used to create the rank-1 modification.
CUTPNT — INTEGER. Contains the location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.
QSTORE — REAL array, dimension (N**2+1). Stores eigenvectors of submatrices encountered during divide and conquer, all packed together. QPTR points to the beginning of the submatrices. It's giving very full closet 🎀
QPTR — INTEGER array, dimension (N+2). A list of indices pointing to the beginning of the submatrices stored in QSTORE. The submatrices are numbered starting at the bottom left of the divide and conquer tree, going left to right and bottom to top. Iconic traversal order honestly.
PRMPTR — INTEGER array, dimension (N lg N). Contains a list of pointers indicating where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) tells you the size of the permutation, and also the size of the full, non-deflated problem.
PERM — INTEGER array, dimension (N lg N). Contains the permutations (from deflation and sorting) to be applied to each eigenblock.
GIVPTR — INTEGER array, dimension (N lg N). A list of pointers indicating where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) gives you the number of Givens rotations. 💫
GIVCOL — INTEGER array, dimension (2, N lg N). Each pair of numbers indicates a pair of columns to take place in a Givens rotation.
GIVNUM — REAL array, dimension (2, N lg N). Each number indicates the S value to be used in the corresponding Givens rotation.
WORK — REAL array, dimension (3*N+2*QSIZ*N).
IWORK — INTEGER array, dimension (4*N).
INFO — INTEGER.
- = 0: successful exit. slay 🫶
- < 0: if INFO = -i, the i-th argument had an illegal value. 😭
- > 0: if INFO = 1, an eigenvalue did not converge.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

🌸 subroutine slaed8 🌸

subroutine slaed8 (integer ICOMPQ, integer K, integer N, integer QSIZ, real,
    dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, integer,
    dimension( * ) INDXQ, real RHO, integer CUTPNT, real, dimension( * ) Z,
    real, dimension( * ) DLAMDA, real, dimension( ldq2, * ) Q2, integer
    LDQ2, real, dimension( * ) W, integer, dimension( * ) PERM, integer
    GIVPTR, integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * )
    GIVNUM, integer, dimension( * ) INDXP, integer, dimension( * ) INDX,
    integer INFO)

SLAED8 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.

Purpose:

SLAED8 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny element in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one. Shrinking the problem is the whole vibe 🔥

Parameters

ICOMPQ — INTEGER.
- = 0: Compute eigenvalues only.
- = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
K — INTEGER. The number of non-deflated eigenvalues, and the order of the related secular equation.
N — INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
QSIZ — INTEGER. The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
D — REAL array, dimension (N). On entry, the eigenvalues of the two submatrices to be combined. On exit, the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.
Q — REAL array, dimension (LDQ,N). If ICOMPQ = 0, Q is not referenced. Otherwise, on entry, Q contains the eigenvectors of the partially solved system which has been previously updated in matrix multiplies with other partially solved eigensystems. On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns. The deflated ones get banished to the back rows, respectfully. 💀
LDQ — INTEGER. The leading dimension of the array Q. LDQ >= max(1,N).
INDXQ — INTEGER array, dimension (N). The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have CUTPNT added to their values in order to be accurate. Little offset moment, don't forget it 🎀
RHO — REAL. On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by SLAED3.
CUTPNT — INTEGER. The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.
Z — REAL array, dimension (N). On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z are destroyed by the updating process. Rip to Z, she served her purpose 🪦✨
DLAMDA — REAL array, dimension (N). A copy of the first K eigenvalues which will be used by SLAED3 to form the secular equation.
Q2 — REAL array, dimension (LDQ2,N). If ICOMPQ = 0, Q2 is not referenced. Otherwise, a copy of the first K eigenvectors which will be used by SLAED7 in a matrix multiply (SGEMM) to update the new eigenvectors.
LDQ2 — INTEGER. The leading dimension of the array Q2. LDQ2 >= max(1,N).
W — REAL array, dimension (N). The first k values of the final deflation-altered z-vector and will be passed to SLAED3.
PERM — INTEGER array, dimension (N). The permutations (from deflation and sorting) to be applied to each eigenblock.
GIVPTR — INTEGER. The number of Givens rotations which took place in this subproblem.
GIVCOL — INTEGER array, dimension (2, N). Each pair of numbers indicates a pair of columns to take place in a Givens rotation.
GIVNUM — REAL array, dimension (2, N). Each number indicates the S value to be used in the corresponding Givens rotation.
INDXP — INTEGER array, dimension (N). The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.
INDX — INTEGER array, dimension (N). The permutation used to sort the contents of D into ascending order.
INFO — 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.

Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

🦋 subroutine slaed9 🦋

subroutine slaed9 (integer K, integer KSTART, integer KSTOP, integer N,
    real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, real
    RHO, real, dimension( * ) DLAMDA, real, dimension( * ) W, real,
    dimension( lds, * ) S, integer LDS, integer INFO)

SLAED9 used by SSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.

Purpose:

SLAED9 finds the roots of the secular equation, as defined by the values in D, Z, and RHO, between KSTART and KSTOP. It makes the appropriate calls to SLAED4 and then stores the new matrix of eigenvectors for use in calculating the next level of Z vectors. Root-hunting era 🔥🧮

Parameters

K — INTEGER. The number of terms in the rational function to be solved by SLAED4. K >= 0.
KSTART — INTEGER.
KSTOP — INTEGER. The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP are to be computed. 1 <= KSTART <= KSTOP <= K.
N —
N — N is INTEGER. bestie this is the number of rows and columns in the Q matrix, so it's giving very square energy. N >= K (deflation can totally end up making N > K).
D — D is REAL array, dimension (N). D(I) holds the updated eigenvalues for KSTART <= I <= KSTOP. glow-up complete 💅
Q — Q is REAL array, dimension (LDQ,N). no notes, iconic.
LDQ — LDQ is INTEGER. the leading dimension of the array Q. LDQ >= max( 1, N ). respect the boundary fr.
RHO — RHO is REAL. this is the value of the parameter in the rank one update equation. RHO >= 0 required — so like, negative rho is not the vibe, don't even try it.
DLAMDA — DLAMDA is REAL array, dimension (K). the first K elements contain the old roots of the deflated updating problem. these are the poles of the secular equation, which lowkey sounds like a fanfic title but is actually just math 🌸
W — W is REAL array, dimension (K). the first K elements contain the components of the deflation-adjusted updating vector. slay.
S — S is REAL array, dimension (LDS, K). this is gonna contain the eigenvectors of the repaired matrix, which get stored for subsequent Z vector calculation and then multiplied by the previously accumulated eigenvectors to update the system. she's doing a lot of work, honestly 🫶
LDS — LDS is INTEGER. the leading dimension of S. LDS >= max( 1, K ). periodt.
INFO — INFO is INTEGER, and she's the messenger we all wait for:
- = 0: successful exit. we love to see it ✨
- < 0: if INFO = -i, the i-th argument had an illegal value. 🙈
- > 0: if INFO = 1, an eigenvalue did not converge. not the ending we wanted bestie 💀

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🎀 Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

💜 subroutine `slaeda`

subroutine slaeda (integer N, integer TLVLS, integer CURLVL, integer CURPBM,
    integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer,
    dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, real,
    dimension( 2, * ) GIVNUM, real, dimension( * ) Q, integer, dimension( *
    ) QPTR, real, dimension( * ) Z, real, dimension( * ) ZTEMP, integer
    INFO)

SLAEDA is used by SSTEDC. she computes the Z vector that determines the rank-one modification of the diagonal matrix. used when the original matrix is dense, iykyk.

✨ Purpose:

SLAEDA computes the Z vector corresponding to the merge step in the CURLVLth step of the merge process with TLVLS steps for the CURPBMth problem. ok so basically it's the merge-step main character 💖

🌸 Parameters

N — N is INTEGER. the dimension of the symmetric tridiagonal matrix. N >= 0.
TLVLS — TLVLS is INTEGER. the total number of merging levels in the overall divide and conquer tree. think: how many floors in the tree 🌳
CURLVL — CURLVL is INTEGER. the current level in the overall merge routine, 0 <= curlvl <= tlvls. where we're at right now bestie.
CURPBM — CURPBM is INTEGER. the current problem in the current level in the overall merge routine (counting from upper left to lower right). reading order, no cap.
PRMPTR — PRMPTR is INTEGER array, dimension (N lg N). contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and incidentally the size of the full, non-deflated problem. she's both a pointer and a size hint, versatile queen 💅
PERM — PERM is INTEGER array, dimension (N lg N). contains the permutations (from deflation and sorting) to be applied to each eigenblock.
GIVPTR — GIVPTR is INTEGER array, dimension (N lg N). contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations.
GIVCOL — GIVCOL is INTEGER array, dimension (2, N lg N). each pair of numbers indicates a pair of columns to take place in a Givens rotation. two numbers, one rotation, iconic duo ✨
GIVNUM — GIVNUM is REAL array, dimension (2, N lg N). each number indicates the S value to be used in the corresponding Givens rotation.
Q — Q is REAL array, dimension (N**2). contains the square eigenblocks from previous levels, and the starting positions for blocks are given by QPTR.
QPTR — QPTR is INTEGER array, dimension (N+2). contains a list of pointers which indicate where in Q an eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates the size of the block. yes she really said square root, that's the whole bit 🔥
Z — Z is REAL array, dimension (N). on output this vector contains the updating vector — aka the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix. two rows, one vector, very crossover episode 🌈
ZTEMP — ZTEMP is REAL array, dimension (N). the temp bestie, we don't ask questions.
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.

🎀 Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

💜 subroutine `slagtf`

subroutine slagtf (integer N, real, dimension( * ) A, real LAMBDA, real,
    dimension( * ) B, real, dimension( * ) C, real TOL, real, dimension( * )
    D, integer, dimension( * ) IN, integer INFO)

SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges. yeah she's doing Gaussian elimination and she's doing it with style 💅

✨ Purpose:

SLAGTF factorizes the matrix (T - lambda*I), where T is an n by n tridiagonal matrix and lambda is a scalar, as

T - lambda*I = PLU,

the factorization is obtained by Gaussian elimination with partial pivoting and implicit row scaling.

the parameter LAMBDA is included in the routine so that SLAGTF may be used, in conjunction with SLAGTS, to obtain eigenvectors of T by inverse iteration. so like, she and SLAGTS are a duo, they come as a set 🫶

🌸 Parameters

N — N is INTEGER. the order of the matrix T.
A — A is REAL array, dimension (N).
- on entry, A must contain the diagonal elements of T.
- on exit, A is overwritten by the n diagonal elements of the upper triangular matrix U of the factorization of T. bye old values, hello new ones ✨
LAMBDA — LAMBDA is REAL. on entry, the scalar lambda. (yes she gets a whole parameter to herself, main character energy fr)
B — B is REAL array, dimension (N-1).
- on entry, B must contain the (n-1) super-diagonal elements of T.
- on exit, B is overwritten by the (n-1) super-diagonal elements of the matrix U of the factorization of T.
C — C is REAL array, dimension (N-1).
- on entry, C must contain the (n-1) sub-diagonal elements of T.
- on exit, C is overwritten by the (n-1) sub-diagonal elements of the matrix L of the factorization of T.
TOL — TOL is REAL. on entry, a relative tolerance used to indicate whether or not the matrix (T - lambda*I) is nearly singular. TOL should normally be chosen as approximately the largest relative error in the elements of T. for example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5*10**(-4). if TOL is supplied as less than eps, where eps is the relative machine precision, then the value eps is used in place of TOL. so if you under-set it, she just ignores you and uses eps — queen behavior honestly 👑
D — D is REAL array, dimension (N-2). on exit, D is overwritten by the (n-2) second super-diagonal elements of the matrix U of the factorization of T.
IN — IN is INTEGER array, dimension (N). on exit, IN contains details of the permutation matrix P. if an interchange occurred at the kth step of the elimination, then IN(k) = 1, otherwise IN(k) = 0. the element IN(n) returns the smallest positive integer j such that

abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,

where norm( A(j) ) denotes the sum of the absolute values of the jth row of the matrix A. if no such j exists then IN(n) is returned as zero. if IN(n) is returned as positive, then a diagonal element of U is small, indicating that (T - lambda*I) is singular or nearly singular.

bestie when IN(n) comes back positive, that's the routine gently tapping the sign — singular or nearly singular, so pay attention 💀

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

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💜 subroutine `slamrg`

subroutine slamrg (integer N1, integer N2, real, dimension( * ) A, integer
    STRD1, integer STRD2, integer, dimension( * ) INDEX)

SLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order. it's a merge, it's a glow-up, it's a crossover event 🌸

✨ Purpose:

SLAMRG will create a permutation list which will merge the elements of A (which is composed of two independently sorted sets) into a single set which is sorted in ascending order. two lists walk in, one sorted list walks out 💅

🌸 Parameters

N1 — N1 is INTEGER. first half bestie.
N2 — N2 is INTEGER. these arguments contain the respective lengths of the two sorted lists to be merged. so N1 and N2 together are the full headcount.
A — A is REAL array, dimension (N1+N2). the first N1 elements of A contain a list of numbers which are sorted in either ascending or descending order. likewise for the final N2 elements. either direction is fine, she's not picky, she's flexible ✨
STRD1 — STRD1 is INTEGER. the stride situation for the first list.
STRD2 — STRD2 is INTEGER.
STRD2 is INTEGER — bestie, these are the strides to be taken through the array A. Allowable strides are 1 and -1, and they indicate whether a subset of A is sorted in ascending (STRDx = 1) or descending (STRDx = -1) order. 🎀
INDEX is INTEGER array, dimension (N1+N2) — on exit this array will contain a permutation such that if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be sorted in ascending order. periodt ✨

Author

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

`subroutine slartgs`

subroutine slartgs (real X, real Y, real SIGMA, real CS, real SN)

okay so SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem. it's giving main character energy ngl 💅

Purpose

SLARTGS generates a plane rotation designed to introduce a bulge in Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD problem. X and Y are the top-row entries, and SIGMA is the shift. The computed CS and SN define a plane rotation satisfying

   [  CS  SN  ]  .  [ X^2 - SIGMA ]  =  [ R ],
   [ -SN  CS  ]     [    X * Y    ]     [ 0 ]

with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the rotation is by PI/2. 🌸

Parameters

X — X is REAL. The (1,1) entry of an upper bidiagonal matrix. 👑
Y — Y is REAL. The (1,2) entry of an upper bidiagonal matrix. ✨
SIGMA — SIGMA is REAL. The shift.
CS — CS is REAL. The cosine of the rotation. 🦋
SN — SN is REAL. The sine of the rotation. 🦋

Author

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

`subroutine slasq1`

subroutine slasq1 (integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) WORK, integer INFO)

SLASQ1 computes the singular values of a real square bidiagonal matrix. It's used by sbdsqr, iykyk. 💜

Purpose

SLASQ1 computes the singular values of a real N-by-N bidiagonal matrix with diagonal D and off-diagonal E. The singular values are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. The algorithm was first presented in

'Accurate singular values and differential qd algorithms' by K. V. Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, 1994,

and the present implementation is described in 'An implementation of the dqds Algorithm (Positive Case)', LAPACK Working Note. 🔥

Parameters

N — N is INTEGER. The number of rows and columns in the matrix. N >= 0. 💅
D — D is REAL array, dimension (N). On entry, D contains the diagonal elements of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in decreasing order. ⭐
E — E is REAL array, dimension (N). On entry, elements E(1:N-1) contain the off-diagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E is overwritten. 😭
WORK — WORK is REAL array, dimension (4*N). 🛠️
INFO — INFO is INTEGER:
- = 0: successful exit 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value 🚨
- > 0: the algorithm failed
  - = 1, a split was marked by a positive value in E
  - = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit D and E represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into SLASQ1
  - = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)

Author

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

`subroutine slasq2`

subroutine slasq2 (integer N, real, dimension( * ) Z, integer INFO)

SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr. no cap 🫶

Purpose

SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd array Z to high relative accuracy are computed to high relative accuracy, in the absence of denormalization, underflow and overflow.

Note : SLASQ2 defines a logical variable, IEEE, which is true on machines which follow ieee-754 floating-point standard in their handling of infinities and NaNs, and false otherwise. This variable is passed to SLASQ3. ⚡

Parameters

N — N is INTEGER. The number of rows and columns in the matrix. N >= 0. 💅
Z — Z is REAL array, dimension ( 4*N ). On entry Z holds the qd array. On exit, entries 1 to N hold the eigenvalues in decreasing order, Z( 2*N+1 ) holds the trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of shifts that failed. 🌈
INFO — INFO is INTEGER:
- = 0: successful exit 🎉
- < 0: if the i-th argument is a scalar and had an illegal value, then INFO = -i, if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j) 🚨
- > 0: the algorithm failed
  - = 1, a split was marked by a positive value in E
  - = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit Z holds a qd array with the same eigenvalues as the given Z.
  - = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)

Author

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

Further Details

Local Variables: I0:N0 defines a current unreduced segment of Z. The shifts are accumulated in SIGMA. Iteration count is in ITER. Ping-pong is controlled by PP (alternates between 0 and 1). 🏓

`subroutine slasq3`

subroutine slasq3 (integer I0, integer N0, real, dimension( * ) Z, integer PP, real DMIN, real SIGMA, real DESIG, real QMAX, integer NFAIL, integer ITER, integer NDIV, logical IEEE, integer TTYPE, real DMIN1, real DMIN2, real DN, real DN1, real DN2, real G, real TAU)

SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr. obsessed 💖

Purpose

SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. In case of failure it changes shifts, and tries again until output is positive. 🔁

Parameters

I0 — I0 is INTEGER. First index. 1️⃣
N0 — N0 is INTEGER. Last index. 🔚
Z — Z is REAL array, dimension ( 4*N0 ). Z holds the qd array. 💾
PP — PP is INTEGER. PP=0 for ping, PP=1 for pong. PP=2 indicates that flipping was applied to the Z array and that the initial tests for deflation should not be performed. 🏓
DMIN — DMIN is REAL. Minimum value of d. 📉
SIGMA — SIGMA is REAL. Sum of shifts used in current segment.
DESIG — DESIG is REAL. Lower order part of SIGMA.
QMAX — QMAX is REAL. Maximum value of q. 📈
NFAIL — NFAIL is INTEGER. Increment NFAIL by 1 each time the shift was too big. 💀
ITER — ITER is INTEGER. Increment ITER by 1 for each iteration.
NDIV — NDIV is INTEGER. Increment NDIV by 1 for each division. ➗
IEEE — IEEE is LOGICAL. Flag for IEEE or non IEEE arithmetic (passed to SLASQ5). ⚡
TTYPE — TTYPE is INTEGER. Shift type.
DMIN1 — DMIN1 is REAL
DMIN2 — DMIN2 is REAL
DN — DN is REAL
DN1 — DN1 is REAL
DN2 — DN2 is REAL
G — G is REAL
TAU — TAU is REAL

These are passed as arguments in order to save their values between calls to SLASQ3. 💾

Author

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

`subroutine slasq4`

subroutine slasq4 (integer I0, integer N0, real, dimension( * ) Z, integer PP, integer N0IN, real DMIN, real DMIN1, real DMIN2, real DN, real DN1, real DN2, real TAU, integer TTYPE, real G)

SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. 🧚

Purpose

SLASQ4 computes an approximation TAU to the smallest eigenvalue using values of d from the previous transform. 🔍

Parameters

I0 — I0 is INTEGER. First index. 1️⃣
N0 — N0 is INTEGER. Last index. 🔚
Z — Z is REAL array, dimension ( 4*N0 ). Z holds the qd array. 💾
PP — PP is INTEGER. PP=0 for ping, PP=1 for pong. 🏓
N0IN — N0IN is INTEGER. The value of N0 at start of EIGTEST.
DMIN — DMIN is REAL. Minimum value of d. 📉
DMIN1 — DMIN1 is REAL. Minimum value of d, excluding D( N0 ).
DMIN2
DMIN2 — REAL. The minimum value of d, but like… excluding D( N0 ) and D( N0-1 ) bestie, we're being picky here 💅
DN — REAL. It's d(N), the moment.
DN1 — REAL. That's d(N-1), one step back ✨
DN2 — REAL. And d(N-2), two steps back, we love a throwback 🎀
TAU — REAL. This is the shift.
TTYPE — INTEGER. Shift type, iykyk.
G — REAL. So G gets passed in as an argument purely so its value survives between calls to SLASQ4. She's carrying the emotional baggage for the group chat, no cap 🫶

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details

CNST1 = 9/16

💖 subroutine `slasq5`

subroutine slasq5 (integer I0, integer N0, real, dimension( * ) Z, integer
    PP, real TAU, real SIGMA, real DMIN, real DMIN1, real DMIN2, real DN,
    real DNM1, real DNM2, logical IEEE, real EPS)

SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. Okay so she's the working girl behind the scenes and we respect her 😭👑

Purpose

SLASQ5 computes one dqds transform in ping-pong form, one version for IEEE machines another for non IEEE machines. Two moods, both valid 🌈

Parameters

I0 — INTEGER. First index. The beginning, the origin story 💫
N0 — INTEGER. Last index. The finale 🎬
Z — REAL array, dimension ( 4*N ). Z holds the qd array. Also EMIN is stored in Z(4*N0) just to avoid an extra argument — efficient queen, we stan 🧚
PP — INTEGER. PP=0 for ping, PP=1 for pong. PING PONG 🏓 literally the whole vibe
TAU — REAL. This is the shift.
SIGMA — REAL. This is the accumulated shift up to this step. She's been through things, she's carrying it all 💜
DMIN — REAL. Minimum value of d.
DMIN1 — REAL. Minimum value of d, excluding D( N0 ).
DMIN2 — REAL. Minimum value of d, excluding D( N0 ) and D( N0-1 ).
DN — REAL. d(N0), the last value of d.
DNM1 — REAL. d(N0-1).
DNM2 — REAL. d(N0-2).
IEEE — LOGICAL. Flag for IEEE or non IEEE arithmetic. Pick your fighter ⚔️
EPS — REAL. This is the value of epsilon used.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 subroutine `slasq6`

subroutine slasq6 (integer I0, integer N0, real, dimension( * ) Z, integer
    PP, real DMIN, real DMIN1, real DMIN2, real DN, real DNM1, real DNM2)

SLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. Sis is back for round two 💅

Purpose

SLASQ6 computes one dqd (shift equal to zero) transform in ping-pong form, with protection against underflow and overflow. She's got a whole security detail, main character energy fr 🔥

Parameters

I0 — INTEGER. First index.
N0 — INTEGER. Last index.
Z — REAL array, dimension ( 4*N ). Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument — same trick, iconic ✨
PP — INTEGER. PP=0 for ping, PP=1 for pong. 🏓
DMIN — REAL. Minimum value of d.
DMIN1 — REAL. Minimum value of d, excluding D( N0 ).
DMIN2 — REAL. Minimum value of d, excluding D( N0 ) and D( N0-1 ).
DN — REAL. d(N0), the last value of d.
DNM1 — REAL. d(N0-1).
DNM2 — REAL. d(N0-2).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

⭐ subroutine `slasrt`

subroutine slasrt (character ID, integer N, real, dimension( * ) D, integer
    INFO)

SLASRT sorts numbers in increasing or decreasing order. Putting things in their place, respectfully 💅

Purpose

Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D').

Use Quick Sort, reverting to Insertion sort on arrays of size <= 20. Dimension of STACK limits N to about 2**32. Bestie that is a LOT of numbers, like genuinely unhinged amounts 📈

Parameters

ID — CHARACTER*1.
- = 'I': sort D in increasing order;
- = 'D': sort D in decreasing order.
N — INTEGER. The length of the array D.
D — REAL array, dimension (N). On entry, the array to be sorted. On exit, D has been sorted into increasing order (D(1) <= ... <= D(N)) or into decreasing order (D(1) >= ... >= D(N)), depending on ID. Glow up complete 💖
INFO — 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 `spttrf`

subroutine spttrf (integer N, real, dimension( * ) D, real, dimension( * )
    E, integer INFO)

SPTTRF — that's literally all she says, mysterious queen, we love a woman of few words 🙈

Purpose

SPTTRF computes the L*D*L**T factorization of a real symmetric positive definite tridiagonal matrix A. The factorization may also be regarded as having the form A = U**T*D*U. Same tea, different cup 🫖

Parameters

N — INTEGER. The order of the matrix A. N >= 0.
D — REAL array, dimension (N). On entry, the n diagonal elements of the tridiagonal matrix A. On exit, the n diagonal elements of the diagonal matrix D from the L*D*L**T factorization of A.
E — REAL array, dimension (N-1). On entry, the (n-1) subdiagonal elements of the tridiagonal matrix A. On exit, the (n-1) subdiagonal elements of the unit bidiagonal factor L from the L*D*L**T factorization of A. E can also be regarded as the superdiagonal of the unit bidiagonal factor U from the U**T*D*U factorization of A. Dual purpose, dual threat, we love to see it ✨
INFO — INTEGER.
- = 0: successful exit
- < 0: if INFO = -k, the k-th argument had an illegal value
- > 0: if INFO = k, the leading minor of order k is not positive definite; if k < N, the factorization could not be completed, while if k = N, the factorization was completed, but D(N) <= 0. Not me being dramatic but this is a whole situation 💀

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

👑 subroutine `sstebz`

subroutine sstebz (character RANGE, character ORDER, integer N, real VL,
    real VU, integer IL, integer IU, real ABSTOL, real, dimension( * ) D,
    real, dimension( * ) E, integer M, integer NSPLIT, real, dimension( * )
    W, integer, dimension( * ) IBLOCK, integer, dimension( * ) ISPLIT, real,
    dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

SSTEBZ — again, no notes, just vibes 😌

Purpose

SSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix T. The user may ask for all eigenvalues, all eigenvalues in the half-open interval (VL, VU], or the IL-th through IU-th eigenvalues. Choose your own adventure, bestie 🌈

To avoid overflow, the matrix must be scaled so that its largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest accuracy, it should not be much smaller than that. Okay so scale responsibly or she WILL overflow and it will be a mess 😭

See W. Kahan 'Accurate Eigenvalues of a Symmetric Tridiagonal Matrix', Report CS41, Computer Science Dept., Stanford University, July 21, 1966. Shoutout to the 1966 legend, real ones know 📚

Parameters

RANGE — CHARACTER*1.
- = 'A': ('All') all eigenvalues will be found.
- = 'V': ('Value') all eigenvalues in the half-open interval (VL, VU] will be found.
- = 'I': ('Index') the IL-th through IU-th eigenvalues (of the entire matrix) will be found.
ORDER — CHARACTER*1.
- = 'B': ('By Block') the eigenvalues will be grouped by split-off block (see IBLOCK, ISPLIT) and ordered from smallest to largest within the block.
- = 'E': ('Entire matrix') the eigenvalues for the entire matrix will be ordered from smallest to largest.
N — INTEGER. The order of the tridiagonal matrix T. N >= 0.
VL — REAL.

If RANGE='V', the lower bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.

VU — REAL.

If RANGE='V', the upper bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.

IL — INTEGER.

If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.

IU — INTEGER.

If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.

ABSTOL — REAL. The absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less. If ABSTOL is less than or equal to zero, then ULP*|T| will be used, where |T| means the 1-norm of T.

Eigenvalues will be computed most accurately when ABSTOL is set to twice the underflow threshold 2*SLAMCH('S'), not zero. Not zero bestie, please, she's begging 🥺

D
D — okay bestie, this one is a REAL array, dimension (N). It's literally the n diagonal elements of the tridiagonal matrix T. periodt. ✨
E — REAL array, dimension (N-1). The (n-1) off-diagonal elements of the tridiagonal matrix T. no cap. 🦋
M — INTEGER. The actual number of eigenvalues found. 0 <= M <= N. (See also the description of INFO=2,3.) iykyk. 👀
NSPLIT — INTEGER. The number of diagonal blocks in the matrix T. 1 <= NSPLIT <= N. slay. 💅
W — REAL array, dimension (N). On exit, the first M elements of W will contain the eigenvalues. (SSTEBZ may use the remaining N-M elements as workspace.) the way that works is lowkey iconic. 🫶
IBLOCK — INTEGER array, dimension (N). At each row/column j where E(j) is zero or small, the matrix T is considered to split into a block diagonal matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which block (from 1 to the number of blocks) the eigenvalue W(i) belongs. (SSTEBZ may use the remaining N-M elements as workspace.) fr fr. 🧩
ISPLIT — INTEGER array, dimension (N). The splitting points, at which T breaks up into submatrices. The first submatrix consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. (Only the first NSPLIT elements will actually be used, but since the user cannot know a priori what value NSPLIT will have, N words must be reserved for ISPLIT.) not me reading this twice 💀 🧠
WORK — REAL array, dimension (4*N). just vibes. ⚙️
IWORK — INTEGER array, dimension (3*N). also vibes. 🔧
INFO — INTEGER. okay so here's the tea:
- = 0: successful exit. we love to see it ✅
- < 0: if INFO = -i, the i-th argument had an illegal value. not cute ❌
- > 0: some or all of the eigenvalues failed to converge or were not computed: 💀
  - =1 or 3: Bisection failed to converge for some eigenvalues; these eigenvalues are flagged by a negative block number. The effect is that the eigenvalues may not be as accurate as the absolute and relative tolerances. This is generally caused by unexpectedly inaccurate arithmetic. oof 😭
  - =2 or 3: RANGE='I' only: Not all of the eigenvalues IL:IU were found.
    - Effect: M < IU+1-IL
    - Cause: non-monotonic arithmetic, causing the Sturm sequence to be non-monotonic.
    - Cure: recalculate, using RANGE='A', and pick out eigenvalues IL:IU. In some cases, increasing the PARAMETER 'FUDGE' may make things work. the cure is giving main character energy 🛠️
  - = 4: RANGE='I', and the Gershgorin interval initially used was too small. No eigenvalues were computed.
    - Probable cause: your machine has sloppy floating-point arithmetic.
    - Cure: Increase the PARAMETER 'FUDGE', recompile, and try again. literally just increase the fudge 🔁

Internal Parameters: ✨

RELFAC — REAL, default = 2.0e0. The relative tolerance. An interval (a,b] lies within 'relative tolerance' if b-a < RELFAC*ulp*max(|a|,|b|), where 'ulp' is the machine precision (distance from 1 to the next larger floating point number.) ngl this is a mouthful but we move 💅
FUDGE — REAL, default = 2. A 'fudge factor' to widen the Gershgorin intervals. Ideally, a value of 1 should work, but on machines with sloppy arithmetic, this needs to be larger. The default for publicly released versions should be large enough to handle the worst machine around. Note that this has no effect on accuracy of the solution. obsessed with a parameter named fudge 🍓

Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

   subroutine sstedc (character COMPZ, integer N, real, dimension( * ) D, real,
       dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real,
       dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK,
       integer LIWORK, integer INFO)

SSTEDC

Purpose:

SSTEDC computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method. The eigenvectors of a full or band real symmetric matrix can also be found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to tridiagonal form. 🌈

This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See SLAED3 for details. 🔥

Parameters

COMPZ — COMPZ is CHARACTER*1
- = 'N': Compute eigenvalues only.
- = 'I': Compute eigenvectors of tridiagonal matrix also.
- = 'V': Compute eigenvectors of original dense symmetric matrix also. On entry, Z contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.
N — N is INTEGER. The dimension of the symmetric tridiagonal matrix. N >= 0.
D — D is REAL array, dimension (N). On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. 📈
E — E is REAL array, dimension (N-1). On entry, the subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. 💥
Z — Z is REAL array, dimension (LDZ,N). On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced. 🦋
LDZ — LDZ is INTEGER. The leading dimension of the array Z. LDZ >= 1. If eigenvectors are desired, then LDZ >= max(1,N).
WORK — WORK is REAL array, dimension (MAX(1,LWORK)). On exit, if INFO = 0, WORK(1) returns the optimal LWORK. ⚙️
LWORK — LWORK is INTEGER. The dimension of the array WORK.
- If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
- If COMPZ = 'V' and N > 1 then LWORK must be at least ( 1 + 3*N + 2*N*lg N + 4*N**2 ), where lg( N ) = smallest integer k such that 2**k >= N.
- If COMPZ = 'I' and N > 1 then LWORK must be at least ( 1 + 4*N + N**2 ).
- Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LWORK need only be max(1,2*(N-1)).
- If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. 🔍
IWORK — IWORK is INTEGER array, dimension (MAX(1,LIWORK)). On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 🔧
LIWORK — LIWORK is INTEGER. The dimension of the array IWORK.
- If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
- If COMPZ = 'V' and N > 1 then LIWORK must be at least ( 6 + 6*N + 5*N*lg N ).
- If COMPZ = 'I' and N > 1 then LIWORK must be at least ( 3 + 5*N ).
- Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LIWORK need only be 1.
- If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the IWORK array, returns this value as the first entry of the IWORK array, and no error message related to LIWORK is issued by XERBLA. 🔎
INFO — INFO is INTEGER
- = 0: successful exit. ✅
- < 0: if INFO = -i, the i-th argument had an illegal value. ❌
- > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1). 😭

Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennessee

   subroutine ssteqr (character COMPZ, integer N, real, dimension( * ) D, real,
       dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real,
       dimension( * ) WORK, integer INFO)

SSTEQR

Purpose:

SSTEQR computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a full or band symmetric matrix can also be found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to tridiagonal form. 🌸

Parameters

COMPZ
COMPZ — CHARACTER*1, and bestie this is the whole personality of the call:
- = 'N': Compute eigenvalues only.
- = 'V': Compute eigenvalues and eigenvectors of the original symmetric matrix. On entry, Z must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form.
- = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.
N — INTEGER. The order of the matrix. N >= 0. that's it that's the parameter
D — REAL array, dimension (N). On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. ✨
E — REAL array, dimension (N-1). On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. gone. deleted. not coming back. 💀
Z — REAL array, dimension (LDZ, N). On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced. she simply does not get invited to the function call, and that's valid 🙈
LDZ — INTEGER. The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).
WORK — REAL array, dimension (max(1,2*N-2)). If COMPZ = 'N', then WORK is not referenced. scratch space that clocks out early when nobody needs it
INFO — INTEGER, the vibe report at the end:
- = 0: successful exit 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value (your fault, not hers)
- > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is orthogonally similar to the original matrix. it tried SO hard, bestie, it really did 😭

🎀 Author 🎀

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

🦋 subroutine ssterf 🦋

subroutine ssterf (integer N, real, dimension( * ) D, real, dimension( * ) E, integer INFO)

SSTERF 👑

Purpose:

SSTERF computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm. no eigenvectors here, she's strictly eigenvalues-only, and honestly? respect the focus 💅

Parameters

N — INTEGER. The order of the matrix. N >= 0.
D — REAL array, dimension (N). On entry, the n diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. ⭐
E — REAL array, dimension (N-1). On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. again. she's gone. we don't talk about it. 💀
INFO — INTEGER:
- = 0: successful exit ✅
- < 0: if INFO = -i, the i-th argument had an illegal value
- > 0: the algorithm failed to find all of the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero. 30*N tries and then she has to tell you the truth 😭

🎀 Author 🎀

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

💖 Author 💖

Generated automatically by Doxygen for LAPACK from the source code. (not a human, just a very diligent robot, and we thank her for her service) 🫶

and that's a wrap on the whole page, bestie!! 🌈 This was Version 3.11.0, dated Sun Nov 27 2022, and it lives in the manual as auxOTHERcomputational(3). go forth and tridiagonalize responsibly 💅✨🎀🦋🌸

📦 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

auxOTHERcomputational(3)

💅✨ auxOTHERcomputational(3) ✨💅

🌸 NAME 🌸

✨ SYNOPSIS ✨

🦋 Functions 🦋

🫶 Detailed Description 🫶

📚 Function Documentation 📚

🌟 character*1 function chla_transtype (integer TRANS) 🌟

🎯 Purpose:

🌸 Author 🌸

🦋 subroutine dbdsdc

💖 Purpose 💖

✨ Parameters ✨

🌈 Author 🌈

🫶 Contributors: 🫶

⭐ subroutine dbdsqr ⭐

💜 Purpose 💜

✨ Parameters ✨

🔧 Internal Parameters:

📝 Note:

✍️ Author

subroutine ddisna (character JOB, integer M, integer N, double precision, dimension( * ) D, double precision, dimension( * ) SEP, integer INFO)

🎯 Purpose:

📋 Parameters

✍️ Author

🎯 Purpose:

📋 Parameters

✍️ Author

🌟 Contributors:

subroutine dlaed1 (integer N, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

🎯 Purpose:

Parameters 💖

Author 🌟

Contributors: 🫶

✨ subroutine dlaed2 ✨

Purpose: 🎯

Parameters 📋

Author 🌟

Contributors: 🫶

✨ subroutine dlaed3 ✨

Purpose: 🎯

Parameters 📋

👑 Author

🫶 Contributors:

🌟 subroutine dlaed4 — the single-root finder ✨

🎯 Purpose:

📋 Parameters

🔍 Internal Parameters:

👑 Author

🫶 Contributors:

🌈 subroutine dlaed5 — the 2-by-2 girlie 💖

🎯 Purpose:

📋 Parameters

👑 Author

🫶 Contributors:

🔥 subroutine dlaed6 — one Newton step, that's the whole job 💀

🎯 Purpose:

📋 Parameters

👑 Author

📜 Further Details:

🫶 Contributors:

👑 subroutine dlaed7 — the updated eigensystem girlie 🦋

🎯 Purpose:

🎀 Parameters 🎀

👑 Author 👑

💖 Contributors 💖

🧚 subroutine dlaed8 🧚

🎀 Purpose 🎀

💅 Parameters 💅

Q2 🫶

Author 💖

Contributors 🌈

🦋 subroutine dlaed9 🦋

Purpose 🎯

Parameters 🌸

Author 💜

Contributors 🌈

🧚 subroutine dlaeda 🧚

Purpose 🎯

Parameters 🌸

🦋 subroutine `dbdsdc`

⭐ subroutine `dbdsqr` ⭐

🌟 subroutine `dlaed4` — the single-root finder ✨

🌈 subroutine `dlaed5` — the 2-by-2 girlie 💖

🔥 subroutine `dlaed6` — one Newton step, that's the whole job 💀

👑 subroutine `dlaed7` — the updated eigensystem girlie 🦋

`Q2` 🫶

🦋 subroutine `dlaed9` 🦋

🧚 subroutine `dlaeda` 🧚

💖 subroutine `dlagtf` 💖

🦋 `subroutine dlasq3` 🦋

🧚 `subroutine dlasq4` 🧚

🌈 `subroutine dlasq5` 🌈

🍓 `subroutine dlasq6` 🍓

👑 `subroutine dlasrt` 👑

💜 `subroutine dstebz` 💜

`dsteqr`

`dsterf`

`iladiag`

`ilaprec`

`ilatrans`

`ilauplo`

`sbdsdc`