OTHERauxiliary(3)

📦 version 3.10.0, in ubuntu jammy

🌸 NAME 🌸

okay so bestie, this is OTHERauxiliary(3) from LAPACK, and it's literally the page for Other Auxiliary Routines. no cap, that's the name. it's giving "the extras that still ate." 💅✨

✨ SYNOPSIS ✨

🎀 Modules 🎀

double
real
complex
complex16

✨ Functions ✨

okay bestie, buckle UP 💖 this is the full roster of little helper routines — the supporting cast, the background dancers, the ones doing the actual math while the main characters get all the credit. every single one of them, no skipping, periodt.

subroutine dcombssq (V1, V2) — DCOMBSSQ is out here adding two scaled sum of squares quantities together. that's the job. it does the job. slay 💅
logical function disnan (DIN) — DISNAN runs a vibe check on its input for NaN. is it a number? is it a NaN? it will let you know 🧐
subroutine dlabad (SMALL, LARGE) — DLABAD. and that's it. that's the whole description. no notes, iconic minimalism 👑
subroutine dlacpy (UPLO, M, N, A, LDA, B, LDB) — DLACPY copies all of a two-dimensional array, or just part of it, into another one. move it, copy it, love it 💖
subroutine dlae2 (A, B, C, RT1, RT2) — DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix. tiny matrix, big brain energy ✨
subroutine dlaebz (IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO) — DLAEBZ computes how many eigenvalues of a real symmetric tridiagonal matrix are less than or equal to a given value, and handles the other tasks that the routine sstebz needs done. it's the whole support system fr fr 🫶
subroutine dlaev2 (A, B, C, RT1, RT2, CS1, SN1) — DLAEV2 computes the eigenvalues and the eigenvectors of a 2-by-2 symmetric/Hermitian matrix. double the work, double the slay 💅
subroutine dlagts (JOB, N, A, B, C, D, IN, Y, TOL, INFO) — DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y, where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf. team effort, bestie 🤝
logical function dlaisnan (DIN1, DIN2) — DLAISNAN tests its input for NaN by comparing two arguments for inequality. sneaky little NaN detector 🕵️
integer function dlaneg (N, D, LLD, SIGMA, PIVMIN, R) — DLANEG computes the Sturm count. just. the Sturm count. 👑
double precision function dlanst (NORM, N, D, E) — DLANST hands you the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value, of a real symmetric tridiagonal matrix. pick your fighter ✨
double precision function dlapy2 (X, Y) — DLAPY2 returns sqrt(x2+y2). pythagoras said what he said 📐
double precision function dlapy3 (X, Y, Z) — DLAPY3 returns sqrt(x2+y2+z2). pythagoras but make it three-dimensional 🦋
subroutine dlarnv (IDIST, ISEED, N, X) — DLARNV returns a vector of random numbers from a uniform or normal distribution. chaos, but make it mathematically rigorous 🎲
subroutine dlarra (N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO) — DLARRA computes the splitting points with the specified threshold. chopping it up, respectfully ✂️
subroutine dlarrb (N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO) — DLARRB provides limited bisection to locate eigenvalues for more accuracy. it's giving precision 💅
subroutine dlarrc (JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO) — DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix. counting, but like, mathematically 🔢
subroutine dlarrd (RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO) — DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. "suitable" is doing a lot of work in that sentence and honestly? relatable 😭
subroutine dlarre (RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO) — DLARRE takes the tridiagonal matrix T, sets the small off-diagonal elements to zero, and then for each unreduced block Ti it finds base representations and eigenvalues. multi-step icon behavior 👑
subroutine dlarrf (N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO) — DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. giving that one eigenvalue its own personal space 🫶
subroutine dlarrj (N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO) — DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T. glow-up arc for your eigenvalues ✨
subroutine dlarrk (N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO) — DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. just one. it's focused. it's got main character energy 💅
subroutine dlarrr (N, D, E, INFO) — DLARRR performs tests to decide whether the symmetric tridiagonal matrix T is worth the expensive computations that guarantee high relative accuracy in the eigenvalues. it's the friend who asks "is this really worth it" before you spend money 🙈

subroutine dlartgp (F, G, CS, SN, R) — DLARTGP generates a plane rotation so that the diagonal is nonnegative. rotate, bestie 💫
subroutine dlaruv (ISEED, N, X) — DLARUV returns a vector of n random real numbers from a uniform distribution. more randomness, now with a lowercase n 🎀
subroutine dlas2 (F, G, H, SSMIN, SSMAX) — DLAS2 computes singular values of a 2-by-2 triangular matrix. tiny triangle, huge significance ✨
subroutine dlascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO) — DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. math said "scale it" and it scaled it 💅
subroutine dlasd0 (N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, WORK, INFO) — DLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. used by sbdsdc. 🦋
subroutine dlasd1 (NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, IDXQ, IWORK, WORK, INFO) — DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. used by sbdsdc. 👑
subroutine dlasd2 (NL, NR, SQRE, K, D, Z, ALPHA, BETA, U, LDU, VT, LDVT, DSIGMA, U2, LDU2, VT2, LDVT2, IDXP, IDX, IDXC, IDXQ, COLTYP, INFO) — DLASD2 merges the two sets of singular values together into a single sorted set. used by sbdsdc. friendship bracelet behavior 🤝
subroutine dlasd3 (NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, INFO) — DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. used by sbdsdc. that's a lot of steps and it ate every single one 🔥
subroutine dlasd4 (N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO) — DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. used by dbdsdc. 💜
subroutine dlasd5 (I, D, Z, DELTA, RHO, DSIGMA, WORK) — DLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. used by sbdsdc. the little sibling of DLASD4 and it's thriving 🌸
subroutine dlasd6 (ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, IWORK, INFO) — DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. used by sbdsdc. 💅
subroutine dlasd7 (ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, INFO) — DLASD7 merges the two sets of singular values together into a single sorted set. then it tries to deflate the size of the problem. used by sbdsdc. ✨
subroutine dlasd8 (ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, DSIGMA, WORK, INFO) — DLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. used by sbdsdc. 🎀
subroutine dlasda (ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO) — DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. used by sbdsdc. the full SVD, the whole package, no half measures 👑
subroutine dlasdq (UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO) — DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. used by sbdsdc. 🦋
subroutine dlasdt (N, LVL, ND, INODE, NDIML, NDIMR, MSUB) — DLASDT creates a tree of subproblems for bidiagonal divide and conquer. used by sbdsdc. it's giving organizational skills 🌳
subroutine dlaset (UPLO, M, N, ALPHA, BETA, A, LDA) — DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values. setting the stage, literally 💅
subroutine dlasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA) — DLASR applies a sequence of plane rotations to a general rectangular matrix. spin it, bestie 💫
subroutine dlasv2 (F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL) — DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix. small triangle, full SVD, obsessed ✨
integer function ieeeck (ISPEC, ZERO, ONE) — IEEECK. that's all we got. mysterious queen, no further comment 👑
integer function iladlc (M, N, A, LDA) — ILADLC scans a matrix for its last non-zero column. finding the end of the line, respectfully 🔍
integer function iladlr (M, N, A, LDA) — ILADLR scans a matrix for its last non-zero row. same energy, other axis 🫶
integer function ilaenv (ISPEC, NAME, OPTS, N1, N2, N3, N4) — ILAENV. no description, just vibes and seven arguments 💜
integer function ilaenv2stage (ISPEC, NAME, OPTS, N1, N2, N3, N4) — ILAENV2STAGE. the sequel. same signature, new era 🎬
integer function iparmq (ISPEC, NAME, OPTS, N, ILO, IHI, LWORK) — IPARMQ. no notes, only mystery and parameters 💅
logical function lsamen (N, CA, CB) — LSAMEN. that's the whole entry. minimalism is a lifestyle ✨
subroutine scombssq (V1, V2) — SCOMBSSQ adds two scaled sum of squares quantities. same job as its double-precision twin, different letter 🎀
logical function sisnan (SIN) — SISNAN tests input for NaN. 🧐
subroutine slabad (SMALL, LARGE) — SLABAD. nothing more to say and honestly, respect 👑
subroutine slacpy (UPLO, M, N, A, LDA, B, LDB) — SLACPY copies all or part of one two-dimensional array to another. 💖
subroutine slae2 (A, B, C, RT1, RT2) — SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix. 🌸
subroutine slaebz (IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO) — SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs the other tasks required by the routine sstebz. 🔥
subroutine slaev2 (A, B, C, RT1, RT2, CS1, SN1) — SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. 💅
subroutine slag2d (M, N, SA, LDSA, A, LDA, INFO) — SLAG2D converts a single precision matrix to a double precision matrix. glow-up, but numerically 🦋
subroutine slagts (JOB, N, A, B, C, D, IN, Y, TOL, INFO) — SLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y, where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf. 🫶
logical function slaisnan (SIN1, SIN2) — SLAISNAN tests input for NaN by comparing two arguments for inequality. 🕵️
integer function slaneg (N, D, LLD, SIGMA, PIVMIN, R) — SLANEG computes the Sturm count. 👑
real function slanst (NORM, N, D, E) — SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix. ✨
real function slapy2 (X, Y) — SLAPY2 returns sqrt(x2+y2). 📐
real function slapy3 (X, Y, Z) — SLAPY3 returns sqrt(x2+y2+z2). 💫
subroutine slarnv (IDIST, ISEED, N, X) — SLARNV returns a vector of random numbers from a uniform or normal distribution. 🎲
subroutine slarra (N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO) — SLARRA computes the splitting points with the specified threshold. ✂️
subroutine slarrb (N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO) — SLARRB
SLARRB — provides limited bisection to locate eigenvalues for more accuracy. 🎀 no cap, she's just out here refining eigenvalues.
slarrc — subroutine slarrc (JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO) — SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix. 💖
slarrd — subroutine slarrd (RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO) — SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. ✨
slarre — subroutine slarre (RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO) — SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues. 💅
slarrf — subroutine slarrf (N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO) — SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. 🦋
slarrj — subroutine slarrj (N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO) — SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T. 🌸
slarrk — subroutine slarrk (N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO) — SLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. 👑
slarrr — subroutine slarrr (N, D, E, INFO) — SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues. 🔥
slartgp — subroutine slartgp (F, G, CS, SN, R) — SLARTGP generates a plane rotation so that the diagonal is nonnegative. 💫
slaruv — subroutine slaruv (ISEED, N, X) — SLARUV returns a vector of n random real numbers from a uniform distribution. 🎲
slas2 — subroutine slas2 (F, G, H, SSMIN, SSMAX) — SLAS2 computes singular values of a 2-by-2 triangular matrix. 🧚
slascl — subroutine slascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO) — SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. 💜
slasd0 — subroutine slasd0 (N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, WORK, INFO) — SLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc. 🌈
slasd1 — subroutine slasd1 (NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, IDXQ, IWORK, WORK, INFO) — SLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc. ⭐
slasd2 — subroutine slasd2 (NL, NR, SQRE, K, D, Z, ALPHA, BETA, U, LDU, VT, LDVT, DSIGMA, U2, LDU2, VT2, LDVT2, IDXP, IDX, IDXC, IDXQ, COLTYP, INFO) — SLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc. 💖
slasd3 — subroutine slasd3 (NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, INFO) — SLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc. 🎀
slasd4 — subroutine slasd4 (N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO) — SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by sbdsdc. ✨
slasd5 — subroutine slasd5 (I, D, Z, DELTA, RHO, DSIGMA, WORK) — SLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. Used by sbdsdc. 🦋
slasd6 — subroutine slasd6 (ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, IWORK, INFO) — SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc. 🌸
slasd7 — subroutine slasd7 (ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, INFO) — SLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc. 💅
slasd8 — subroutine slasd8 (ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, DSIGMA, WORK, INFO) — SLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc. 👑
slasda — subroutine slasda (ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO) — SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. 🔥
slasdq — subroutine slasdq (UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO) — SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. 💫
slasdt — subroutine slasdt (N, LVL, ND, INODE, NDIML, NDIMR, MSUB) — SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc. 🧚
slaset — subroutine slaset (UPLO, M, N, ALPHA, BETA, A, LDA) — SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values. 💜
slasr — subroutine slasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA) — SLASR applies a sequence of plane rotations to a general rectangular matrix. 🌈
slasv2 — subroutine slasv2 (F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL) — SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix. ⭐
xerbla — subroutine xerbla (SRNAME, INFO) — XERBLA. 💀
xerbla_array — subroutine xerbla_array (SRNAME_ARRAY, SRNAME_LEN, INFO) — XERBLA_ARRAY. 🙈

✨ Detailed Description ✨

This is the group of Other Auxiliary routines. 💅

🌸 Function Documentation 🌸

💖 `subroutine dcombssq (double precision, dimension( 2 ) V1, double precision, dimension( 2 ) V2)`

DCOMBSSQ adds two scaled sum of squares quantities. 🎀

Purpose:

DCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2. That is,

   V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq
                           + V2_scale**2 * V2_sumsq

Parameters

V1 — V1 is DOUBLE PRECISION array, dimension (2). The first scaled sum. V1(1) = V1_scale, V1(2) = V1_sumsq. ✨
V2 — V2 is DOUBLE PRECISION array, dimension (2). The second scaled sum. V2(1) = V2_scale, V2(2) = V2_sumsq. 💖

Author

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

🦋 `logical function disnan (double precision, intent(in) DIN)`

DISNAN tests input for NaN. 👑

Purpose:

DISNAN returns .TRUE. if its argument is NaN, and .FALSE. otherwise. To be replaced by the Fortran 2003 intrinsic in the future. 🔮

Parameters

DIN — DIN is DOUBLE PRECISION. Input to test for NaN. 🧚

Author

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

🌈 `subroutine dlabad (double precision SMALL, double precision LARGE)`

DLABAD. 💅

Purpose:

DLABAD takes as input the values computed by DLAMCH for underflow and overflow, and returns the square root of each of these values if the log of LARGE is sufficiently large. This subroutine is intended to identify machines with a large exponent range, such as the Crays, and redefine the underflow and overflow limits to be the square roots of the values computed by DLAMCH. This subroutine is needed because DLAMCH does not compensate for poor arithmetic in the upper half of the exponent range, as is found on a Cray. 😭

Parameters

SMALL — SMALL is DOUBLE PRECISION. On entry, the underflow threshold as computed by DLAMCH. On exit, if LOG10(LARGE) is sufficiently large, the square root of SMALL, otherwise unchanged. 🌸
LARGE — LARGE is DOUBLE PRECISION. On entry, the overflow threshold as computed by DLAMCH. On exit, if LOG10(LARGE) is sufficiently large, the square root of LARGE, otherwise unchanged. 💜

Author

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

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

DLACPY copies all or part of one two-dimensional array to another. 🫶

Purpose:

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

Parameters

UPLO — UPLO is CHARACTER*1. Specifies the part of the matrix A to be copied to B. = 'U': Upper triangular part. = 'L': Lower triangular part. Otherwise: All of the matrix A. 💖
M — M is INTEGER. The number of rows of the matrix A. M >= 0. 🎀
N — N is INTEGER. The number of columns of the matrix A. N >= 0. 🌸
A — A is DOUBLE PRECISION array, dimension (LDA,N). The m by n matrix A. If UPLO = 'U', only the upper triangle or trapezoid is accessed; if UPLO = 'L', only the lower triangle or trapezoid is accessed. 🦋
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M). 💅
B — B is DOUBLE PRECISION array, dimension (LDB,N). On exit, B = A in the locations specified by UPLO. 👑
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M). 🔥

Author

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

💜 `subroutine dlae2 (double precision A, double precision B, double precision C, double precision RT1, double precision RT2)`

DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix. 💫

Purpose:

DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix

   [  A   B  ]
   [  B   C  ].

On return, RT1 is the eigenvalue of larger absolute value, and RT2 is the eigenvalue of smaller absolute value. 🌈

Parameters

A — A is DOUBLE PRECISION. The (1,1) element of the 2-by-2 matrix. 🍓
B — B is DOUBLE PRECISION. it's the (1,2) and (2,1) elements of the 2-by-2 matrix. 💅
C — C is DOUBLE PRECISION. it's the (2,2) element of the 2-by-2 matrix. 🎀
RT1 — RT1 is DOUBLE PRECISION. it's the eigenvalue of larger absolute value. 👑
RT2 — RT2 is DOUBLE PRECISION. it's the eigenvalue of smaller absolute value. ✨

Author 🫶

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details: 💖

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

subroutine dlaebz

   subroutine dlaebz (integer IJOB, integer NITMAX, integer N, integer MMAX,
       integer MINP, integer NBMIN, double precision ABSTOL, double precision
       RELTOL, double precision PIVMIN, double precision, dimension( * ) D,
       double precision, dimension( * ) E, double precision, dimension( * ) E2,
       integer, dimension( * ) NVAL, double precision, dimension( mmax, * ) AB,
       double precision, dimension( * ) C, integer MOUT, integer, dimension(
       mmax, * ) NAB, double precision, dimension( * ) WORK, integer,
       dimension( * ) IWORK, integer INFO)

DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz. no cap, it's giving very specific eigenvalue energy. 🦋

Purpose: 🎯

DLAEBZ contains the iteration loops which compute and use the function N(w), which is the count of eigenvalues of a symmetric tridiagonal matrix T less than or equal to its argument w. It performs a choice of two types of loops:

IJOB=1, followed by IJOB=2: It takes as input a list of intervals and returns a list of sufficiently small intervals whose union contains the same eigenvalues as the union of the original intervals. The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. The output interval (AB(j,1),AB(j,2)] will contain eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. 🌸

IJOB=3: It performs a binary search in each input interval (AB(j,1),AB(j,2)] for a point w(j) such that N(w(j))=NVAL(j), and uses C(j) as the starting point of the search. If such a w(j) is found, then on output AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output (AB(j,1),AB(j,2)] will be a small interval containing the point where N(w) jumps through NVAL(j), unless that point lies outside the initial interval. 🔥

Note that the intervals are in all cases half-open intervals, i.e., of the form (a,b], which includes b but not a. 🎀

To avoid underflow, the matrix should be scaled so that its largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value. To assure the most accurate computation of small eigenvalues, the matrix should be scaled to be not much smaller than that, either. ✨

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966 📚

Note: the arguments are, in general, not checked for unreasonable values. 💀

Parameters

IJOB — IJOB is INTEGER. Specifies what is to be done:
- = 1: Compute NAB for the initial intervals.
- = 2: Perform bisection iteration to find eigenvalues of T.
- = 3: Perform bisection iteration to invert N(w), i.e., to find a point which has a specified number of eigenvalues of T to its left.
- Other values will cause DLAEBZ to return with INFO=-1.
NITMAX — NITMAX is INTEGER. The maximum number of "levels" of bisection to be performed, i.e., an interval of width W will not be made smaller than 2^(-NITMAX) * W. If not all intervals have converged after NITMAX iterations, then INFO is set to the number of non-converged intervals.
N — N is INTEGER. The dimension n of the tridiagonal matrix T. It must be at least 1.
MMAX — MMAX is INTEGER. The maximum number of intervals. If more than MMAX intervals are generated, then DLAEBZ will quit with INFO=MMAX+1.
MINP — MINP is INTEGER. The initial number of intervals. It may not be greater than MMAX.
NBMIN — NBMIN is INTEGER. The smallest number of intervals that should be processed using a vector loop. If zero, then only the scalar loop will be used.
ABSTOL — ABSTOL is DOUBLE PRECISION. The minimum (absolute) width of an interval. When an interval is narrower than ABSTOL, or than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged. This must be at least zero.
RELTOL — RELTOL is DOUBLE PRECISION. The minimum relative width of an interval. When an interval is narrower than ABSTOL, or than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged. Note: this should always be at least radix*machine epsilon.
PIVMIN — PIVMIN is DOUBLE PRECISION. The minimum absolute value of a "pivot" in the Sturm sequence loop. This must be at least max |e(j)**2|*safe_min and at least safe_min, where safe_min is at least the smallest number that can divide one without overflow.
D — D is DOUBLE PRECISION array, dimension (N). The diagonal elements of the tridiagonal matrix T.
E — E is DOUBLE PRECISION array, dimension (N). The offdiagonal elements of the tridiagonal matrix T in positions 1 through N-1. E(N) is arbitrary.
E2 — E2 is DOUBLE PRECISION array, dimension (N). The squares of the offdiagonal elements of the tridiagonal matrix T. E2(N) is ignored.
NVAL — NVAL is INTEGER array, dimension (MINP). If IJOB=1 or 2, not referenced. If IJOB=3, the desired values of N(w). The elements of NVAL will be reordered to correspond with the intervals in AB. Thus, NVAL(j) on output will not, in general be the same as NVAL(j) on input, but it will correspond with the interval (AB(j,1),AB(j,2)] on output.
AB — AB is DOUBLE PRECISION array, dimension (MMAX,2). The endpoints of the intervals. AB(j,1) is a(j), the left endpoint of the j-th interval, and AB(j,2) is b(j), the right endpoint of the j-th interval. The input intervals will, in general, be modified, split, and reordered by the calculation.
C — C is DOUBLE PRECISION array, dimension (MMAX). If IJOB=1, ignored. If IJOB=2, workspace. If IJOB=3, then on input C(j) should be initialized to the first search point in the binary search.
MOUT — MOUT is INTEGER. If IJOB=1, the number of eigenvalues in the intervals. If IJOB=2 or 3, the number of intervals output. If IJOB=3, MOUT will equal MINP.
NAB — NAB is INTEGER array, dimension (MMAX,2). If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). If IJOB=2, then on input, NAB(i,j) should be set. It must satisfy the condition: N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), which means that in interval i only eigenvalues NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with IJOB=1. On output, NAB(i,j) will contain max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of the input interval that the output interval (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the the input values of NAB(k,1) and NAB(k,2). If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), unless N(w) > NVAL(i) for all search points w, in which case NAB(i,1) will not be modified, i.e., the output value will be the same as the input value (modulo reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) for all search points w, in which case NAB(i,2) will not be modified. Normally, NAB should be set to some distinctive value(s) before DLAEBZ is called.
WORK — WORK is DOUBLE PRECISION array, dimension (MMAX). Workspace.
IWORK — IWORK is INTEGER array, dimension (MMAX). Workspace.
INFO — INFO is INTEGER.
- = 0: All intervals converged.
- = 1--MMAX: The last INFO intervals did not converge.
- = MMAX+1: More than MMAX intervals were generated.

Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details: 🌈

This routine is intended to be called only by other LAPACK routines, thus the interface is less user-friendly. It is intended for two purposes:

(a) finding eigenvalues. In this case, DLAEBZ should have one or more initial intervals set up in AB, and DLAEBZ should be called with IJOB=1. This sets up NAB, and also counts the eigenvalues. Intervals with no eigenvalues would usually be thrown out at this point. Also, if not all the eigenvalues in an interval i are desired, NAB(i,1) can be increased or NAB(i,2) decreased. For example, set NAB(i,1)=NAB(i,2)-1 to get the largest eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX no smaller than the value of MOUT returned by the call with IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the tolerance specified by ABSTOL and RELTOL. 💖
(b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). In this case, start with a Gershgorin interval (a,b). Set up AB to contain 2 search intervals, both initially (a,b). One NVAL element should contain f-1 and the other should contain l, while C should contain a and b, resp. NAB(i,1) should be -1 and NAB(i,2) should be N+1, to flag an error if the desired interval does not lie in (a,b). DLAEBZ is then called with IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and w(l-r)=...=w(l+k) are handled similarly. 🦋

💅 `subroutine dlaev2`

subroutine dlaev2 (double precision A, double precision B, double precision
    C, double precision RT1, double precision RT2, double precision CS1,
    double precision SN1)

DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. bestie this one is TINY and it's giving little matrix energy ✨

Purpose:

DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix

   [  A   B  ]
   [  B   C  ].

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

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

Parameters

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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

RT1 is accurate to a few ulps barring over/underflow.

RT2 may be inaccurate if there is massive cancellation in the determinant A*C-B*B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases. 😭

CS1 and SN1 are accurate to a few ulps barring over/underflow.

Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds

   underflow_threshold / macheps.

🌈 `subroutine dlagts`

subroutine dlagts (integer JOB, integer N, double precision, dimension( * )
    A, double precision, dimension( * ) B, double precision, dimension( * )
    C, double precision, dimension( * ) D, integer, dimension( * ) IN,
    double precision, dimension( * ) Y, double precision TOL, integer INFO)

DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y, where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf. slay, it's the follow-up act 🎤

Purpose:

DLAGTS may be used to solve one of the systems of equations

   (T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,

where T is an n by n tridiagonal matrix, for x, following the factorization of (T - lambda*I) as

   (T - lambda*I) = P*L*U ,

by routine DLAGTF. The choice of equation to be solved is controlled by the argument JOB, and in each case there is an option to perturb zero or very small diagonal elements of U, this option being intended for use in applications such as inverse iteration. 💅

Parameters

JOB JOB is INTEGER Specifies the job to be performed by DLAGTS as follows:
- = 1: The equations (T - lambda*I)x = y are to be solved, but diagonal elements of U are not to be perturbed.
- = -1: The equations (T - lambda*I)x = y are to be solved and, if overflow would otherwise occur, the diagonal elements of U are to be perturbed. See argument TOL below.
- = 2: The equations (T - lambda*I)**Tx = y are to be solved, but diagonal elements of U are not to be perturbed.
- = -2: The equations (T - lambda*I)**Tx = y are to be solved and, if overflow would otherwise occur, the diagonal elements of U are to be perturbed. See argument TOL below. 🦋
N N is INTEGER The order of the matrix T.
A A is DOUBLE PRECISION array, dimension (N) On entry, A must contain the diagonal elements of U as returned from DLAGTF.
B B is DOUBLE PRECISION array, dimension (N-1) On entry, B must contain the first super-diagonal elements of U as returned from DLAGTF.
C C is DOUBLE PRECISION array, dimension (N-1) On entry, C must contain the sub-diagonal elements of L as returned from DLAGTF.
D D is DOUBLE PRECISION array, dimension (N-2) On entry, D must contain the second super-diagonal elements of U as returned from DLAGTF.
IN IN is INTEGER array, dimension (N) On entry, IN must contain details of the matrix P as returned from DLAGTF. 🎀
Y Y is DOUBLE PRECISION array, dimension (N) On entry, the right hand side vector y. On exit, Y is overwritten by the solution vector x. periodt.
TOL TOL is DOUBLE PRECISION On entry, with JOB < 0, TOL should be the minimum perturbation to be made to very small diagonal elements of U. TOL should normally be chosen as about eps*norm(U), where eps is the relative machine precision, but if TOL is supplied as non-positive, then it is reset to eps*max( abs( u(i,j) ) ). If JOB > 0 then TOL is not referenced.

On exit, TOL is changed as described above, only if TOL is non-positive on entry. Otherwise TOL is unchanged.

INFO INFO is INTEGER
- = 0: successful exit
- < 0: if INFO = -i, the i-th argument had an illegal value
- > 0: overflow would occur when computing the INFO(th) element of the solution vector x. This can only occur when JOB is supplied as positive and either means that a diagonal element of U is very small, or that the elements of the right-hand side vector y are very large. 😭

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ `logical function dlaisnan`

logical function dlaisnan (double precision, intent(in) DIN1, double
    precision, intent(in) DIN2)

DLAISNAN tests input for NaN by comparing two arguments for inequality. not the glow-up anyone asked for but here we are 🙈

Purpose:

This routine is not for general use. It exists solely to avoid over-optimization in DISNAN.

DLAISNAN checks for NaNs by comparing its two arguments for inequality. NaN is the only floating-point value where NaN != NaN returns .TRUE.. To check for NaNs, pass the same variable as both arguments.

A compiler must assume that the two arguments are not the same variable, and the test will not be optimized away. Interprocedural or whole-program optimization may delete this test. The ISNAN functions will be replaced by the correct Fortran 03 intrinsic once the intrinsic is widely available. 🔥

Parameters

DIN1 DIN1 is DOUBLE PRECISION
DIN2 DIN2 is DOUBLE PRECISION Two numbers to compare for inequality. 💜

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 `integer function dlaneg`

integer function dlaneg (integer N, double precision, dimension( * ) D,
    double precision, dimension( * ) LLD, double precision SIGMA, double
    precision PIVMIN, integer R)

DLANEG computes the Sturm count. iconic name, iconic job 👑

Purpose:

DLANEG computes the Sturm count, the number of negative pivots encountered while factoring tridiagonal T - sigma I = L D L^T. This implementation works directly on the factors without forming the tridiagonal matrix T. The Sturm count is also the number of eigenvalues of T less than sigma.

This routine is called from DLARRB.

The current routine does not use the PIVMIN parameter but rather requires IEEE-754 propagation of Infinities and NaNs. This routine also has no input range restrictions but does require default exception handling such that x/0 produces Inf when x is non-zero, and Inf/Inf produces NaN. For more information, see:

  Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
  Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
  Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
  (Tech report version in LAWN 172 with the same title.)

Parameters

N

okay bestie we're back in it, mid-parameter-list and thriving 💅 here's the rest of the args for the function we were just doing:

N — INTEGER. The order of the matrix. slay, that's it, that's the whole job.
D — DOUBLE PRECISION array, dimension (N). The N diagonal elements of the diagonal matrix D. the main characters of the diagonal, honestly.
LLD — DOUBLE PRECISION array, dimension (N-1). The (N-1) elements L(i)*L(i)*D(i). not L(i)*D(i), not L(i), it's the whole squared-and-scaled situation, no cap.
SIGMA — DOUBLE PRECISION. Shift amount in T - sigma I = L D L^T. the lil nudge that shifts everything, iconic.
PIVMIN — DOUBLE PRECISION. The minimum pivot in the Sturm sequence. May be used when zero pivots are encountered on non-IEEE-754 architectures. 🫶
R — INTEGER. The twist index for the twisted factorization that is used for the negcount. twisted factorization?? yes and it ate, fr fr.

🌸 Author 🌸

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

🎀 Contributors 🎀

Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA
Jason Riedy, University of California, Berkeley, USA

👑 `dlanst` — the norm queen

double precision function dlanst (character NORM, integer N, double
    precision, dimension( * ) D, double precision, dimension( * ) E)

DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.

Purpose:

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

Returns:

DLANST gives you:

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

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

Parameters:

NORM — CHARACTER*1. Specifies the value to be returned in DLANST as described above. pick your fighter, bestie.
N — INTEGER. The order of the matrix A. N >= 0. When N = 0, DLANST is set to zero. ✨
D — DOUBLE PRECISION array, dimension (N). The diagonal elements of A.
E — DOUBLE PRECISION array, dimension (N-1). The (n-1) sub-diagonal or super-diagonal elements of A.

Author:

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

💜 `dlapy2` — pythagoras but make it safe

double precision function dlapy2 (double precision X, double precision Y)

DLAPY2 returns sqrt(x2+y2). yes, the whole entire point is the hypotenuse. 💅

Purpose:

DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary overflow and unnecessary underflow. it's the cautious queen of the diagonal — no drama, no overflow, no underflow, just vibes and accuracy. 🌈

Parameters:

X — DOUBLE PRECISION.
Y — DOUBLE PRECISION. X and Y specify the values x and y.

Author:

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

🦋 `dlapy3` — pythagoras but in 3D, main character energy

double precision function dlapy3 (double precision X, double precision Y,
    double precision Z)

DLAPY3 returns sqrt(x2+y2+z2). three coordinates, one iconic diagonal. 🔥

Purpose:

DLAPY3 returns sqrt(x**2+y**2+z**2), taking care not to cause unnecessary overflow and unnecessary underflow. same energy as dlapy2 but she brought a friend. 🫶

Parameters:

X — DOUBLE PRECISION.
Y — DOUBLE PRECISION.
Z — DOUBLE PRECISION. X, Y and Z specify the values x, y and z.

Author:

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

✨ `dlarnv` — the random number generator, no cap

subroutine dlarnv (integer IDIST, integer, dimension( 4 ) ISEED, integer N,
    double precision, dimension( * ) X)

DLARNV returns a vector of random numbers from a uniform or normal distribution. it's giving chaos, but make it deterministic. 🎲

Purpose:

DLARNV returns a vector of n random real numbers from a uniform or normal distribution.

Parameters:

IDIST — INTEGER. Specifies the distribution of the random numbers:
- = 1: uniform (0,1)
- = 2: uniform (-1,1)
- = 3: normal (0,1)

choose your fate, bestie. 🍓

ISEED — INTEGER array, dimension (4). On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. On exit, the seed is updated. 🫶 yes the seed changes on the way out, that's the whole point, iykyk.

bestie the rules here are real: elements 0–4095, and ISEED(4) has to be odd. don't skip that, it's not a suggestion. 💀

N — INTEGER. The number of random numbers to be generated.
X — DOUBLE PRECISION array, dimension (N). The generated random numbers.

Author:

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

Further Details:

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

🧚 `dlarra` — she finds the splitting points

subroutine dlarra (integer N, double precision, dimension( * ) D, double
    precision, dimension( * ) E, double precision, dimension( * ) E2, double
    precision SPLTOL, double precision TNRM, integer NSPLIT, integer,
    dimension( * ) ISPLIT, integer INFO)

DLARRA computes the splitting points with the specified threshold. 🌸

Purpose:

Compute the splitting points with threshold SPLTOL. DLARRA sets any "small" off-diagonal elements to zero. tiny numbers get yeeted to zero, respectfully. 🙈

Parameters:

N — INTEGER. The order of the matrix. N > 0. (note: strictly greater, no zero allowed here bestie)
D — DOUBLE PRECISION array, dimension (N). On entry, the N diagonal elements of the tridiagonal matrix T.
E — DOUBLE PRECISION array, dimension (N). On entry, the first (N-1) entries contain the subdiagonal elements of the tridiagonal matrix T; E(N) need not be set. On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, are set to zero, the other entries of E are untouched. the untouched ones stay untouched, we don't do unnecessary edits here. 💅
E2 — DOUBLE PRECISION array, dimension (N). On entry, the first (N-1) entries contain the SQUARES of the subdiagonal elements of the tridiagonal matrix T; E2(N) need not be set. On exit, the entries E2( ISPLIT( I ) ), 1 <= I <= NSPLIT, have been set to zero. 🔥 squared, not the raw values — read that twice.
SPLTOL — DOUBLE PRECISION. The threshold for splitting. Two criteria can be used:
- SPLTOL<0 : criterion based on absolute off-diagonal value
- SPLTOL>0 : criterion that preserves relative accuracy

the sign literally picks the strategy, it's giving "your vibe check is your algorithm." ⭐

TNRM — DOUBLE PRECISION. The norm of the matrix.
NSPLIT — INTEGER. The number of blocks T splits into. 1 <= NSPLIT <= N.
ISPLIT — INTEGER array, dimension (N). The splitting points, at which T breaks up into blocks. The first block consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. it's a whole seating chart for your matrix, iconic. 🎀
INFO — INTEGER.
- = 0: successful exit ✅

Author:

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

Contributors:

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

💖 `dlarrb` — limited bisection, maximum accuracy

subroutine dlarrb (integer N, double precision, dimension( * ) D, double
    precision, dimension( * ) LLD, integer IFIRST, integer ILAST, double
    precision RTOL1, double precision RTOL2, integer OFFSET, double
    precision, dimension( * ) W, double precision, dimension( * ) WGAP,
    double precision, dimension( * ) WERR, double precision, dimension( * )
    WORK, integer, dimension( * ) IWORK, double precision PIVMIN, double
    precision SPDIAM, integer TWIST, integer INFO)

DLARRB provides limited bisection to locate eigenvalues for more accuracy. 🌈 she's not doing the whole thing, she's just refining. limited bisection, unlimited slay.

Purpose:

Given the relatively robust representation (RRR) L D L^T, DLARRB does "limited" bisection to refine the eigenvalues of L D L^T, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial guesses for these eigenvalues are input in W, the corresponding estimate of the error in these guesses and their gaps are input in WERR and WGAP, respectively. During bisection, intervals [left, right] are maintained by storing their mid-points and semi-widths in the arrays W and WERR respectively. 🫶 mid-points and semi-widths, that's the whole storage scheme, no cap.

Parameters:

N — INTEGER. The order of the matrix.
D — DOUBLE PRECISION array, dimension (N). The N diagonal elements of the diagonal matrix D.
LLD — DOUBLE PRECISION array, dimension (N-1). The (N-1) elements L(i)*L(i)*D(i).
IFIRST — INTEGER. The index of the first eigenvalue to be computed.
ILAST — INTEGER. The index of the last eigenvalue to be computed.
RTOL1 — DOUBLE PRECISION.
RTOL2 — DOUBLE PRECISION. Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) where GAP is the (estimated) distance to the nearest eigenvalue. that's the convergence test, memorize it or perish. 💅
OFFSET — INTEGER. Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are to be used.
W
W — DOUBLE PRECISION array, dimension (N). On the way in, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are your best guesses at the eigenvalues of L D L^T, indexed IFIRST through ILAST. On the way out? bestie those estimates are ✨refined✨.
WGAP — DOUBLE PRECISION array, dimension (N-1). On input, the (estimated) gaps between consecutive eigenvalues of L D L^T, like WGAP(I-OFFSET) is the gap between eigenvalues I and I+1. Heads up tho: if IFIRST = ILAST then WGAP(IFIRST-OFFSET) must be set to ZERO. On output, these gaps get refined too. 💅
WERR — DOUBLE PRECISION array, dimension (N). On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined. slay.
WORK — DOUBLE PRECISION array, dimension (2*N). Workspace, aka the lil scratch pad 🎀.
IWORK — INTEGER array, dimension (2*N). Workspace again, integer flavored.
PIVMIN — DOUBLE PRECISION. The minimum pivot in the Sturm sequence. iconic little number.
SPDIAM — DOUBLE PRECISION. The spectral diameter of the matrix. 🌈
TWIST — INTEGER. The twist index for the twisted factorization that's used for the negcount:
- TWIST = N: compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
- TWIST = 1: compute negcount from L D L^T - LAMBDA I = U- D- U-^T
- TWIST = R: compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
INFO — INTEGER. Error flag. that's it that's the whole job 🫶

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors:

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

🦋 subroutine dlarrc

subroutine dlarrc (character JOBT, integer N, double precision VL, double
    precision VU, double precision, dimension( * ) D, double precision,
    dimension( * ) E, double precision PIVMIN, integer EIGCNT, integer LCNT,
    integer RCNT, integer INFO)

DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix. no cap, that's the whole vibe check 📊

Purpose:

Find the number of eigenvalues of the symmetric tridiagonal matrix T that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T if JOBT = 'L'.

Parameters

JOBT — CHARACTER*1
- = 'T': Compute Sturm count for matrix T.
- = 'L': Compute Sturm count for matrix L D L^T.
N — INTEGER. The order of the matrix. N > 0.
VL — DOUBLE PRECISION. The lower bound for the eigenvalues.
VU — DOUBLE PRECISION. The upper bound for the eigenvalues.
D — DOUBLE PRECISION array, dimension (N).
- JOBT = 'T': The N diagonal elements of the tridiagonal matrix T.
- JOBT = 'L': The N diagonal elements of the diagonal matrix D.
E — DOUBLE PRECISION array, dimension (N).
- JOBT = 'T': The N-1 offdiagonal elements of the matrix T.
- JOBT = 'L': The N-1 offdiagonal elements of the matrix L.
PIVMIN — DOUBLE PRECISION. The minimum pivot in the Sturm sequence for T.
EIGCNT — INTEGER. The number of eigenvalues of the symmetric tridiagonal matrix T that are in the interval (VL,VU].
LCNT — INTEGER.
RCNT — INTEGER. The left and right negcounts of the interval.
INFO — INTEGER.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

💜 subroutine dlarrd

subroutine dlarrd (character RANGE, character ORDER, integer N, double
    precision VL, double precision VU, integer IL, integer IU, double
    precision, dimension( * ) GERS, double precision RELTOL, double
    precision, dimension( * ) D, double precision, dimension( * ) E, double
    precision, dimension( * ) E2, double precision PIVMIN, integer NSPLIT,
    integer, dimension( * ) ISPLIT, integer M, double precision, dimension(
    * ) W, double precision, dimension( * ) WERR, double precision WL,
    double precision WU, integer, dimension( * ) IBLOCK, integer, dimension(
    * ) INDEXW, double precision, dimension( * ) WORK, integer, dimension( *
    ) IWORK, integer INFO)

DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. 💖

Purpose:

DLARRD computes the eigenvalues of a symmetric tridiagonal matrix T to suitable accuracy. This is an auxiliary code to be called from DSTEMR.

The user may ask for all eigenvalues, all eigenvalues in the half-open interval (VL, VU], or the IL-th through IU-th eigenvalues.

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. ⚠️

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966.

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 — DOUBLE PRECISION. 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 — DOUBLE PRECISION. 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'.
GERS — DOUBLE PRECISION array, dimension (2*N). The N Gerschgorin intervals (the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i))).
RELTOL — DOUBLE PRECISION. The minimum relative width of an interval. When an interval is narrower than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged. Note: this should always be at least radix*machine epsilon.
D — DOUBLE PRECISION array, dimension (N). The n diagonal elements of the tridiagonal matrix T.
E — DOUBLE PRECISION array, dimension (N-1). The (n-1) off-diagonal elements of the tridiagonal matrix T.
E2 — DOUBLE PRECISION array, dimension (N-1). The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
PIVMIN — DOUBLE PRECISION. The minimum pivot allowed in the Sturm sequence for T.
NSPLIT — INTEGER. The number of diagonal blocks in the matrix T. 1 <= NSPLIT <= N.
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.) 🧚
M — INTEGER. The actual number of eigenvalues found. 0 <= M <= N. (See also the description of INFO=2,3.)
W

`W` 💖

W is DOUBLE PRECISION array, dimension (N).

On exit, the first M elements of W will contain the eigenvalue approximations. DLARRD computes an interval I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue approximation is given as the interval midpoint:

W(j)= ( a_j + b_j)/2

The corresponding error is bounded by:

WERR(j) = abs( a_j - b_j)/2

`WERR` 🌸

WERR is DOUBLE PRECISION array, dimension (N).

The error bound on the corresponding eigenvalue approximation in W. bestie this is literally the "how sure are we" number, no cap.

`WL` 🦋

WL is DOUBLE PRECISION.

`WU` ✨

WU is DOUBLE PRECISION.

The interval (WL, WU] contains all the wanted eigenvalues.

If RANGE='V', then WL=VL and WU=VU.
If RANGE='A', then WL and WU are the global Gerschgorin bounds on the spectrum.
If RANGE='I', then WL and WU are computed by DLAEBZ from the index range specified.

`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. (DLARRD may use the remaining N-M elements as workspace.)

`INDEXW` 🎀

INDEXW is INTEGER array, dimension (N).

The indices of the eigenvalues within each block (submatrix); for example, INDEXW(i)= j and IBLOCK(i)=k imply that the i-th eigenvalue W(i) is the j-th eigenvalue in block k. it's giving "which one in which block" energy, periodt.

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

not the sloppy floating-point arithmetic jump scare 😭 if you see INFO = 4 your machine is lowkey embarrassing you fr fr.

✨ Internal Parameters ✨

`FUDGE` 🍓

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 — it's just here to keep the vibes (and the intervals) wide enough.

🌟 Contributors 🌟

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

👑 Author 👑

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

🦋 `subroutine dlarre` 🦋

subroutine dlarre (character RANGE, integer N, double precision VL, double
precision VU, integer IL, integer IU, double precision, dimension( * )
D, double precision, dimension( * ) E, double precision, dimension( * )
E2, double precision RTOL1, double precision RTOL2, double precision
SPLTOL, integer NSPLIT, integer, dimension( * ) ISPLIT, integer M,
double precision, dimension( * ) W, double precision, dimension( * )
WERR, double precision, dimension( * ) WGAP, integer, dimension( * )
IBLOCK, integer, dimension( * ) INDEXW, double precision, dimension( * )
GERS, double precision PIVMIN, double precision, dimension( * ) WORK,
integer, dimension( * ) IWORK, integer INFO)

DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.

🎯 Purpose 🎯

To find the desired eigenvalues of a given real symmetric tridiagonal matrix T, DLARRE sets any "small" off-diagonal elements to zero, and for each unreduced block T_i, it finds:

a suitable shift at one end of the block's spectrum,
the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
eigenvalues of each L_i D_i L_i^T.

The representations and eigenvalues found are then used by DSTEMR to compute the eigenvectors of T.

The accuracy varies depending on whether bisection is used to find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to conpute all and then discard any unwanted one.

As an added benefit, DLARRE also outputs the n Gerschgorin intervals for the matrices L_i D_i L_i^T.

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

`N`

N is INTEGER

The order of the matrix. N > 0.

`VL`

VL is DOUBLE PRECISION

If RANGE='V', the lower bound for the eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU.

If RANGE='I' or ='A', DLARRE computes bounds on the desired part of the spectrum.

`VU`

VU is DOUBLE PRECISION

If RANGE='V', the upper bound for the eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU.

If RANGE='I' or ='A', DLARRE computes bounds on the desired part of the spectrum.

`IL`

IL is INTEGER

If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N.

`IU`

IU is INTEGER

If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N.

`D`

D is DOUBLE PRECISION array, dimension (N)

On entry, the N diagonal elements of the tridiagonal matrix T.

On exit, the N diagonal elements of the diagonal matrices D_i.

`E`

E is DOUBLE PRECISION array, dimension (N)

On entry, the first (N-1) entries contain the subdiagonal elements of the tridiagonal matrix T; E(N) need not be set.

On exit, E contains the subdiagonal elements of the unit bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, contain the base points sigma_i on output.

`E2`

E2 is DOUBLE PRECISION array, dimension (N)

On entry, the first (N-1) entries contain the SQUARES of the subdiagonal elements of the tridiagonal matrix T; E2(N) need not be set.

On exit, the entries E2( ISPLIT( I ) ), 1 <= I <= NSPLIT, have been set to zero.

`RTOL1`

RTOL1 is DOUBLE PRECISION

`RTOL2`

RTOL2 is DOUBLE PRECISION

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

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

`SPLTOL`

SPLTOL is DOUBLE PRECISION

The threshold for splitting. if it's under this, it's getting split, sorry bestie 🫶

`NSPLIT`

NSPLIT is INTEGER

The number of blocks T splits into. 1 <= NSPLIT <= N.

`ISPLIT`

ISPLIT is INTEGER array, dimension (N)

The splitting points, at which T breaks up into blocks. The first block consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.

`M`

M is INTEGER

The total number of eigenvalues (of all L_i D_i L_i^T) found.

`W`

W is DOUBLE PRECISION array, dimension (N)

The first M elements contain the eigenvalues. The eigenvalues of each of the blocks, L_i D_i L_i^T, are sorted in ascending order (DLARRE may use the remaining N-M elements as workspace).

`WERR`

WERR is DOUBLE PRECISION array, dimension (N)

The error bound on the corresponding eigenvalue in W.

`WGAP`

WGAP — DOUBLE PRECISION array, dimension (N)

This one is the separation from the right neighbor eigenvalue in W. 💫

Real talk tho, the gap is only measured against eigenvalues of the same block, because every single block has its own representation tree. Exception bestie: at the right end of a block, we store the left gap instead. Asymmetric but that's the vibe.

IBLOCK — INTEGER array, dimension (N)

The indices of the blocks (submatrices) associated with the corresponding eigenvalues in W. So IBLOCK(i)=1 if eigenvalue W(i) belongs to the first block from the top, =2 if W(i) belongs to the second block, etc. Counting from the top, always. 👑

INDEXW — INTEGER array, dimension (N)

The indices of the eigenvalues within each block (submatrix). For example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the i-th eigenvalue W(i) is the 10-th eigenvalue in block 2. Two coordinates, one eigenvalue, it's giving address label energy. 🎯

GERS — DOUBLE PRECISION array, dimension (2*N)

The N Gerschgorin intervals, and the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i)). 📏

PIVMIN — DOUBLE PRECISION

The minimum pivot in the Sturm sequence for T.

WORK — DOUBLE PRECISION array, dimension (6*N)

Workspace. That's it, that's the whole job description.

IWORK — INTEGER array, dimension (5*N)

Workspace. Same energy, integer edition.

INFO — INTEGER
- = 0: successful exit 🎉
- > 0: a problem occurred in DLARRE.
- < 0: one of the called subroutines signaled an internal problem. You're gonna need to inspect the corresponding parameter IINFO for further information.
- =-1: problem in DLARRD.
- = 2: no base representation could be found in MAXTRY iterations. Increasing MAXTRY and recompilation might be a remedy.
- =-3: problem in DLARRB when computing the refined root representation for DLASQ2.
- =-4: problem in DLARRB when performing bisection on the desired part of the spectrum.
- =-5: problem in DLASQ2.
- =-6: problem in DLASQ2.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 Further Details 🌈

The base representations are required to suffer very little element growth, and consequently they define all their eigenvalues to high relative accuracy. That's the whole point, bestie — accuracy is the main character here.

🫶 Contributors 🫶

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

🧚 subroutine `dlarrf` 🧚

subroutine dlarrf (integer N, double precision, dimension( * ) D, double
    precision, dimension( * ) L, double precision, dimension( * ) LD,
    integer CLSTRT, integer CLEND, double precision, dimension( * ) W,
    double precision, dimension( * ) WGAP, double precision, dimension( * )
    WERR, double precision SPDIAM, double precision CLGAPL, double precision
    CLGAPR, double precision PIVMIN, double precision SIGMA, double
    precision, dimension( * ) DPLUS, double precision, dimension( * ) LPLUS,
    double precision, dimension( * ) WORK, integer INFO)

DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. She's a glow-up subroutine. ✨

🌸 Purpose 🌸

Given the initial representation L D L^T and its cluster of close eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... W( CLEND ), DLARRF finds a new relatively robust representation L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the eigenvalues of L(+) D(+) L(+)^T is relatively isolated. One is enough, and honestly? One is iconic.

💅 Parameters 💅

N — INTEGER. The order of the matrix (subblock, if the matrix split).
D — DOUBLE PRECISION array, dimension (N). The N diagonal elements of the diagonal matrix D.
L — DOUBLE PRECISION array, dimension (N-1). The (N-1) subdiagonal elements of the unit bidiagonal matrix L.
LD — DOUBLE PRECISION array, dimension (N-1). The (N-1) elements L(i)*D(i).
CLSTRT — INTEGER. The index of the first eigenvalue in the cluster.
CLEND — INTEGER. The index of the last eigenvalue in the cluster.
W — DOUBLE PRECISION array, dimension and the dimension is >= (CLEND-CLSTRT+1). The eigenvalue APPROXIMATIONS of L D L^T in ascending order. W( CLSTRT ) through W( CLEND ) form the cluster of relatively close eigenvalues.
WGAP — DOUBLE PRECISION array, dimension and the dimension is >= (CLEND-CLSTRT+1). The separation from the right neighbor eigenvalue in W.
WERR — DOUBLE PRECISION array, dimension and the dimension is >= (CLEND-CLSTRT+1). WERR contains the semiwidth of the uncertainty interval of the corresponding eigenvalue APPROXIMATION in W. It's the little bubble of "we're pretty sure it's in here somewhere". 🫧
SPDIAM — DOUBLE PRECISION. Estimate of the spectral diameter obtained from the Gerschgorin intervals.
CLGAPL — DOUBLE PRECISION. (That's the whole spec, no notes. She's mysterious. 💀)
CLGAPR — DOUBLE PRECISION. Absolute gap on each end of the cluster. Set by the calling routine to protect against shifts too close to eigenvalues outside the cluster. Protective queen behavior.
PIVMIN — DOUBLE PRECISION. The minimum pivot allowed in the Sturm sequence.
SIGMA — DOUBLE PRECISION. The shift used to form L(+) D(+) L(+)^T.
DPLUS — DOUBLE PRECISION array, dimension (N). The N diagonal elements of the diagonal matrix D(+).
LPLUS — DOUBLE PRECISION array, dimension (N-1). The first (N-1) elements of LPLUS contain the subdiagonal elements of the unit bidiagonal matrix L(+).
WORK — DOUBLE PRECISION array, dimension (2*N). Workspace.
INFO — INTEGER. Signals processing OK (=0) or failure (=1).

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors 🫶

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

🦋 subroutine `dlarrj` 🦋

subroutine dlarrj (integer N, double precision, dimension( * ) D, double
    precision, dimension( * ) E2, integer IFIRST, integer ILAST, double
    precision RTOL, integer OFFSET, double precision, dimension( * ) W,
    double precision, dimension( * ) WERR, double precision, dimension( * )
    WORK, integer, dimension( * ) IWORK, double precision PIVMIN, double
    precision SPDIAM, integer INFO)

DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T. Glow-up for numbers. 🔥

🌸 Purpose 🌸

Given the initial eigenvalue approximations of T, DLARRJ does bisection to refine the eigenvalues of T, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial guesses for these eigenvalues are input in W, and the corresponding estimate of the error in these guesses goes in WERR. During bisection, intervals [left, right] are maintained by storing their mid-points and semi-widths in the arrays W and WERR respectively. Squeeze, refine, repeat. 💅

💅 Parameters 💅

N — INTEGER. The order of the matrix.
D — DOUBLE PRECISION array, dimension (N). The N diagonal elements of T.
E2 — DOUBLE PRECISION array, dimension (N-1). The squares of the (N-1) subdiagonal elements of T.
IFIRST — INTEGER. The index of the first eigenvalue to be computed.
ILAST — INTEGER. The index of the last eigenvalue to be computed.
RTOL — DOUBLE PRECISION. Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|). That's the finish line, no cap.
OFFSET — INTEGER. Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are to be used.
W — DOUBLE PRECISION array, dimension (N). On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. On output, these estimates are refined. ✨
WERR — DOUBLE PRECISION array, dimension (N). On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined.
WORK — DOUBLE PRECISION array, dimension (2*N). Workspace.
IWORK — INTEGER array, dimension (2*N). Workspace.
PIVMIN — DOUBLE PRECISION. The minimum pivot in the Sturm sequence for T.
SPDIAM — DOUBLE PRECISION. The spectral diameter of T.
INFO — INTEGER. Error flag.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🫶 Contributors 🫶

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

💫 `dlarrk`

subroutine dlarrk (integer N, integer IW, double precision GL, double precision GU, double precision, dimension( * ) D, double precision, dimension( * ) E2, double precision PIVMIN, double precision RELTOL, double precision W, double precision WERR, integer INFO)

DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. ✨

Purpose 💖

okay so bestie, here's the tea: DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. It's an auxiliary code, so it's not out here living its main character life on its own — it gets called from DSTEMR. 🫶

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966. (an actual icon from 1966, we stan)

Parameters 🎀

N — N is INTEGER. The order of the tridiagonal matrix T. N >= 0.
IW — IW is INTEGER. The index of the eigenvalues to be returned.
GL — GL is DOUBLE PRECISION
GU — GU is DOUBLE PRECISION. An upper and a lower bound on the eigenvalue.
D — D is DOUBLE PRECISION array, dimension (N). The n diagonal elements of the tridiagonal matrix T.
E2 — E2 is DOUBLE PRECISION array, dimension (N-1). The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
PIVMIN — PIVMIN is DOUBLE PRECISION. The minimum pivot allowed in the Sturm sequence for T.
RELTOL — RELTOL is DOUBLE PRECISION. The minimum relative width of an interval. When an interval is narrower than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged. Note: this should always be at least radix*machine epsilon.
W — W is DOUBLE PRECISION
WERR — WERR is DOUBLE PRECISION. The error bound on the corresponding eigenvalue approximation in W.
INFO — INFO is INTEGER
- = 0: Eigenvalue converged
- = -1: Eigenvalue did NOT converge (not me getting the sad ending 😭)

Internal Parameters 🌈

FUDGE — DOUBLE PRECISION, default = 2. A "fudge factor" to widen the Gershgorin intervals. (the fudge factor!! the moment!!)

Author 👑

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

🦋 `dlarrr`

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

DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues.

Purpose 💅

Perform tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues. (basically it's the vibe check before we spend the bag 💸)

Parameters ✨

N — N is INTEGER. The order of the matrix. N > 0.
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). On entry, the first (N-1) entries contain the subdiagonal elements of the tridiagonal matrix T; E(N) is set to ZERO.
INFO — INFO is INTEGER
- INFO = 0(default) : the matrix warrants computations preserving relative accuracy.
- INFO = 1 : the matrix warrants computations guaranteeing only absolute accuracy.

Author 👑

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

Contributors 🫶

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

🌸 `dlartgp`

subroutine dlartgp (double precision F, double precision G, double precision CS, double precision SN, double precision R)

DLARTGP generates a plane rotation so that the diagonal is nonnegative.

Purpose 💖

DLARTGP generates a plane rotation so that

   [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1.
   [ -SN  CS  ]     [ G ]     [ 0 ]

This is a slower, more accurate version of the Level 1 BLAS routine DROTG, with the following other differences:

F and G are unchanged on return. (they're just standing there, unbothered, moisturized ✨)
If G=0, then CS=(+/-)1 and SN=0.
If F=0 and (G .ne. 0), then CS=0 and SN=(+/-)1.

The sign is chosen so that R >= 0.

Parameters 🎀

F — F is DOUBLE PRECISION. The first component of vector to be rotated.
G — G is DOUBLE PRECISION. The second component of vector to be rotated.
CS — CS is DOUBLE PRECISION. The cosine of the rotation.
SN — SN is DOUBLE PRECISION. The sine of the rotation.
R — R is DOUBLE PRECISION. The nonzero component of the rotated vector.

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

Author 👑

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

💜 `dlaruv`

subroutine dlaruv (integer, dimension( 4 ) ISEED, integer N, double precision, dimension( n ) X)

DLARUV returns a vector of n random real numbers from a uniform distribution.

Purpose ⭐

DLARUV returns a vector of n random real numbers from a uniform (0,1) distribution (n <= 128).

This is an auxiliary routine called by DLARNV and ZLARNV.

Parameters 🍓

ISEED — ISEED is INTEGER array, dimension (4). On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. On exit, the seed is updated.
N — N is INTEGER. The number of random numbers to be generated. N <= 128.
X — X is DOUBLE PRECISION array, dimension (N). The generated random numbers.

Author 👑

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

Further Details 🔥

This routine uses a multiplicative congruential method with modulus 2**48 and multiplier 33952834046453 (see G.S.Fishman, 'Multiplicative congruential random number generators with modulus 2**b: an exhaustive analysis for b = 32 and a partial analysis for b = 48', Math. Comp. 189, pp 331-344, 1990).

48-bit integers are stored in 4 integer array elements with 12 bits per element. Hence the routine is portable across machines with integers of 32 bits or more. (portability ate and left no crumbs 💅)

🧚 `dlas2`

subroutine dlas2 (double precision F, double precision G, double precision H, double precision SSMIN, double precision SSMAX)

DLAS2 computes singular values of a 2-by-2 triangular matrix.

Purpose 🌈

DLAS2 computes the singular values of the 2-by-2 matrix

   [  F   G  ]
   [  0   H  ].

On return, SSMIN is the smaller singular value and SSMAX is the larger singular value.

Parameters 💫

F — F is DOUBLE PRECISION. The (1,1) element of the 2-by-2 matrix.
G — G is DOUBLE PRECISION. The (1,2) element of the 2-by-2 matrix.
H — H is DOUBLE PRECISION. The (2,2) element of the 2-by-2 matrix.
SSMIN — SSMIN is DOUBLE PRECISION. The smaller singular value.
SSMAX — SSMAX is DOUBLE PRECISION. The larger singular value.

Author 👑

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

Further Details 🔍

Barring over/underflow, all output quantities are correct to within a few units in the last place (ulps), even in the absence of a guard digit in addition/subtraction.

In IEEE arithmetic, the code works correctly if one matrix element is infinite. (infinite? bestie the code said "I can handle that" 😌)

Overflow will not occur unless the largest singular value itself overflows, or is within a few ulps of overflow. (On machines with partial overflow, like the Cray, overflow may occur if the largest singular value is within a factor of 2 of overflow.)

Underflow is harmless if underflow is gradual. Otherwise, results may correspond to a matrix modified by perturbations of size near the underflow threshold.

👑 `dlascl`

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

DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.

Purpose 💖

DLASCL multiplies the M by N real matrix A by the real scalar CTO/CFROM. This is done without over/underflow as long as the final result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that A may be full, upper triangular, lower triangular, upper Hessenberg, or banded. (it's giving "one scalar, so many storage formats" 🎀)

Parameters ✨

TYPE — TYPE is CHARACTER*1. TYPE indices the storage type of the input matrix.
- = 'G': A is a full matrix.
- = 'L': A is a lower triangular matrix.
- = 'U': A is an upper triangular matrix.
- = 'H': A is an upper Hessenberg matrix.
- = 'B': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the lower half stored.
- = 'Q': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the upper half stored.
- = 'Z': A is a band matrix with lower bandwidth KL and upper bandwidth KU. See DGBTRF for storage details.
KL
KL is INTEGER — the lower bandwidth of A. Only referenced if TYPE = 'B', 'Q' or 'Z'. 💅
KU is INTEGER — the upper bandwidth of A. Only referenced if TYPE = 'B', 'Q' or 'Z'. same energy bestie 🫶
CFROM is DOUBLE PRECISION
CTO is DOUBLE PRECISION

The matrix A is multiplied by CTO/CFROM. A(I,J) is computed without over/underflow if the final result CTO*A(I,J)/CFROM can be represented without over/underflow. CFROM must be nonzero. like fr fr do not pass zero, she will simply not survive 😭

M is INTEGER — the number of rows of the matrix A. M >= 0.
N is INTEGER — the number of columns of the matrix A. N >= 0.
A is DOUBLE PRECISION array, dimension (LDA,N) — the matrix to be multiplied by CTO/CFROM. See TYPE for the storage type. ✨
LDA is INTEGER — the leading dimension of the array A.
- If TYPE = 'G', 'L', 'U', 'H', LDA >= max(1,M);
- TYPE = 'B', LDA >= KL+1;
- TYPE = 'Q', LDA >= KU+1;
- TYPE = 'Z', LDA >= 2*KL+KU+1.

pick your lane and meet the minimum, no cutting corners bestie 💅

INFO is INTEGER
- 0 — successful exit, we ate 💖
- <0 — if INFO = -i, the i-th argument had an illegal value. not iconic, fix ur input 😭

✨ Author ✨

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🦋 `subroutine dlasd0`

subroutine dlasd0 (integer N, integer SQRE, double precision, dimension( * )
    D, double precision, dimension( * ) E, double precision, dimension( ldu,
    * ) U, integer LDU, double precision, dimension( ldvt, * ) VT, integer
    LDVT, integer SMLSIZ, integer, dimension( * ) IWORK, double precision,
    dimension( * ) WORK, integer INFO)

DLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc. 💫

🎀 Purpose:

Using a divide and conquer approach, DLASD0 computes the singular value decomposition (SVD) of a real upper bidiagonal N-by-M matrix B with diagonal D and offdiagonal E, where M = N + SQRE. The algorithm computes orthogonal matrices U and VT such that B = U * S * VT. The singular values S are overwritten on D. no cap, D gets recycled into S, that's the glow-up ✨

A related subroutine, DLASDA, computes only the singular values, and optionally, the singular vectors in compact form. 👑

🌸 Parameters

N is INTEGER

On entry, the row dimension of the upper bidiagonal matrix. This is also the dimension of the main diagonal array D.

SQRE is INTEGER

Specifies the column dimension of the bidiagonal matrix.

= 0: The bidiagonal matrix has column dimension M = N;
= 1: The bidiagonal matrix has column dimension M = N+1;

D is DOUBLE PRECISION array, dimension (N)

On entry D contains the main diagonal of the bidiagonal matrix. On exit D, if INFO = 0, contains its singular values. 💅

E is DOUBLE PRECISION array, dimension (M-1)

Contains the subdiagonal entries of the bidiagonal matrix. On exit, E has been destroyed. rip E, gone but not forgotten 💀

U is DOUBLE PRECISION array, dimension (LDU, N)

On exit, U contains the left singular vectors.

LDU is INTEGER

On entry, leading dimension of U.

VT is DOUBLE PRECISION array, dimension (LDVT, M)

On exit, VT**T contains the right singular vectors.

LDVT is INTEGER

On entry, leading dimension of VT.

SMLSIZ is INTEGER

On entry, maximum size of the subproblems at the bottom of the computation tree. 🌳

IWORK is INTEGER array, dimension (8*N)
WORK is DOUBLE PRECISION array, dimension (3*M**2+2*M)
INFO is INTEGER
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, a singular value did not converge. she's struggling, not the vibe 😭

✨ 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 dlasd1`

subroutine dlasd1 (integer NL, integer NR, integer SQRE, double precision,
    dimension( * ) D, double precision ALPHA, double precision BETA, double
    precision, dimension( ldu, * ) U, integer LDU, double precision,
    dimension( ldvt, * ) VT, integer LDVT, integer, dimension( * ) IDXQ,
    integer, dimension( * ) IWORK, double precision, dimension( * ) WORK,
    integer INFO)

DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc. 💫

🎀 Purpose:

DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.

A related subroutine DLASD7 handles the case in which the singular values (and the singular vectors in factored form) are desired. 👑

DLASD1 computes the SVD as follows:

              ( D1(in)    0    0       0 )
  B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
              (   0       0   D2(in)   0 )

    = U(out) * ( D(out) 0) * VT(out)

where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros elsewhere; and the entry b is empty if SQRE = 0. 🧚

The left singular vectors of the original matrix are stored in U, and the transpose of the right singular vectors are stored in VT, and the singular values are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem when there are multiple singular values or when there are zeros 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 DLASD2.
The second stage consists of calculating the updated singular values. This is done by finding the square roots of the roots of the secular equation via the routine DLASD4 (as called by DLASD3). This routine also calculates the singular vectors of the current problem.
The final stage consists of computing the updated singular vectors directly using the updated singular values. The singular vectors for the current problem are multiplied with the singular vectors from the overall problem.

three stages, three vibes, all iconic 💅✨

🌸 Parameters

NL is INTEGER

The row dimension of the upper block. NL >= 1.

NR is INTEGER

The row dimension of the lower block. NR >= 1.

SQRE is INTEGER
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has row dimension N = NL + NR + 1, and column dimension M = N + SQRE.

D is DOUBLE PRECISION array, dimension (N = NL+NR+1).

On entry D(1:NL,1:NL) contains the singular values of the upper block; and D(NL+2:N) contains the singular values of the lower block. On exit D(1:N) contains the singular values of the modified matrix. glow-up complete ✨

ALPHA is DOUBLE PRECISION

Contains the diagonal element associated with the added row.

BETA is DOUBLE PRECISION

Contains the off-diagonal element associated with the added row.

U is DOUBLE PRECISION array, dimension (LDU,N)

On entry U(1:NL, 1:NL) contains the left singular vectors of the upper block; U(NL+2:N, NL+2:N) contains the left singular vectors of the lower block. On exit U contains the left singular vectors of the bidiagonal matrix.

LDU is INTEGER

The leading dimension of the array U. LDU >= max( 1, N ).

VT is DOUBLE PRECISION array, dimension (LDVT,M) where M = N + SQRE.

On entry VT(1:NL+1, 1:NL+1)**T contains the right singular vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains the right singular vectors of the lower block. On exit VT**T contains the right singular vectors of the bidiagonal matrix.

LDVT is INTEGER

The leading dimension of the array VT. LDVT >= max( 1, M ).

IDXQ is INTEGER array, dimension (N)

This contains the permutation which will reintegrate the subproblem just solved back into sorted order, i.e. D( IDXQ( I = 1, N ) ) will be in ascending order. sorting era 🔥

IWORK is INTEGER array, dimension ( 4 * N )
WORK is DOUBLE PRECISION array, dimension ( 3*M**2 + 2*M )
INFO is INTEGER
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, a singular value did not converge. she's struggling, not the vibe 😭

✨ 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

🦋 `dlasd2` — the merge-and-deflate bestie 🦋

   subroutine dlasd2 (integer NL, integer NR, integer SQRE, integer K, double
       precision, dimension( * ) D, double precision, dimension( * ) Z, double
       precision ALPHA, double precision BETA, double precision, dimension(
       ldu, * ) U, integer LDU, double precision, dimension( ldvt, * ) VT,
       integer LDVT, double precision, dimension( * ) DSIGMA, double precision,
       dimension( ldu2, * ) U2, integer LDU2, double precision, dimension(
       ldvt2, * ) VT2, integer LDVT2, integer, dimension( * ) IDXP, integer,
       dimension( * ) IDX, integer, dimension( * ) IDXC, integer, dimension( *
       ) IDXQ, integer, dimension( * ) COLTYP, integer INFO)

DLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc.

✨ Purpose ✨

So like, DLASD2 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem — yes bestie, deflate, it's giving "make it smaller, we don't need all that". 💅

There are two ways deflation can happen:

when two or more singular values 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. Periodt.

DLASD2 is called from DLASD1.

🌸 Parameters 🌸

NL — INTEGER. The row dimension of the upper block. NL >= 1.
NR — INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — INTEGER. Okay so:
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns. Slay.

K — INTEGER. Contains the dimension of the non-deflated matrix. This is the order of the related secular equation. 1 <= K <=N.
D — DOUBLE PRECISION array, dimension(N). On entry D contains the singular values of the two submatrices to be combined. On exit D contains the trailing (N-K) updated singular values (those which were deflated) sorted into increasing order. Iconic behavior honestly.
Z — DOUBLE PRECISION array, dimension(N). On exit Z contains the updating row vector in the secular equation.
ALPHA — DOUBLE PRECISION. Contains the diagonal element associated with the added row.
BETA — DOUBLE PRECISION. Contains the off-diagonal element associated with the added row.
U — DOUBLE PRECISION array, dimension(LDU,N). On entry U contains the left singular vectors of two submatrices in the two square blocks with corners at (1,1), (NL, NL), and (NL+2, NL+2), (N,N). On exit U contains the trailing (N-K) updated left singular vectors (those which were deflated) in its last N-K columns. 🫶
LDU — INTEGER. The leading dimension of the array U. LDU >= N.
VT — DOUBLE PRECISION array, dimension(LDVT,M). On entry VT**T contains the right singular vectors of two submatrices in the two square blocks with corners at (1,1), (NL+1, NL+1), and (NL+2, NL+2), (M,M). On exit VT**T contains the trailing (N-K) updated right singular vectors (those which were deflated) in its last N-K columns. In case SQRE =1, the last row of VT spans the right null space. Not me reading VT**T and forgetting it's the transpose again 😭
LDVT — INTEGER. The leading dimension of the array VT. LDVT >= M.
DSIGMA — DOUBLE PRECISION array, dimension (N). Contains a copy of the diagonal elements (K-1 singular values and one zero) in the secular equation.
U2 — DOUBLE PRECISION array, dimension(LDU2,N). Contains a copy of the first K-1 left singular vectors which will be used by DLASD3 in a matrix multiply (DGEMM) to solve for the new left singular vectors. U2 is arranged into four blocks:
- The first block contains a column with 1 at NL+1 and zero everywhere else;
- the second block contains non-zero entries only at and above NL;
- the third contains non-zero entries only below NL+1;
- and the fourth is dense.

Four blocks. Very structured. We love a girl with a plan. ✨

LDU2 — INTEGER. The leading dimension of the array U2. LDU2 >= N.
VT2 — DOUBLE PRECISION array, dimension(LDVT2,N). VT2**T contains a copy of the first K right singular vectors which will be used by DLASD3 in a matrix multiply (DGEMM) to solve for the new right singular vectors. VT2 is arranged into three blocks:
- The first block contains a row that corresponds to the special 0 diagonal element in SIGMA;
- the second block contains non-zeros only at and before NL +1;
- the third block contains non-zeros only at and after NL +2.
LDVT2 — INTEGER. The leading dimension of the array VT2. LDVT2 >= M.
IDXP — INTEGER array, dimension(N). This will contain the permutation used to place deflated values of D at the end of the array. On output IDXP(2:K) points to the nondeflated D-values and IDXP(K+1:N) points to the deflated singular values.
IDX — INTEGER array, dimension(N). This will contain the permutation used to sort the contents of D into ascending order.
IDXC — INTEGER array, dimension(N). This will contain the permutation used to arrange the columns of the deflated U matrix into three groups:
- the first group contains non-zero entries only at and above NL,
- the second contains non-zero entries only below NL+2,
- and the third is dense.
IDXQ — INTEGER array, dimension(N). This contains the permutation which separately sorts the two sub-problems in D into ascending order. Note that entries in the first half of this permutation must first be moved one position backward; and entries in the second half must first have NL+1 added to their values. Two-step prep, iykyk. 💅
COLTYP — INTEGER array, dimension(N). As workspace, this will contain a label which will indicate which of the following types a column in the U2 matrix or a row in the VT2 matrix is:
1. non-zero in the upper half only
2. non-zero in the lower half only
3. dense
4. deflated

On exit, it is an array of dimension 4, with COLTYP(I) being the dimension of the I-th type columns.

INFO — INTEGER. The receipts:
- = 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 🎀

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

🌈 `dlasd3` — square roots of the secular equation, but make it fashion 🌈

   subroutine dlasd3 (integer NL, integer NR, integer SQRE, integer K, double
       precision, dimension( * ) D, double precision, dimension( ldq, * ) Q,
       integer LDQ, double precision, dimension( * ) DSIGMA, double precision,
       dimension( ldu, * ) U, integer LDU, double precision, dimension( ldu2, *
       ) U2, integer LDU2, double precision, dimension( ldvt, * ) VT, integer
       LDVT, double precision, dimension( ldvt2, * ) VT2, integer LDVT2,
       integer, dimension( * ) IDXC, integer, dimension( * ) CTOT, double
       precision, dimension( * ) Z, integer INFO)

DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc.

✨ Purpose ✨

DLASD3 finds all the square roots of the roots of the secular equation, as defined by the values in D and Z. It makes the appropriate calls to DLASD4 and then updates the singular vectors by matrix multiplication. She's a caller AND a doer, main character energy fr fr. 👑

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 XMP, Cray YMP, Cray C 90, or Cray 2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. So like… probably fine bestie, no cap. 🙈

DLASD3 is called from DLASD1.

🌸 Parameters 🌸

NL — INTEGER. The row dimension of the upper block. NL >= 1.
NR — INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — INTEGER. Here we go again:
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns.

K — INTEGER. The size of the secular equation, 1 =< K = < N. (Yes, that's how it's written in the source. We keep it. 🫶)
D — DOUBLE PRECISION array, dimension(K). On exit the square roots of the roots of the secular equation, in ascending order. Sorted, tidy, organized — the dream. ✨
Q — DOUBLE PRECISION array, dimension (LDQ,K). That's it, that's the whole note. 💅
LDQ — INTEGER. The leading dimension of the array Q. LDQ >= K.
DSIGMA — 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.
U — DOUBLE PRECISION array, dimension (LDU, N). The last N - K columns of this matrix contain the deflated left singular vectors.
LDU — INTEGER. The leading dimension of the array U. LDU >= N.
U2 — DOUBLE PRECISION array, dimension (LDU2, N). The first K columns of this matrix contain the non-deflated left singular vectors for the split problem.
LDU2 —
LDU2 — an INTEGER, bestie 💅 and it's the leading dimension of the array U2. So like, LDU2 >= N. no cap.
VT — DOUBLE PRECISION array, dimension (LDVT, M). The last M - K columns of VT**T are holding the deflated right singular vectors, and honestly? they ate. ✨
LDVT — INTEGER. The leading dimension of the array VT. LDVT >= N, periodt.
VT2 — DOUBLE PRECISION array, dimension (LDVT2, N). The first K columns of VT2**T are carrying the non-deflated right singular vectors for the split problem 🎀
LDVT2 — INTEGER. The leading dimension of the array VT2. LDVT2 >= N.
IDXC — INTEGER array, dimension ( N ). It's the permutation used to arrange the columns of U (and rows of VT) into three groups, and here's the tea:
- the first group has non-zero entries only at and above (or before) NL +1
- the second has non-zero entries only at and below (or after) NL+2
- the third is dense, no notes 🌈

The first column of U and the row of VT are treated separately, however 🙈

Also, and this part matters: the rows of the singular vectors found by DLASD4 must be likewise permuted before the matrix multiplies can take place. iykyk.

CTOT — INTEGER array, dimension ( 4 ). It's a count of the total number of the various types of columns in U (or rows in VT), as described in IDXC. The fourth column type is any column which has been deflated. 💜
Z — DOUBLE PRECISION array, dimension (K). The first K elements of this array contain the components of the deflation-adjusted updating row vector. icon behaviour fr.
INFO — INTEGER, and here's the vibe check:
- = 0: successful exit 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value 💀
- > 0: if INFO = 1, a singular value did not converge 😭

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

subroutine dlasd4 (integer N, integer I, double precision, dimension( * ) D,
    double precision, dimension( * ) Z, double precision, dimension( * )
    DELTA, double precision RHO, double precision SIGMA, double precision,
    dimension( * ) WORK, integer INFO)

DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. It's used by dbdsdc 🦋

✨ Purpose ✨

This subroutine computes the square root of the I-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix whose entries are given as the squares of the corresponding entries in the array d, and that

       0 <= 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 ) * 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. cute, efficient, slay.

✨ 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 DOUBLE PRECISION array, dimension ( N ). The original eigenvalues. It is assumed that they are in order, 0 <= D(I) < D(J) for I < J.
Z — Z is DOUBLE PRECISION array, dimension ( N ). The components of the updating vector.
DELTA — DELTA is DOUBLE PRECISION array, dimension ( N ). If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th component. If N = 1, then DELTA(1) = 1. The vector DELTA contains the information necessary to construct the (singular) eigenvectors.
RHO — RHO is DOUBLE PRECISION. The scalar in the symmetric updating formula.
SIGMA — SIGMA is DOUBLE PRECISION. The computed sigma_I, the I-th updated eigenvalue.
WORK — WORK is DOUBLE PRECISION array, dimension ( N ). If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th component. If N = 1, then WORK( 1 ) = 1.
INFO — INFO is INTEGER, and the tea is:
- = 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 dlasd5 ✨

subroutine dlasd5 (integer I, double precision, dimension( 2 ) D, double
    precision, dimension( 2 ) Z, double precision, dimension( 2 ) DELTA,
    double precision RHO, double precision DSIGMA, double precision,
    dimension( 2 ) WORK)

DLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. It's used by sbdsdc 🧚

✨ Purpose ✨

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

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

The diagonal entries in the array D are assumed to satisfy

                   0 <= 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 DOUBLE PRECISION array, dimension ( 2 ). The original eigenvalues. We assume 0 <= D(1) < D(2).
Z — Z is DOUBLE PRECISION array, dimension ( 2 ). The components of the updating vector.
DELTA — DELTA is DOUBLE PRECISION array, dimension ( 2 ). Contains (D(j) - sigma_I) in its j-th component. The vector DELTA contains the information necessary to construct the eigenvectors.
RHO — RHO is DOUBLE PRECISION. The scalar in the symmetric updating formula.
DSIGMA — DSIGMA is DOUBLE PRECISION. The computed sigma_I, the I-th updated eigenvalue.
WORK — WORK is DOUBLE PRECISION array, dimension ( 2 ). WORK contains (D(j) + sigma_I) in its j-th component.

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

subroutine dlasd6 (integer ICOMPQ, integer NL, integer NR, integer SQRE,
    double precision, dimension( * ) D, double precision, dimension( * ) VF,
    double precision, dimension( * ) VL, double precision ALPHA, double
    precision BETA, integer, dimension( * ) IDXQ, integer, dimension( * )
    PERM, integer GIVPTR, integer, dimension( ldgcol, * ) GIVCOL, integer
    LDGCOL, double precision, dimension( ldgnum, * ) GIVNUM, integer LDGNUM,
    double precision, dimension( ldgnum, * ) POLES, double precision,
    dimension( * ) DIFL, double precision, dimension( * ) DIFR, double
    precision, dimension( * ) Z, integer K, double precision C, double
    precision S, double precision, dimension( * ) WORK, integer, dimension(
    * ) IWORK, integer INFO)

DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc ✨

✨ Purpose ✨

DLASD6 computes the SVD of an updated upper bidiagonal matrix B obtained by merging two smaller ones by appending a row. This routine is used only for the problem which requires all singular values and optionally singular vector matrices in factored form. B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. A related subroutine, DLASD1, handles the case in which all singular values and singular vectors of the bidiagonal matrix are desired.

DLASD6 computes the SVD as follows:

                      ( D1(in)    0    0       0 )
          B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
                      (   0       0   D2(in)   0 )

            = U(out) * ( D(out) 0) * VT(out)

where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros elsewhere; and the entry b is empty if SQRE = 0.

The singular values of B can be computed using D1, D2, the first components of all the right singular vectors of the lower block, and the last components of all the right singular vectors of the upper block. These components are stored and updated in VF and VL, respectively, in DLASD6. Hence U and VT are not explicitly referenced.

The singular values are stored in D. The algorithm consists of two stages:

The first stage consists of deflating the size of the problem when there are multiple singular values 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 DLASD7. 💅
The second stage consists of calculating the updated singular values. This is done by finding the roots of the secular equation via the routine DLASD4 (as called by DLASD8). This routine also updates VF and VL and computes the distances between the updated singular values and the old singular values. 🔥

DLASD6 is called from DLASDA.

✨ Parameters ✨

ICOMPQ
ICOMPQ — it's an INTEGER that decides whether singular vectors are computed in factored form, bestie. = 0: compute singular values only. = 1: compute singular vectors in factored form as well. Choose your fighter. ✨
NL — INTEGER, the row dimension of the upper block. NL >= 1. Respect the minimum. 💅
NR — INTEGER, the row dimension of the lower block. NR >= 1. Same energy. 🦋
SQRE — INTEGER. = 0: the lower block is an NR-by-NR square matrix. = 1: the lower block is an NR-by-(NR+1) rectangular matrix. The bidiagonal matrix has row dimension N = NL + NR + 1, and column dimension M = N + SQRE. Math is mathing. 📐
D — DOUBLE PRECISION array, dimension ( NL+NR+1 ). On entry, D(1:NL,1:NL) contains the singular values of the upper block, and D(NL+2:N) contains the singular values of the lower block. On exit, D(1:N) contains the singular values of the modified matrix. She's updated. 💖
VF — DOUBLE PRECISION array, dimension ( M ). On entry, VF(1:NL+1) contains the first components of all right singular vectors of the upper block; and VF(NL+2:M) contains the first components of all right singular vectors of the lower block. On exit, VF contains the first components of all right singular vectors of the bidiagonal matrix. Firsties! ✨
VL — DOUBLE PRECISION array, dimension ( M ). On entry, VL(1:NL+1) contains the last components of all right singular vectors of the upper block; and VL(NL+2:M) contains the last components of all right singular vectors of the lower block. On exit, VL contains the last components of all right singular vectors of the bidiagonal matrix. Last but not least. 🎀
ALPHA — DOUBLE PRECISION. Contains the diagonal element associated with the added row. Main character diagonal. 👑
BETA — DOUBLE PRECISION. Contains the off-diagonal element associated with the added row. The sidekick. 💜
IDXQ — INTEGER array, dimension ( N ). This contains the permutation which will reintegrate the subproblem just solved back into sorted order, i.e. D( IDXQ( I = 1, N ) ) will be in ascending order. Sorting? We love to see it. 🔥
PERM — INTEGER array, dimension ( N ). The permutations (from deflation and sorting) to be applied to each block. Not referenced if ICOMPQ = 0. So if you're not doing vectors, she's ghost. 👻
GIVPTR — INTEGER. The number of Givens rotations which took place in this subproblem. Not referenced if ICOMPQ = 0. Count them if you care. 💅
GIVCOL — INTEGER array, dimension ( LDGCOL, 2 ). Each pair of numbers indicates a pair of columns to take place in a Givens rotation. Not referenced if ICOMPQ = 0. Two columns, one rotation, iconic. 🌀
LDGCOL — INTEGER. Leading dimension of GIVCOL, must be at least N. Don't lowball it. 🚫
GIVNUM — DOUBLE PRECISION array, dimension ( LDGNUM, 2 ). Each number indicates the C or S value to be used in the corresponding Givens rotation. Not referenced if ICOMPQ = 0. Trig values said present. 📐
LDGNUM — INTEGER. The leading dimension of GIVNUM and POLES, must be at least N. Same rule, same energy. ✨
POLES — DOUBLE PRECISION array, dimension ( LDGNUM, 2 ). On exit, POLES(1,*) is an array containing the new singular values obtained from solving the secular equation, and POLES(2,*) is an array containing the poles in the secular equation. Not referenced if ICOMPQ = 0. Secular equation? She's solving. 🧮
DIFL — DOUBLE PRECISION array, dimension ( N ). On exit, DIFL(I) is the distance between the I-th updated (undeflated) singular value and the I-th (undeflated) old singular value. Distance makes the heart grow fonder. 💔
DIFR — DOUBLE PRECISION array, dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and dimension ( K ) if ICOMPQ = 0. On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not defined and will not be referenced. If ICOMPQ = 1, DIFR(1:K,2) is an array containing the normalizing factors for the right singular vector matrix. See DLASD8 for details on DIFL and DIFR. Details in the group chat. 📚
Z — DOUBLE PRECISION array, dimension ( M ). The first elements of this array contain the components of the deflation-adjusted updating row vector. Glow-up row vector. ✨
K — INTEGER. Contains the dimension of the non-deflated matrix. This is the order of the related secular equation. 1 <= K <=N. Keep it in range, bestie. 📏
C — DOUBLE PRECISION. C contains garbage if SQRE =0 and the C-value of a Givens rotation related to the right null space if SQRE = 1. Garbage? We said what we said. 🗑️
S — DOUBLE PRECISION. S contains garbage if SQRE =0 and the S-value of a Givens rotation related to the right null space if SQRE = 1. Same garbage energy. 🗑️
WORK — DOUBLE PRECISION array, dimension ( 4 * M ). Workspace. Do not question her. 💅
IWORK — INTEGER array, dimension ( 3 * N ). Integer workspace. Slay. 🔥
INFO — INTEGER. = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, a singular value did not converge. The receipts. 🧾

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 dlasd7

   subroutine dlasd7 (integer ICOMPQ, integer NL, integer NR, integer SQRE,
       integer K, double precision, dimension( * ) D, double precision,
       dimension( * ) Z, double precision, dimension( * ) ZW, double precision,
       dimension( * ) VF, double precision, dimension( * ) VFW, double
       precision, dimension( * ) VL, double precision, dimension( * ) VLW,
       double precision ALPHA, double precision BETA, double precision,
       dimension( * ) DSIGMA, integer, dimension( * ) IDX, integer, dimension(
       * ) IDXP, integer, dimension( * ) IDXQ, integer, dimension( * ) PERM,
       integer GIVPTR, integer, dimension( ldgcol, * ) GIVCOL, integer LDGCOL,
       double precision, dimension( ldgnum, * ) GIVNUM, integer LDGNUM, double
       precision C, double precision S, integer INFO)

DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc. She's a merger and a deflator. 💅

Purpose:

So like, DLASD7 merges the two sets of singular values 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 singular values 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. DLASD7 is called from DLASD6. Periodt. ✨

Parameters

ICOMPQ — INTEGER. Specifies whether singular vectors are to be computed in compact form, as follows: = 0: compute singular values only. = 1: compute singular vectors of upper bidiagonal matrix in compact form. Compact queen. 💅
NL — INTEGER. The row dimension of the upper block. NL >= 1. Respect. 🦋
NR — INTEGER. The row dimension of the lower block. NR >= 1. Same. 🦋
SQRE — INTEGER. = 0: the lower block is an NR-by-NR square matrix. = 1: the lower block is an NR-by-(NR+1) rectangular matrix. The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns. Math again. 📐
K — INTEGER. Contains the dimension of the non-deflated matrix, this is the order of the related secular equation. 1 <= K <=N. Keep it cute. 📏
D — DOUBLE PRECISION array, dimension ( N ). On entry D contains the singular values of the two submatrices to be combined. On exit D contains the trailing (N-K) updated singular values (those which were deflated) sorted into increasing order. Sorted and updated. ✨
Z — DOUBLE PRECISION array, dimension ( M ). On exit Z contains the updating row vector in the secular equation. Row vector glow. 💖
ZW — DOUBLE PRECISION array, dimension ( M ). Workspace for Z. Supportive bestie. 🫶
VF — DOUBLE PRECISION array, dimension ( M ). On entry, VF(1:NL+1) contains the first components of all right singular vectors of the upper block; and VF(NL+2:M) contains the first components of all right singular vectors of the lower block. On exit, VF contains the first components of all right singular vectors of the bidiagonal matrix. Firsties again. ✨
VFW — DOUBLE PRECISION array, dimension ( M ). Workspace for VF. Supportive bestie. 🫶
VL — DOUBLE PRECISION array, dimension ( M ). On entry, VL(1:NL+1) contains the last components of all right singular vectors of the upper block; and VL(NL+2:M) contains the last components of all right singular vectors of the lower block. On exit, VL contains the last components of all right singular vectors of the bidiagonal matrix. Lasties. 🎀
VLW — DOUBLE PRECISION array, dimension ( M ). Workspace for VL. Supportive bestie. 🫶
ALPHA — DOUBLE PRECISION. Contains the diagonal element associated with the added row. Main character. 👑
BETA — DOUBLE PRECISION. Contains the off-diagonal element associated with the added row. Sidekick. 💜
DSIGMA — DOUBLE PRECISION array, dimension ( N ). Contains a copy of the diagonal elements (K-1 singular values and one zero) in the secular equation. Copy that. 📋
IDX
IDX — INTEGER array, dimension ( N ). This one is holding the receipts: it contains the permutation that got used to sort whatever's in D into ascending order. 📋
IDXP — INTEGER array, dimension ( N ). The placement planner ✨ It contains the permutation used to place the deflated values of D at the end of the array. On output, IDXP(2:K) points to the nondeflated D-values, and IDXP(K+1:N) points to the deflated singular values.
IDXQ — INTEGER array, dimension ( N ). This holds the permutation that separately sorts the two sub-problems in D into ascending order. Heads up tho: entries in the first half of this permutation must first be moved one position backward, and entries in the second half must first have NL+1 added to their values. Do not skip that step, bestie, it is not optional 🙈
PERM — INTEGER array, dimension ( N ). The permutations (from deflation and sorting) that get applied to each singular block. Not referenced if ICOMPQ = 0.
GIVPTR — INTEGER. The number of Givens rotations that took place in this subproblem. Not referenced if ICOMPQ = 0. 🔄
GIVCOL — INTEGER array, dimension ( LDGCOL, 2 ). Every pair of numbers in here points at a pair of columns where a Givens rotation is gonna take place. Not referenced if ICOMPQ = 0.
LDGCOL — INTEGER. The leading dimension of GIVCOL, and it must be at least N.
GIVNUM — DOUBLE PRECISION array, dimension ( LDGNUM, 2 ). Each number in here is the C or S value that gets used in the corresponding Givens rotation. Not referenced if ICOMPQ = 0.
LDGNUM — INTEGER. The leading dimension of GIVNUM, must be at least N.
C — DOUBLE PRECISION. C contains garbage if SQRE = 0, and the C-value of a Givens rotation related to the right null space if SQRE = 1. 🗑️ (yes, actual garbage, the manpage said it, not me)
S — DOUBLE PRECISION. Same energy: S contains garbage if SQRE = 0, and the S-value of a Givens rotation related to the right null space if SQRE = 1.
INFO — INTEGER. The vibe check at the end:
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value. 🚨 illegal argument alert 🚨

🌸 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 `dlasd8` 🦋

subroutine dlasd8 (integer ICOMPQ, integer K, double precision, dimension( * ) D, double precision, dimension( * ) Z, double precision, dimension( * ) VF, double precision, dimension( * ) VL, double precision, dimension( * ) DIFL, double precision, dimension( lddifr, * ) DIFR, integer LDDIFR, double precision, dimension( * ) DSIGMA, double precision, dimension( * ) WORK, integer INFO)

DLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc. Like, it's giving extremely specific mathematician energy and honestly? We respect it. 💅

✨ Purpose ✨

DLASD8 finds the square roots of the roots of the secular equation, as defined by the values in DSIGMA and Z. It makes the appropriate calls to DLASD4, and stores, for each element in D, the distance to its two nearest poles (elements in DSIGMA). It also updates the arrays VF and VL, the first and last components of all the right singular vectors of the original bidiagonal matrix.

DLASD8 is called from DLASD6. 🔗

🌈 Parameters 🌈

ICOMPQ — INTEGER. Says whether singular vectors are to be computed in factored form in the calling routine:
- = 0: Compute singular values only.
- = 1: Compute singular vectors in factored form as well.
K — INTEGER. The number of terms in the rational function to be solved by DLASD4. K >= 1. 🔢
D — DOUBLE PRECISION array, dimension ( K ). On output, D contains the updated singular values. ✨
Z — DOUBLE PRECISION array, dimension ( K ). On entry, the first K elements of this array contain the components of the deflation-adjusted updating row vector. On exit, Z is updated.
VF — DOUBLE PRECISION array, dimension ( K ). On entry, VF contains information passed through DBEDE8. On exit, VF contains the first K components of the first components of all right singular vectors of the bidiagonal matrix.
VL — DOUBLE PRECISION array, dimension ( K ). On entry, VL contains information passed through DBEDE8. On exit, VL contains the first K components of the last components of all right singular vectors of the bidiagonal matrix.
DIFL — DOUBLE PRECISION array, dimension ( K ). On exit, DIFL(I) = D(I) - DSIGMA(I).
DIFR — DOUBLE PRECISION array, dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and dimension ( K ) if ICOMPQ = 0. On exit, DIFR(I,1) = D(I) - DSIGMA(I+1). DIFR(K,1) is not defined and will not be referenced — do not go looking for it, bestie, it is simply not there 🙈. If ICOMPQ = 1, DIFR(1:K,2) is an array containing the normalizing factors for the right singular vector matrix.
LDDIFR — INTEGER. The leading dimension of DIFR, must be at least K.
DSIGMA — DOUBLE PRECISION array, dimension ( K ). On entry, the first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation. On exit, the elements of DSIGMA may be very slightly altered in value. 👀 very slightly, bestie. Barely. A whisper of a change.
WORK — DOUBLE PRECISION array, dimension (3*K).
INFO — INTEGER:
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value. 🚨
- > 0: if INFO = 1, a singular value did not converge. 💀 it did NOT converge

🌸 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 `dlasda` 👑

subroutine dlasda (integer ICOMPQ, integer SMLSIZ, integer N, integer SQRE, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldu, * ) U, integer LDU, double precision, dimension( ldu, * ) VT, integer, dimension( * ) K, double precision, dimension( ldu, * ) DIFL, double precision, dimension( ldu, * ) DIFR, double precision, dimension( ldu, * ) Z, double precision, dimension( ldu, * ) POLES, integer, dimension( * ) GIVPTR, integer, dimension( ldgcol, * ) GIVCOL, integer LDGCOL, integer, dimension( ldgcol, * ) PERM, double precision, dimension( ldu, * ) GIVNUM, double precision, dimension( * ) C, double precision, dimension( * ) S, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. She's the main character of the divide-and-conquer bidiagonal saga, fr fr. 🌟

✨ Purpose ✨

Using a divide and conquer approach, DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal N-by-M matrix B with diagonal D and offdiagonal E, where M = N + SQRE. The algorithm computes the singular values in the SVD B = U * S * VT. The orthogonal matrices U and VT are optionally computed in compact form.

A related subroutine, DLASD0, computes the singular values and the singular vectors in explicit form. 🔗✨

🌈 Parameters 🌈

ICOMPQ — INTEGER. Says whether singular vectors are to be computed in compact form, like this:
- = 0: Compute singular values only.
- = 1: Compute singular vectors of upper bidiagonal matrix in compact form.
SMLSIZ — INTEGER. The maximum size of the subproblems at the bottom of the computation tree. 🌳
N — INTEGER. The row dimension of the upper bidiagonal matrix. This is also the dimension of the main diagonal array D.
SQRE — INTEGER. Specifies the column dimension of the bidiagonal matrix.
- = 0: The bidiagonal matrix has column dimension M = N;
- = 1: The bidiagonal matrix has column dimension M = N + 1.
D — DOUBLE PRECISION array, dimension ( N ). On entry D contains the main diagonal of the bidiagonal matrix. On exit D, if INFO = 0, contains its singular values. ✨
E — DOUBLE PRECISION array, dimension ( M-1 ). Contains the subdiagonal entries of the bidiagonal matrix. On exit, E has been destroyed. 💀 gone. reduced to atoms. do not expect to read it back.
U — DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left singular vector matrices of all subproblems at the bottom level. 💅
LDU — INTEGER, LDU = > N. The leading dimension of arrays U, VT, DIFL, DIFR, POLES, GIVNUM, and Z.
VT — DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT**T contains the right singular vector matrices of all subproblems at the bottom level.
K — INTEGER array, dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th secular equation on the computation tree. 🌳
DIFL — DOUBLE PRECISION array, dimension ( LDU, NLVL ), where NLVL = floor(log_2 (N/SMLSIZ))).
DIFR — okay bestie the DIFR tea is literally right there and we are spilling it the second this part picks back up 💅✨
DIFR — a DOUBLE PRECISION array, and bestie its shape depends on ICOMPQ, so pay attention 💅 dimension ( LDU, 2 * NLVL ) when ICOMPQ = 1, and dimension ( N ) when ICOMPQ = 0. If ICOMPQ = 1, then on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1) record the distances between singular values on the I-th level and the singular values on the (I -1)-th level, and DIFR(1:N, 2 * I ) holds the normalizing factors for the right singular vector matrix. See DLASD8 for the deets. 🦋
Z — a DOUBLE PRECISION array, dimension ( LDU, NLVL ) if ICOMPQ = 1 and dimension ( N ) if ICOMPQ = 0. The first K elements of Z(1, I) contain the components of the deflation-adjusted updating row vector for subproblems on the I-th level. That's it, that's the tea. 🍓
POLES — a DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced at all if ICOMPQ = 0 (she's not invited to that party 🙈). If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and POLES(1, 2*I) contain the new and old singular values involved in the secular equations on the I-th level. ✨
GIVPTR — an INTEGER array, dimension ( N ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records the number of Givens rotations performed on the I-th problem on the computation tree. Counting rotations, so main character of her 🌈
GIVCOL — an INTEGER array, dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I, GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations of Givens rotations performed on the I-th level on the computation tree. Location drop 📍
LDGCOL — an INTEGER, LDGCOL = > N. The leading dimension of arrays GIVCOL and PERM. She's the structural support, lowkey essential 🫶
PERM — an INTEGER array, dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records permutations done on the I-th level of the computation tree. Rearranging, as she does 💅
GIVNUM — a DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I, GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-values of Givens rotations performed on the I-th level on the computation tree. The actual numbers, iconic 🌟
C — a DOUBLE PRECISION array, dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1 and the I-th subproblem is not square, on exit, C( I ) contains the C-value of a Givens rotation related to the right null space of the I-th subproblem. 🎀
S — a DOUBLE PRECISION array, dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1 and the I-th subproblem is not square, on exit, S( I ) contains the S-value of a Givens rotation related to the right null space of the I-th subproblem. C and S, a matched set, we love to see it 💖
WORK — a DOUBLE PRECISION array, dimension (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)). The workspace. Give her the room she asked for. 🧚
IWORK — an INTEGER array, dimension (7*N). Same energy, but integer flavored. 🔢
INFO — an INTEGER, and this is the vibe check right here:
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value. Not her fault, blame the caller 😭
- > 0: if INFO = 1, a singular value did not converge. 💀

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

subroutine dlasdq (character UPLO, integer SQRE, 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)

DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. No cap, that's the whole job description. 💼

✨ Purpose:

DLASDQ computes the singular value decomposition (SVD) of a real (upper or lower) bidiagonal matrix with diagonal D and offdiagonal E, accumulating the transformations if desired. Letting B denote the input bidiagonal matrix, the algorithm computes orthogonal matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose of P). The singular values S are overwritten on D.

The input matrix U is changed to U * Q if desired.

The input matrix VT is changed to P**T * VT if desired.

The input matrix C is changed to Q**T * C if desired.

See "Computing Small Singular Values of Bidiagonal Matrices With Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, LAPACK Working Note #3, for a detailed description of the algorithm. Go read the lore bestie, it's giving foundational 📚

🦋 Parameters

UPLO — UPLO is CHARACTER*1. On entry, UPLO specifies whether the input bidiagonal matrix is upper or lower bidiagonal, and whether it is square are not.
- UPLO = 'U' or 'u' → B is upper bidiagonal.
- UPLO = 'L' or 'l' → B is lower bidiagonal.
SQRE — SQRE is INTEGER.
- = 0: then the input matrix is N-by-N.
- = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and (N+1)-by-N if UPLU = 'L'.

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns.

N — N is INTEGER. On entry, N specifies the number of rows and columns in the matrix. N must be at least 0.
NCVT — NCVT is INTEGER. On entry, NCVT specifies the number of columns of the matrix VT. NCVT must be at least 0.
NRU — NRU is INTEGER. On entry, NRU specifies the number of rows of the matrix U. NRU must be at least 0.
NCC — NCC is INTEGER. On entry, NCC specifies the number of columns of the matrix C. NCC must be at least 0.
D — D is DOUBLE PRECISION array, dimension (N). On entry, D contains the diagonal entries of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in ascending order. Sorted, organized, slaying 💅
E — E is DOUBLE PRECISION array. dimension is (N-1) if SQRE = 0 and N if SQRE = 1. On entry, the entries of E contain the offdiagonal entries of the bidiagonal matrix whose SVD is desired. On normal exit, E will contain 0. If the algorithm does not converge, D and E will contain the diagonal and superdiagonal entries of a bidiagonal matrix orthogonally equivalent to the one given as input.
VT — VT is DOUBLE PRECISION array, dimension (LDVT, NCVT). On entry, contains a matrix which on exit has been premultiplied by P**T, dimension N-by-NCVT if SQRE = 0 and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
LDVT — LDVT is INTEGER. On entry, LDVT specifies the leading dimension of VT as declared in the calling (sub) program. LDVT must be at least 1. If NCVT is nonzero LDVT must also be at least N.
U — U is DOUBLE PRECISION array, dimension (LDU, N). On entry, contains a matrix which on exit has been postmultiplied by Q, dimension NRU-by-N if SQRE = 0 and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
LDU — LDU is INTEGER. On entry, LDU specifies the leading dimension of U as declared in the calling (sub) program. LDU must be at least max( 1, NRU ). Respect the max, bestie 🫡
C — C is DOUBLE PRECISION array, dimension (LDC, NCC). On entry, contains an N-by-NCC matrix which on exit has been premultiplied by Q**T dimension N-by-NCC if SQRE = 0 and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
LDC — LDC is INTEGER. On entry, LDC specifies the leading dimension of C as declared in the calling (sub) program. LDC must be at least 1. If NCC is nonzero, LDC must also be at least N.
WORK — WORK is DOUBLE PRECISION array, dimension (4*N). Workspace. Only referenced if one of NCVT, NRU, or NCC is nonzero, and if N is at least 2. Selective, and honestly? Valid. 💅
INFO — INFO is INTEGER. On exit, a value of 0 indicates a successful exit. If INFO < 0, argument number -INFO is illegal. If INFO > 0, the algorithm did not converge, and INFO specifies how many superdiagonals did not converge. That's your count of things that said "no thanks" 😭

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

subroutine dlasdt (integer N, integer LVL, integer ND, integer, dimension( *
) INODE, integer, dimension( * ) NDIML, integer, dimension( * ) NDIMR,
integer MSUB)

DLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc. She builds the whole family tree 🌳

✨ Purpose:

DLASDT creates a tree of subproblems for bidiagonal divide and conquer.

🦋 Parameters

N — N is INTEGER. On entry, the number of diagonal elements of the bidiagonal matrix.
LVL — LVL is INTEGER. On exit, the number of levels on the computation tree.
ND — ND is INTEGER. On exit, the number of nodes on the tree.
INODE — INODE is INTEGER array, dimension ( N ). On exit, centers of subproblems.
NDIML — NDIML is INTEGER array, dimension ( N ). On exit, row dimensions of left children.
NDIMR
NDIMR — INTEGER array, dimension ( N ). On exit? It's holding the row dimensions of the right children. that's it, that's the whole job. 💖
MSUB — INTEGER. On entry, this is the max row dimension that each subproblem down at the bottom of the tree is allowed to have. like a height restriction but for subproblems.

💖 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 dlaset`

subroutine dlaset (character UPLO, integer M, integer N, double precision
    ALPHA, double precision BETA, double precision, dimension( lda, * ) A,
    integer LDA)

DLASET fills in the off-diagonal elements and the diagonal elements of a matrix with whatever values you hand it. bestie it is literally the "set it and forget it" of matrix init. ✨

🦋 Purpose:

DLASET initializes an m-by-n matrix A to BETA on the diagonal and ALPHA on the offdiagonals. that's the whole vibe.

🌈 Parameters

UPLO — UPLO is CHARACTER*1. Tells the routine which part of the matrix A is getting set, because sometimes you only want half the drama.
- = 'U': the upper triangular part gets set; the strictly lower triangular part of A is not changed. (don't touch her. 🙈)
- = 'L': the lower triangular part gets set; the strictly upper triangular part of A is not changed.
- Otherwise: the whole entire matrix A is set. all of it. no crumbs.
M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
ALPHA — ALPHA is DOUBLE PRECISION. The constant that the offdiagonal elements are being set to.
BETA — BETA is DOUBLE PRECISION. The constant that the diagonal elements are being set to.
A — A is DOUBLE PRECISION array, dimension (LDA,N). On exit, the leading m-by-n submatrix of A is set like this:

if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n,
if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n,
otherwise,     A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j,

and, for all UPLO, A(i,i) = BETA, 1<=i<=min(m,n).

not the i.ne.j sneaking in there to say "diagonal stays out of this" 💅

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

💜 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🧚 `subroutine dlasr`

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

DLASR applies a sequence of plane rotations to a general rectangular matrix. it's giving "spin the whole thing, but make it orthogonal." 🔥

⭐ Purpose:

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

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

A := P*A

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

A := A*P**T

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

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

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

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

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

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

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

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

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

where R(k) appears as a rank-2 modification to the identity matrix in rows and columns k and k+1. iconic behavior tbh

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

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

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

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

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

where R(k) appears in rows and columns k and z. The rotations are performed without ever forming P(k) explicitly. no notes, that's just efficiency, periodt.

💖 Parameters

SIDE — SIDE is CHARACTER*1. Says whether the plane rotation matrix P gets applied to A on the left or the right.
- = 'L': Left, compute A := P*A
- = 'R': Right, compute A:= A*P**T
PIVOT — PIVOT is CHARACTER*1. Specifies the plane for which P(k) is a plane rotation matrix.
- = 'V': Variable pivot, the plane (k,k+1)
- = 'T': Top pivot, the plane (1,k+1)
- = 'B': Bottom pivot, the plane (k,z)
DIRECT — DIRECT is CHARACTER*1. Says whether P is a forward or backward sequence of plane rotations.
- = 'F': Forward, P = P(z-1)*...*P(2)*P(1)
- = 'B': Backward, P = P(1)*P(2)*...*P(z-1)
M — M is INTEGER. The number of rows of the matrix A. If m <= 1, an immediate return is effected. like it just leaves. 🫶
N — N is INTEGER. The number of columns of the matrix A. If n <= 1, an immediate return is effected.
C — C is DOUBLE PRECISION array, dimension

(M-1) if SIDE = 'L'
(N-1) if SIDE = 'R'

The cosines c(k) of the plane rotations.

S — S is DOUBLE PRECISION array, dimension

(M-1) if SIDE = 'L'
(N-1) if SIDE = 'R'

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

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

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

💜 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🎀 `subroutine dlasv2`

subroutine dlasv2 (double precision F, double precision G, double precision
    H, double precision SSMIN, double precision SSMAX, double precision SNR,
    double precision CSR, double precision SNL, double precision CSL)

DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix. tiny matrix, huge main character energy. 👑

🌸 Purpose:

DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix

[  F   G  ]
[  0   H  ].

On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and right singular vectors for abs(SSMAX), giving the decomposition

[ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
[-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].

look at her go 💅✨

🦋 Parameters

F — F is DOUBLE PRECISION. The (1,1) element of the 2-by-2 matrix.
G — G is DOUBLE PRECISION. The (1,2) element of the 2-by-2 matrix.
H — H is DOUBLE PRECISION. The (2,2) element of the 2-by-2 matrix.
SSMIN — SSMIN is DOUBLE PRECISION. abs(SSMIN) is the smaller singular value.
SSMAX — SSMAX is DOUBLE PRECISION. abs(SSMAX) is the larger singular value.
SNL — SNL is DOUBLE PRECISION.
CSL — CSL is DOUBLE PRECISION. The vector (CSL, SNL) is a unit left singular vector for the singular value abs(SSMAX).
SNR — SNR is DOUBLE PRECISION.
CSR — CSR is DOUBLE PRECISION. The vector (CSR, SNR) is a unit right singular vector for the singular value abs(SSMAX).

💜 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 Further Details:

Any input parameter may be aliased with any output parameter. chaotic but valid, we love a flexible queen
Barring over/underflow and assuming a guard digit in subtraction, all output quantities are correct to within a few units in the last place (ulps).
In IEEE arithmetic, the code works correctly if one matrix element is infinite.
Overflow will not occur unless the largest singular value itself overflows or is within a few ulps of overflow. (On machines with partial overflow, like the Cray, overflow may occur if the largest singular value is within a factor of 2 of overflow.)
Underflow is harmless if underflow is gradual. Otherwise, results may correspond to a matrix modified by perturbations of size near the underflow threshold.

bestie if underflow is not gradual, the answers you get back are basically for a slightly different matrix. not a vibe, just a fact. 💀

🦋 `integer function ieeeck (integer ISPEC, real ZERO, real ONE)`

IEEECK — the literal vibe-check for your floating point arithmetic, bestie ✨

Purpose:

IEEECK gets called from ILAENV to verify that Infinity and possibly NaN arithmetic is safe — meaning it won't trap on you 💀

Parameters

ISPEC — ISPEC is INTEGER. Says whether we're testing just infinity arithmetic, or infinity and NaN arithmetic.
- = 0: Verify infinity arithmetic only.
- = 1: Verify infinity and NaN arithmetic.
ZERO — ZERO is REAL. Must contain the value 0.0. It's passed in so the compiler can't be shady and optimize this whole code block away 🙈
ONE — ONE is REAL. Must contain the value 1.0. Same energy — passed in so the compiler doesn't optimize this code away.

RETURN VALUE: INTEGER

= 0: Arithmetic failed to produce the correct answers 😭
= 1: Arithmetic produced the correct answers 🔥

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌸 `integer function iladlc (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA)`

ILADLC scans a matrix for its last non-zero column. Obsessed, honestly 🎀

Purpose:

ILADLC scans A for its last non-zero column. That's it. That's the whole job 💅

Parameters

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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🧚 `integer function iladlr (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA)`

ILADLR scans a matrix for its last non-zero row 🦋

Purpose:

ILADLR scans A for its last non-zero row. Periodt. 👑

Parameters

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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 `integer function ilaenv (integer ISPEC, character( ) NAME, character( ) OPTS, integer N1, integer N2, integer N3, integer N4)`

ILAENV — the main character of tuning parameters, no cap ✨

Purpose:

ILAENV is called from the LAPACK routines to choose problem-dependent parameters for the local environment. Peep ISPEC for a whole description of the parameters 👀

ILAENV returns an INTEGER:

if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC
if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value 💀

This version provides a set of parameters which should give good, but not optimal, performance on many of the currently available computers. Users are encouraged to modify this subroutine to set the tuning parameters for their particular machine using the option and problem size information in the arguments 🛠️

This routine will not function correctly if it is converted to all lower case. Converting it to all upper case is allowed.

Parameters

ISPEC — ISPEC is INTEGER. Specifies the parameter to be returned as the value of ILAENV.
- = 1: the optimal blocksize; if this value is 1, an unblocked algorithm will give the best performance.
- = 2: the minimum block size for which the block routine should be used; if the usable block size is less than this value, an unblocked routine should be used.
- = 3: the crossover point (in a block routine, for N less than this value, an unblocked routine should be used)
- = 4: the number of shifts, used in the nonsymmetric eigenvalue routines (DEPRECATED)
- = 5: the minimum column dimension for blocking to be used; rectangular blocks must have dimension at least k by m, where k is given by ILAENV(2,...) and m by ILAENV(5,...)
- = 6: the crossover point for the SVD (when reducing an m by n matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds this value, a QR factorization is used first to reduce the matrix to a triangular form.)
- = 7: the number of processors
- = 8: the crossover point for the multishift QR method for nonsymmetric eigenvalue problems (DEPRECATED)
- = 9: maximum size of the subproblems at the bottom of the computation tree in the divide-and-conquer algorithm (used by xGELSD and xGESDD)
- =10: ieee infinity and NaN arithmetic can be trusted not to trap
- =11: infinity arithmetic can be trusted not to trap
- 12 <= ISPEC <= 17: xHSEQR or related subroutines, see IPARMQ for detailed explanation
NAME — NAME is CHARACTER*(*). The name of the calling subroutine, in either upper case or lower case.
OPTS — OPTS is CHARACTER*(*). The character options to the subroutine NAME, concatenated into a single character string. For example, UPLO = 'U', TRANS = 'T', and DIAG = 'N' for a triangular routine would be specified as OPTS = 'UTN'.
N1 — N1 is INTEGER
N2 — N2 is INTEGER
N3 — N3 is INTEGER
N4 — N4 is INTEGER. Problem dimensions for the subroutine NAME; these may not all be required.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

The following conventions have been used when calling ILAENV from the LAPACK routines:

OPTS is a concatenation of all of the character options to subroutine NAME, in the same order that they appear in the argument list for NAME, even if they are not used in determining the value of the parameter specified by ISPEC.
The problem dimensions N1, N2, N3, N4 are specified in the order that they appear in the argument list for NAME. N1 is used first, N2 second, and so on, and unused problem dimensions are passed a value of -1.
The parameter value returned by ILAENV is checked for validity in the calling subroutine. For example, ILAENV is used to retrieve the optimal blocksize for STRTRI as follows:

NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
IF( NB.LE.1 ) NB = MAX( 1, N )

💜 `integer function ilaenv2stage (integer ISPEC, character( ) NAME, character( ) OPTS, integer N1, integer N2, integer N3, integer N4)`

ILAENV2STAGE — the two-stage sequel we've all been waiting for 🎀

Purpose:

ILAENV2STAGE is called from the LAPACK routines to choose problem-dependent parameters for the local environment. See ISPEC for a description of the parameters. It sets problem and machine dependent parameters useful for *_2STAGE and related subroutines 🦋

ILAENV2STAGE returns an INTEGER:

if ILAENV2STAGE >= 0: ILAENV2STAGE returns the value of the parameter specified by ISPEC
if ILAENV2STAGE < 0: if ILAENV2STAGE = -k, the k-th argument had an illegal value 💀

This version provides a set of parameters which should give good, but not optimal, performance on many of the currently available computers for the 2-stage solvers. Users are encouraged to modify this subroutine to set the tuning parameters for their particular machine using the option and problem size information in the arguments 🛠️

This routine will not function correctly if it is converted to all lower case. Converting it to all upper case is allowed.

Parameters

ISPEC — ISPEC is INTEGER. Specifies the parameter to be returned as the value of ILAENV2STAGE.
- = 1: the optimal blocksize nb for the reduction to BAND
- = 2: the optimal blocksize ib for the eigenvectors singular vectors update routine
- = 3: The length of the array that store the Housholder representation for the second stage Band to Tridiagonal or Bidiagonal
- = 4: The workspace needed for the routine in input.
- = 5: For future release.
NAME — NAME is CHARACTER*(*). The name of the calling subroutine, in either upper case or lower case.
OPTS — OPTS is CHARACTER*(*). The character options to the subroutine NAME, concatenated into a single character string. For example, UPLO = 'U', TRANS = 'T', and DIAG = 'N' for a triangular routine would be specified as OPTS = 'UTN'.
N1 — N1 is INTEGER
N2 — N2 is INTEGER
N3 — N3 is INTEGER
N4 — N4 is INTEGER. Problem dimensions for the subroutine NAME; these may not all be required.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Nick R. Papior

Further Details:

okay so before we dip into the next function, here are the conventions LAPACK uses when it calls ILAENV2STAGE from the LAPACK routines 💅:

OPTS is a concatenation of all of the character options to subroutine NAME, in the same order that they appear in the argument list for NAME, even if they are not used in determining the value of the parameter specified by ISPEC. periodt.
The problem dimensions N1, N2, N3, N4 are specified in the order that they appear in the argument list for NAME. N1 is used first, N2 second, and so on, and unused problem dimensions are passed a value of -1. no cap.
The parameter value returned by ILAENV2STAGE is checked for validity in the calling subroutine. validity check, slay.

integer function `iparmq`

   integer function iparmq (integer ISPEC, character, dimension( * ) NAME,
       character, dimension( * ) OPTS, integer N, integer ILO, integer IHI,
       integer LWORK)

IPARMQ ✨ the moment ✨

Purpose:

This program sets problem and machine dependent parameters useful for xHSEQR and related subroutines for eigenvalue problems. It is called whenever IPARMQ is called with 12 <= ISPEC <= 16. iykyk.

Parameters

ISPEC — ISPEC is INTEGER. ISPEC specifies which tunable parameter IPARMQ should return. choose your fighter 💖
- ISPEC=12: (INMIN) Matrices of order nmin or less are sent directly to xLAHQR, the implicit double shift QR algorithm. NMIN must be at least 11.
- ISPEC=13: (INWIN) Size of the deflation window. This is best set greater than or equal to the number of simultaneous shifts NS. Larger matrices benefit from larger deflation windows.
- ISPEC=14: (INIBL) Determines when to stop nibbling and invest in an (expensive) multi-shift QR sweep. If the aggressive early deflation subroutine finds LD converged eigenvalues from an order NW deflation window and LD > (NW*NIBBLE)/100, then the next QR sweep is skipped and early deflation is applied immediately to the remaining active diagonal block. Setting IPARMQ(ISPEC=14) = 0 causes TTQRE to skip a multi-shift QR sweep whenever early deflation finds a converged eigenvalue. Setting IPARMQ(ISPEC=14) greater than or equal to 100 prevents TTQRE from skipping a multi-shift QR sweep.
- ISPEC=15: (NSHFTS) The number of simultaneous shifts in a multi-shift QR iteration.
- ISPEC=16: (IACC22) IPARMQ is set to 0, 1 or 2 with the following meanings.
  - 0: During the multi-shift QR/QZ sweep, blocked eigenvalue reordering, blocked Hessenberg-triangular reduction, reflections and/or rotations are not accumulated when updating the far-from-diagonal matrix entries.
  - 1: During the multi-shift QR/QZ sweep, blocked eigenvalue reordering, blocked Hessenberg-triangular reduction, reflections and/or rotations are accumulated, and matrix-matrix multiplication is used to update the far-from-diagonal matrix entries.
  - 2: During the multi-shift QR/QZ sweep, blocked eigenvalue reordering, blocked Hessenberg-triangular reduction, reflections and/or rotations are accumulated, and 2-by-2 block structure is exploited during matrix-matrix multiplies.
  - (If xTRMM is slower than xGEMM, then IPARMQ(ISPEC=16)=1 may be more efficient than IPARMQ(ISPEC=16)=2 despite the greater level of arithmetic work implied by the latter choice.)
- ISPEC=17: (ICOST) An estimate of the relative cost of flops within the near-the-diagonal shift chase compared to flops within the BLAS calls of a QZ sweep.
NAME — NAME is CHARACTER string. Name of the calling subroutine.
OPTS — OPTS is CHARACTER string. This is a concatenation of the string arguments to TTQRE.
N — N is INTEGER. N is the order of the Hessenberg matrix H.
ILO — ILO is INTEGER.
IHI — IHI is INTEGER. It is assumed that H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N.
LWORK — LWORK is INTEGER. The amount of workspace available.

Author

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

Further Details:

Little is known about how best to choose these parameters. It is possible to use different values of the parameters for each of CHSEQR, DHSEQR, SHSEQR and ZHSEQR. it's giving mystery.
It is probably best to choose different parameters for different matrices and different parameters at different times during the iteration, but this has not been implemented --- yet. not yet, bestie.
The best choices of most of the parameters depend in an ill-understood way on the relative execution rate of xLAQR3 and xLAQR5 and on the nature of each particular eigenvalue problem. Experiment may be the only practical way to determine which choices are most effective. vibe check via experiment.
Following is a list of default values supplied by IPARMQ. These defaults may be adjusted in order to attain better performance in any particular computational environment.
IPARMQ(ISPEC=12) The xLAHQR vs xLAQR0 crossover point. Default: 75. (Must be at least 11.)
IPARMQ(ISPEC=13) Recommended deflation window size. This depends on ILO, IHI and NS, the number of simultaneous shifts returned by IPARMQ(ISPEC=15). The default for (IHI-ILO+1) <= 500 is NS. The default for (IHI-ILO+1) > 500 is 3*NS/2.
IPARMQ(ISPEC=14) Nibble crossover point. Default: 14.
IPARMQ(ISPEC=15) Number of simultaneous shifts, NS. a multi-shift QR iteration.

                                   If IHI-ILO+1 is ...

                                   greater than      ...but less    ... the
                                   or equal to ...      than        default is

                                           0               30       NS =   2+
                                          30               60       NS =   4+
                                          60              150       NS =  10
                                         150              590       NS =  **
                                         590             3000       NS =  64
                                        3000             6000       NS = 128
                                        6000             infinity   NS = 256

(+) By default matrices of this order are passed to the implicit double shift routine xLAHQR. See IPARMQ(ISPEC=12) above. These values of NS are used only in case of a rare xLAHQR failure.
(**) The asterisks (**) indicate an ad-hoc function increasing from 10 to 64.
IPARMQ(ISPEC=16) Select structured matrix multiply. (See ISPEC=16 above for details.) Default: 3.
IPARMQ(ISPEC=17) Relative cost heuristic for blocksize selection. Expressed as a percentage. Default: 10.

logical function `lsamen`

   logical function lsamen (integer N, character*( * ) CA, character*( * ) CB)

LSAMEN the case-insensitive bestie

Purpose:

LSAMEN tests if the first N letters of CA are the same as the first N letters of CB, regardless of case. LSAMEN returns .TRUE. if CA and CB are equivalent except for case and .FALSE. otherwise. LSAMEN also returns .FALSE. if LEN( CA ) or LEN( CB ) is less than N. slay comparison.

Parameters

N — N is INTEGER. The number of characters in CA and CB to be compared.
CA — CA is CHARACTER*(*).
CB — CB is CHARACTER*(*). CA and CB specify two character strings of length at least N. Only the first N characters of each string will be accessed.

Author

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

subroutine `scombssq`

   subroutine scombssq (real, dimension( 2 ) V1, real, dimension( 2 ) V2)

SCOMBSSQ adds two scaled sum of squares quantities math moment 🧮

Purpose:

SCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2. That is,

   V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq
                           + V2_scale**2 * V2_sumsq

Parameters

V1 — V1 is REAL array, dimension (2). The first scaled sum. V1(1) = V1_scale, V1(2) = V1_sumsq.
V2 — V2 is REAL array, dimension (2). The second scaled sum. V2(1) = V2_scale, V2(2) = V2_sumsq.

Author

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

logical function `sisnan`

   logical function sisnan (real, intent(in) SIN)

SISNAN tests input for NaN.

Purpose:

SISNAN</p> returns .TRUE. if its argument is NaN, and .FALSE. otherwise. Heads up bestie: this one is slated to be replaced by the Fortran 2003 intrinsic somewhere down the line. 🫶

Parameters

SIN — SIN is REAL. It's the input we're testing for NaN. 💅

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `slabad` (real `SMALL`, real `LARGE`)

SLABAD

Purpose ✨

SLABAD takes as input the values computed by SLAMCH for underflow and overflow, and hands you back the square root of each of those values if the log of LARGE is sufficiently large. No cap.

So like, this subroutine exists to spot machines with a huuuge exponent range — think Crays — and redefine the underflow and overflow limits to be the square roots of whatever SLAMCH computed. Why tho? Because SLAMCH does not compensate for poor arithmetic in the upper half of the exponent range, which is exactly the kind of thing you find on a Cray. 💀

Parameters

SMALL — SMALL is REAL. On entry, it's the underflow threshold as computed by SLAMCH. On exit, if LOG10(LARGE) is sufficiently large, it becomes the square root of SMALL — otherwise it's left untouched, bestie. 🌸
LARGE — LARGE is REAL. On entry, it's the overflow threshold as computed by SLAMCH. On exit, if LOG10(LARGE) is sufficiently large, it becomes the square root of LARGE — otherwise, unchanged. Periodt.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

SLACPY copies all or part of one two-dimensional array to another. It's giving copy-paste but make it linear algebra. ✨

Purpose ✨

SLACPY copies all or part of a two-dimensional matrix A to another matrix B. Simple as that, and honestly? Iconic.

Parameters

UPLO — UPLO is CHARACTER*1. Tells us which part of matrix A gets copied over to B:
- = 'U': Upper triangular part 👑
- = 'L': Lower triangular part 🦋
- Otherwise: All of the matrix A
M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
A — A is REAL array, dimension (LDA,N). The m by n matrix A. If UPLO = 'U', only the upper triangle or trapezoid is accessed; if UPLO = 'L', only the lower triangle or trapezoid is accessed. She's selective, and that's valid. 💅
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). On exit, B = A in the locations specified by UPLO. 🎀
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,M).

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

subroutine `slae2` (real `A`, real `B`, real `C`, real `RT1`, real `RT2`)

SLAE2</p> computes the eigenvalues of a 2-by-2 symmetric matrix. Tiny matrix, big main character energy. 🔥

Purpose ✨

SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix

   [  A   B  ]
   [  B   C  ].

On return, RT1 is the eigenvalue of larger absolute value, and RT2 is the eigenvalue of smaller absolute value. Boom. 🌈

Parameters

A — A is REAL. The (1,1) element of the 2-by-2 matrix.
B — B is REAL. The (1,2) and (2,1) elements of the 2-by-2 matrix.
C — C is REAL. The (2,2) element of the 2-by-2 matrix.
RT1 — RT1 is REAL. The eigenvalue of larger absolute value. 👑
RT2 — RT2 is REAL. The eigenvalue of smaller absolute value.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

RT1 is accurate to a few ulps barring over/underflow. Slay.

RT2 may be inaccurate if there's massive cancellation in the determinant A*C-B*B. To compute RT2 accurately in all cases you'd need higher precision, or correctly rounded, or correctly truncated arithmetic. It's giving "we did our best". 😭

Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds

   underflow_threshold / macheps.

subroutine `slaebz` (integer `IJOB`, integer `NITMAX`, integer `N`, integer `MMAX`, integer `MINP`, integer `NBMIN`, real `ABSTOL`, real `RELTOL`, real `PIVMIN`, real, dimension( * ) `D`, real, dimension( * ) `E`, real, dimension( * ) `E2`, integer, dimension( * ) `NVAL`, real, dimension( mmax, * ) `AB`, real, dimension( * ) `C`, integer `MOUT`, integer, dimension( mmax, * ) `NAB`, real, dimension( * ) `WORK`, integer, dimension( * ) `IWORK`, integer `INFO`)

SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz. She's the whole backstage crew, fr fr. 🫶

Purpose ✨

SLAEBZ contains the iteration loops which compute and use the function N(w), which is the count of eigenvalues of a symmetric tridiagonal matrix T less than or equal to its argument w. It performs a choice of two types of loops:

Note that the intervals are in all cases half-open intervals, i.e., of the form (a,b], which includes b but not a. Inclusive on one side only, she's complicated like that. 💅

To avoid underflow, the matrix should be scaled so that its largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value. And to assure the most accurate computation of small eigenvalues, the matrix should be scaled to be not much smaller than that, either. Basically, right in that Goldilocks zone. 🎀

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966.

Note: the arguments are, in general, not checked for unreasonable values. So like… please be normal about your inputs, bestie. 🙈

Parameters

IJOB — IJOB is INTEGER. Specifies what is to be done:
- = 1: Compute NAB for the initial intervals.
- = 2: Perform bisection iteration to find eigenvalues of T.
- = 3: Perform bisection iteration to invert N(w), i.e., to find a point which has a specified number of eigenvalues of T to its left.
- Other values will cause SLAEBZ to return with INFO=-1. 💀
NITMAX — NITMAX is INTEGER. The maximum number of "levels" of bisection to be performed, i.e., an interval of width W will not be made smaller than 2^(-NITMAX) * W. If not all intervals have converged after NITMAX iterations, then INFO is set to the number of non-converged intervals. She keeps receipts. 🧾
N — N is INTEGER. The dimension n of the tridiagonal matrix T. It must be at least 1.
MMAX — MMAX is INTEGER. The maximum number of intervals. If more than MMAX intervals are generated, then SLAEBZ will quit with INFO=MMAX+1. 😭
MINP — MINP is INTEGER. The initial number of intervals. It may not be greater than MMAX.
NBMIN — NBMIN is INTEGER. The smallest number of intervals that should be processed using a vector loop. If zero, then only the scalar loop will be used. Choose your fighter. ✨
ABSTOL — ABSTOL is REAL. The minimum (absolute) width of an interval. When an interval is narrower than ABSTOL, or than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged. This must be at least zero.
RELTOL — RELTOL is REAL. The minimum relative width of an interval. When an interval is narrower than ABSTOL, or than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged.

Note: this should always be at least radix*machine epsilon.

PIVMIN — PIVMIN is REAL. The minimum absolute value of a "pivot" in the Sturm sequence loop. This must be at least max |e(j)**2|*safe_min and at least safe_min, where safe_min is at least the smallest number that can divide one without overflow.
D — D is REAL array, dimension (N). The diagonal elements of the tridiagonal matrix T. The main characters of the diagonal. 👑
E — E is REAL array, dimension (N). The offdiagonal elements of the tridiagonal matrix T in positions 1 through N-1. E(N) is arbitrary. Do whatever with it, she doesn't care. 🙈
E2 — E2 is REAL array, dimension (N). The squares of the offdiagonal elements of the tridiagonal matrix T. E2(N) is ignored.
NVAL — …

`NVAL` 🌸

NVAL is INTEGER array, dimension (MINP)
- If IJOB=1 or 2, not referenced. 🙈
- If IJOB=3, the desired values of N(w). The elements of NVAL will be reordered to correspond with the intervals in AB. So like, NVAL(j) on output will not, in general, be the same as NVAL(j) on input, but it will correspond with the interval (AB(j,1),AB(j,2)] on output. bestie the array is not where you left it, ok?? 💅

`AB` ✨

AB is REAL array, dimension (MMAX,2)
- The endpoints of the intervals. AB(j,1) is a(j), the left endpoint of the j-th interval, and AB(j,2) is b(j), the right endpoint of the j-th interval. The input intervals will, in general, be modified, split, and reordered by the calculation. it's giving chaos but make it math 🦋

`C` 🎀

C is REAL array, dimension (MMAX)
- If IJOB=1, ignored.
- If IJOB=2, workspace.
- If IJOB=3, then on input C(j) should be initialized to the first search point in the binary search.

`MOUT` 💖

MOUT is INTEGER
- If IJOB=1, the number of eigenvalues in the intervals.
- If IJOB=2 or 3, the number of intervals output.
- If IJOB=3, MOUT will equal MINP.

`NAB` 👑

NAB is INTEGER array, dimension (MMAX,2)
- If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
- If IJOB=2, then on input, NAB(i,j) should be set. It must satisfy the condition:

N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),

which means that in interval i only eigenvalues NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with IJOB=1. On output, NAB(i,j) will contain max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of the input interval that the output interval (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the the input values of NAB(k,1) and NAB(k,2).

If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), unless N(w) > NVAL(i) for all search points w, in which case NAB(i,1) will not be modified, i.e., the output value will be the same as the input value (modulo reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) for all search points w, in which case NAB(i,2) will not be modified. Normally, NAB should be set to some distinctive value(s) before SLAEBZ is called. 😭

`WORK` 🧚

WORK is REAL array, dimension (MMAX)
- Workspace.

`IWORK` 🌈

IWORK is INTEGER array, dimension (MMAX)
- Workspace.

`INFO` ⭐

INFO is INTEGER
- = 0: All intervals converged. slay 💅
- = 1--MMAX: The last INFO intervals did not converge. 💀
- = MMAX+1: More than MMAX intervals were generated.

✨ Author ✨

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💜 Further Details: 💜

This routine is intended to be called only by other LAPACK routines, thus the interface is less user-friendly. It is intended for two purposes:

(a) finding eigenvalues. In this case, SLAEBZ should have one or more initial intervals set up in AB, and SLAEBZ should be called with IJOB=1. This sets up NAB, and also counts the eigenvalues. Intervals with no eigenvalues would usually be thrown out at this point. Also, if not all the eigenvalues in an interval i are desired, NAB(i,1) can be increased or NAB(i,2) decreased. For example, set NAB(i,1)=NAB(i,2)-1 to get the largest eigenvalue. SLAEBZ is then called with IJOB=2 and MMAX no smaller than the value of MOUT returned by the call with IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the tolerance specified by ABSTOL and RELTOL.

(b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). In this case, start with a Gershgorin interval (a,b). Set up AB to contain 2 search intervals, both initially (a,b). One NVAL element should contain f-1 and the other should contain l, while C should contain a and b, resp. NAB(i,1) should be -1 and NAB(i,2) should be N+1, to flag an error if the desired interval does not lie in (a,b). SLAEBZ is then called with IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and w(l-r)=...=w(l+k) are handled similarly. 🔥

🦋 subroutine slaev2 🦋

subroutine slaev2 (real A, real B, real C, real RT1, real RT2, real CS1,
    real SN1)

SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. iconic behavior fr fr 💅

Purpose:

SLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix

[  A   B  ]
[  B   C  ].

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

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

Parameters 🌸

`A`

A is REAL
- The (1,1) element of the 2-by-2 matrix.

`B`

B is REAL
- The (1,2) element and the conjugate of the (2,1) element of the 2-by-2 matrix.

`C`

C is REAL
- The (2,2) element of the 2-by-2 matrix.

`RT1`

RT1 is REAL
- The eigenvalue of larger absolute value.

`RT2`

RT2 is REAL
- The eigenvalue of smaller absolute value.

`CS1`

CS1 is REAL

`SN1`

SN1 is REAL
- The vector (CS1, SN1) is a unit right eigenvector for RT1.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

RT1 is accurate to a few ulps barring over/underflow. 🫶

CS1 and SN1 are accurate to a few ulps barring over/underflow.

Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps.

💜 subroutine slag2d 💜

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

SLAG2D converts a single precision matrix to a double precision matrix. glow-up but for your data ✨

Purpose:

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

Note that while it is possible to overflow while converting from double to single, it is not possible to overflow when converting from single to double. 🙌

This is an auxiliary routine so there is no argument checking. bestie, it trusts you completely, don't abuse that 😭

Parameters 🎀

`M`

M is INTEGER
- The number of lines of the matrix A. M >= 0.

`N`

N is INTEGER
- The number of columns of the matrix A. N >= 0.

`SA`

SA is REAL array, dimension (LDSA,N)
- On entry, the M-by-N coefficient matrix SA.

`LDSA`

LDSA is INTEGER
- The leading dimension of the array SA. LDSA >= max(1,M).

`A`

A is DOUBLE PRECISION array, dimension (LDA,N)
- On exit, the M-by-N coefficient matrix A.

`LDA`

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

`INFO`

INFO is INTEGER
- = 0: successful exit

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 subroutine slagts 🌈

subroutine slagts (integer JOB, integer N, real, dimension( * ) A, real,
    dimension( * ) B, real, dimension( * ) C, real, dimension( * ) D,
    integer, dimension( * ) IN, real, dimension( * ) Y, real TOL, integer
    INFO)

SLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y, where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf. no cap, it just solves it 💅

Purpose:

SLAGTS may be used to solve one of the systems of equations

(T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,

where T is an n by n tridiagonal matrix, for x, following the factorization of (T - lambda*I) as

(T - lambda*I) = P*L*U ,

by routine SLAGTF. The choice of equation to be solved is controlled by the argument JOB, and in each case there is an option to perturb zero or very small diagonal elements of U, this option being intended for use in applications such as inverse iteration. 🔥

Parameters 👑

`JOB`

JOB is INTEGER
- Specifies the job to be performed by SLAGTS as follows:
- = 1: The equations (T - lambda*I)x = y are to be solved, but diagonal elements of U are not to be perturbed.
- = -1: The equations (T - lambda*I)x = y are to be solved and, if overflow would otherwise occur, the diagonal elements of U are to be perturbed. See argument TOL below.
- = 2: The equations (T - lambda*I)**Tx = y are to be solved, but diagonal elements of U are not to be perturbed.
- = -2: The equations (T - lambda*I)**Tx = y are to be solved and, if overflow would otherwise occur, the diagonal elements of U are to be perturbed. See argument TOL below.

`N`

N is INTEGER
- The order of the matrix T.

`A`

A — okay so A is a REAL array, dimension (N). On entry, bestie, A must contain the diagonal elements of U as returned from SLAGTF. slay 💅
B — B is a REAL array, dimension (N-1). On entry, B must contain the first super-diagonal elements of U as returned from SLAGTF. ✨
C — C is a REAL array, dimension (N-1). On entry, C must contain the sub-diagonal elements of L as returned from SLAGTF. 🎀
D — D is a REAL array, dimension (N-2). On entry, D must contain the second super-diagonal elements of U as returned from SLAGTF. 🦋
IN — IN is an INTEGER array, dimension (N). On entry, IN must contain details of the matrix P as returned from SLAGTF. 🌸
Y — Y is a REAL array, dimension (N). On entry, the right hand side vector y. On exit, Y is overwritten by the solution vector x. 💖
TOL — TOL is REAL. On entry, with JOB < 0, TOL should be the minimum perturbation to be made to very small diagonal elements of U. TOL should normally be chosen as about eps*norm(U), where eps is the relative machine precision, but if TOL is supplied as non-positive, then it is reset to eps*max( abs( u(i,j) ) ). If JOB > 0 then TOL is not referenced. On exit, TOL is changed as described above, only if TOL is non-positive on entry. Otherwise TOL is unchanged. 🔥
INFO — INFO is INTEGER.
- = 0: successful exit
- < 0: if INFO = -i, the i-th argument had an illegal value
- > 0: overflow would occur when computing the INFOth element of the solution vector x. This can only occur when JOB is supplied as positive and either means that a diagonal element of U is very small, or that the elements of the right-hand side vector y are very large. 💀

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 `logical function slaisnan (real, intent(in) SIN1, real, intent(in) SIN2)` 🌟

SLAISNAN tests input for NaN by comparing two arguments for inequality. no cap.

✨ Purpose ✨

This routine is not for general use. It exists solely to avoid over-optimization in SISNAN. iykyk

SLAISNAN checks for NaNs by comparing its two arguments for inequality. NaN is the only floating-point value where NaN != NaN returns .TRUE.. To check for NaNs, pass the same variable as both arguments.

🎀 Parameters 🎀

SIN1 — SIN1 is REAL.
SIN2 — SIN2 is REAL. Two numbers to compare for inequality.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🔥 `integer function slaneg (integer N, real, dimension( * ) D, real, dimension( * ) LLD, real SIGMA, real PIVMIN, integer R)` 🔥

SLANEG computes the Sturm count.

✨ Purpose ✨

SLANEG computes the Sturm count, the number of negative pivots encountered while factoring tridiagonal T - sigma I = L D L^T. This implementation works directly on the factors without forming the tridiagonal matrix T. The Sturm count is also the number of eigenvalues of T less than sigma.

This routine is called from SLARRB.

Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 (Tech report version in LAWN 172 with the same title.)

🦋 Parameters 🦋

N — N is INTEGER. The order of the matrix.
D — D is REAL array, dimension (N). The N diagonal elements of the diagonal matrix D.
LLD — LLD is REAL array, dimension (N-1). The (N-1) elements L(i)*L(i)*D(i).
SIGMA — SIGMA is REAL. Shift amount in T - sigma I = L D L^T.
PIVMIN — PIVMIN is REAL. The minimum pivot in the Sturm sequence. May be used when zero pivots are encountered on non-IEEE-754 architectures.
R — R is INTEGER. The twist index for the twisted factorization that is used for the negcount.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 Contributors 🌈

Osni Marques, LBNL/NERSC, USA

Christof Voemel, University of California, Berkeley, USA

Jason Riedy, University of California, Berkeley, USA

💜 `real function slanst (character NORM, integer N, real, dimension( * ) D, real, dimension( * ) E)` 💜

SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.

✨ Purpose ✨

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

🎀 Returns 🎀

SLANST =

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

🌸 Parameters 🌸

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

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ `real function slapy2 (real X, real Y)` ✨

SLAPY2 returns sqrt(x2+y2).

💖 Purpose 💖

SLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary overflow and unnecessary underflow.

🦋 Parameters 🦋

X — X is REAL.
Y — Y is REAL. X and Y specify the values x and y.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌟 `real function slapy3 (real X, real Y, real Z)` 🌟

SLAPY3 returns sqrt(x2+y2+z2).

💖 Purpose 💖

SLAPY3 returns sqrt(x**2+y**2+z**2), taking care not to cause unnecessary overflow and unnecessary underflow.

🦋 Parameters 🦋

X — X is REAL.
Y — Y is REAL.
Z — Z is REAL. X, Y and Z specify the values x, y and z.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🎀 `subroutine slarnv (integer IDIST, integer, dimension( 4 ) ISEED, integer N, real, dimension( * ) X)` 🎀

SLARNV returns a vector of random numbers from a uniform or normal distribution.

✨ Purpose ✨

SLARNV returns a vector of n random real numbers from a uniform or normal distribution.

🌸 Parameters 🌸

IDIST — IDIST is INTEGER. Specifies the distribution of the random numbers:
- = 1: uniform (0,1)
- = 2: uniform (-1,1)
- = 3: normal (0,1)
ISEED — ISEED is INTEGER array, dimension (4). On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. On exit, the seed is updated.
N — N is INTEGER. The number of random numbers to be generated.
X — X is REAL array, dimension (N). The generated random numbers.

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🔥 Further Details 🔥

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

💜 `subroutine slarra (integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) E2, real SPLTOL, real TNRM, integer NSPLIT, integer, dimension( * ) ISPLIT, integer INFO)` 💜

SLARRA computes the splitting points with the specified threshold.

✨ Purpose ✨

Compute the splitting points with threshold SPLTOL. SLARRA sets any "small" off-diagonal elements to zero.

🎀 Parameters 🎀

N — N is INTEGER. The order of the matrix. N > 0.
D — D is REAL array, dimension (N). On entry, the N diagonal elements of the tridiagonal matrix T.
E —
E — REAL array, dimension (N). On entry, the first (N-1) entries hold the subdiagonal elements of the tridiagonal matrix T; E(N) doesn't even need to be set, no pressure 🫶. On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, get set to zero, and the other entries of E are left completely untouched. 💅
E2 — REAL array, dimension (N). On entry, the first (N-1) entries hold the SQUARES of the subdiagonal elements of the tridiagonal matrix T; E2(N) need not be set. On exit, the entries E2( ISPLIT( I ) ), 1 <= I <= NSPLIT, have been set to zero. 🎀
SPLTOL — REAL. The threshold for splitting, bestie. Two criteria can be used: SPLTOL<0 gives you a criterion based on absolute off-diagonal value, and SPLTOL>0 gives you a criterion that preserves relative accuracy. Choose wisely ✨
TNRM — REAL. The norm of the matrix. That's it, that's the whole vibe.
NSPLIT — INTEGER. The number of blocks T splits into. 1 <= NSPLIT <= N.
ISPLIT — INTEGER array, dimension (N). The splitting points, aka where T breaks up into blocks. The first block is rows/columns 1 to ISPLIT(1), the second is rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th one is rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. It's giving very organized bestie energy 🦋
INFO — INTEGER. = 0: successful exit. Periodt. 💖

Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors 💜

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

💖 `subroutine slarrb` 💖

   subroutine slarrb (integer N, real, dimension( * ) D, real, dimension( * )
       LLD, integer IFIRST, integer ILAST, real RTOL1, real RTOL2, integer
       OFFSET, real, dimension( * ) W, real, dimension( * ) WGAP, real,
       dimension( * ) WERR, real, dimension( * ) WORK, integer, dimension( * )
       IWORK, real PIVMIN, real SPDIAM, integer TWIST, integer INFO)

SLARRB provides limited bisection to locate eigenvalues for more accuracy. Bestie it is literally out here doing the detective work on your eigenvalues 🔥

Purpose

Given the relatively robust representation (RRR) L D L^T, SLARRB does "limited" bisection to refine the eigenvalues of L D L^T, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial guesses for these eigenvalues come in through W, and the corresponding estimate of the error in these guesses plus their gaps come in through WERR and WGAP, respectively. During bisection, intervals [left, right] are maintained by storing their mid-points and semi-widths in the arrays W and WERR respectively. So cute, so precise 🎀

Parameters

N — INTEGER. The order of the matrix.
D — REAL array, dimension (N). The N diagonal elements of the diagonal matrix D.
LLD — REAL array, dimension (N-1). The (N-1) elements L(i)*L(i)*D(i).
IFIRST — INTEGER. The index of the first eigenvalue to be computed.
ILAST — INTEGER. The index of the last eigenvalue to be computed.
RTOL1 — REAL. (Yeah bestie, that's literally all the page tells us about this one — no cap 🙈)
RTOL2 — REAL. Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if

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

where GAP is the (estimated) distance to the nearest eigenvalue. Iconic convergence behaviour honestly 💅

OFFSET — INTEGER. Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are the ones being used.
W — REAL array, dimension (N). On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. On output, these estimates are refined. Glow up ✨
WGAP — REAL array, dimension (N-1). On input, the (estimated) gaps between consecutive eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between eigenvalues I and I+1. Note that if IFIRST = ILAST then WGAP(IFIRST-OFFSET) must be set to ZERO — do not skip this bestie, it matters 😭. On output, these gaps are refined.
WERR — REAL array, dimension (N). On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined.
WORK — REAL array, dimension (2*N). Workspace. The little helper doing unpaid labour 🫶
IWORK — INTEGER array, dimension (2*N). Workspace.
PIVMIN — REAL. The minimum pivot in the Sturm sequence.
SPDIAM — REAL. The spectral diameter of the matrix.
TWIST — INTEGER. The twist index for the twisted factorization that is used for the negcount.
- TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
- TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
- TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)

The way this factorization twists?? main character behaviour 👑

INFO — INTEGER. Error flag.

Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors 💜

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

🌈 `subroutine slarrc` 🌈

   subroutine slarrc (character JOBT, integer N, real VL, real VU, real,
       dimension( * ) D, real, dimension( * ) E, real PIVMIN, integer EIGCNT,
       integer LCNT, integer RCNT, integer INFO)

SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix. It's a counter bestie, it counts 🔥

Purpose

Find the number of eigenvalues of the symmetric tridiagonal matrix T that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T if JOBT = 'L'.

Parameters

JOBT — CHARACTER*1.
- = 'T': Compute Sturm count for matrix T.
- = 'L': Compute Sturm count for matrix L D L^T.
N — INTEGER. The order of the matrix. N > 0. Gotta be positive, no zero vibes allowed 💅
VL — REAL. The lower bound for the eigenvalues.
VU — REAL. The upper bound for the eigenvalues.
D — REAL array, dimension (N).
- JOBT = 'T': The N diagonal elements of the tridiagonal matrix T.
- JOBT = 'L': The N diagonal elements of the diagonal matrix D.
E — REAL array, dimension (N).
- JOBT = 'T': The N-1 offdiagonal elements of the matrix T.
- JOBT = 'L': The N-1 offdiagonal elements of the matrix L.
PIVMIN — REAL. The minimum pivot in the Sturm sequence for T.
EIGCNT — INTEGER. The number of eigenvalues of the symmetric tridiagonal matrix T that are in the interval (VL,VU]. There it is, the whole point 🎀
LCNT — INTEGER. (That's all it says bestie, it's a mystery 🙈)
RCNT — INTEGER. The left and right negcounts of the interval.
INFO — INTEGER. (Also just vibing here with no description, iconic honestly ✨)

Author 🌸

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors 💜

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

⭐ `subroutine slarrd` ⭐

   subroutine slarrd (character RANGE, character ORDER, integer N, real VL,
       real VU, integer IL, integer IU, real, dimension( * ) GERS, real RELTOL,
       real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) E2,
       real PIVMIN, integer NSPLIT, integer, dimension( * ) ISPLIT, integer M,
       real, dimension( * ) W, real, dimension( * ) WERR, real WL, real WU,
       integer, dimension( * ) IBLOCK, integer, dimension( * ) INDEXW, real,
       dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)

SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. She's precise, she's powerful, she's the moment 👑

Purpose

SLARRD computes the eigenvalues of a symmetric tridiagonal matrix T to suitable accuracy. This is an auxiliary code to be called from SSTEMR. The user may ask for all eigenvalues, all eigenvalues in the half-open interval (VL, VU], or the IL-th through IU-th eigenvalues.

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. Bestie this is not optional, please scale responsibly 😭

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966. Reading homework, but make it chic 📚

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 —
ORDER — this one's a CHARACTER*1, bestie, and it tells us how the eigenvalues should be sorted. Choices:
- 'B': ("By Block") 💅 the eigenvalues get grouped by their split-off block (peep IBLOCK, ISPLIT) and then ordered from smallest to largest within the block.
- 'E': ("Entire matrix") 👑 the eigenvalues for the whole entire matrix get ordered from smallest to largest.
N — N is INTEGER. It's the order of the tridiagonal matrix T. And it's gotta satisfy N >= 0. Periodt.
VL — VL is REAL. If RANGE='V', this is the lower bound of the interval we're searching for eigenvalues in. 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 REAL. If RANGE='V', this is the upper bound of the interval we're searching for eigenvalues in. 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', it's the index of the smallest eigenvalue to be returned. Condition: 1 <= IL <= IU <= N, if N > 0; and IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.
IU — IU is INTEGER. If RANGE='I', it's the index of the largest eigenvalue to be returned. Same condition: 1 <= IL <= IU <= N, if N > 0; and IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.
GERS — GERS is a REAL array, dimension (2*N). It holds the N Gerschgorin intervals (and the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i))). 🫶
RELTOL — RELTOL is REAL. It's the minimum relative width of an interval. When an interval is narrower than RELTOL times the larger (in magnitude) endpoint, then it's considered sufficiently small, i.e. converged. 💫 Heads up though: this should always be at least radix*machine epsilon.
D — D is a REAL array, dimension (N). The n diagonal elements of the tridiagonal matrix T.
E — E is a REAL array, dimension (N-1). The (n-1) off-diagonal elements of the tridiagonal matrix T.
E2 — E2 is a REAL array, dimension (N-1). The (n-1) squared off-diagonal elements of the tridiagonal matrix T. ✨
PIVMIN — PIVMIN is REAL. The minimum pivot allowed in the Sturm sequence for T. No cap.
NSPLIT — NSPLIT is INTEGER. The number of diagonal blocks in the matrix T, with 1 <= NSPLIT <= N.
ISPLIT — ISPLIT is an INTEGER array, dimension (N). These are the splitting points, where T breaks up into submatrices. The first submatrix is rows/columns 1 to ISPLIT(1), the second is rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th is rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. (Only the first NSPLIT elements actually get used, but since you literally cannot know a priori what value NSPLIT will have, N words must be reserved for ISPLIT.) 🧚
M — M is INTEGER. The actual number of eigenvalues found, 0 <= M <= N. (See also the description of INFO=2,3.)
W — W is a REAL array, dimension (N). On exit, the first M elements of W will contain the eigenvalue approximations. SLARRD computes an interval I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue approximation is given as the interval midpoint W(j)= ( a_j + b_j)/2. The corresponding error is bounded by WERR(j) = abs( a_j - b_j)/2. 🌸
WERR — WERR is a REAL array, dimension (N). The error bound on the corresponding eigenvalue approximation in W.
WL — WL is REAL. 💜
WU — WU is REAL. The interval (WL, WU] contains all the wanted eigenvalues. Here's the tea:
- If RANGE='V', then WL=VL and WU=VU.
- If RANGE='A', then WL and WU are the global Gerschgorin bounds on the spectrum.
- If RANGE='I', then WL and WU are computed by SLAEBZ from the index range specified.
IBLOCK — IBLOCK is an 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 which block (from 1 to the number of blocks) the eigenvalue W(i) belongs to. (And psst — SLARRD may use the remaining N-M elements as workspace.)
INDEXW — INDEXW is an INTEGER array, dimension (N). The indices of the eigenvalues within each block (submatrix); for example, INDEXW(i)= j and IBLOCK(i)=k imply that the i-th eigenvalue W(i) is the j-th eigenvalue in block k. iykyk 🎀
WORK — WORK is a REAL array, dimension (4*N).
IWORK — IWORK is an INTEGER array, dimension (3*N).
INFO — INFO is INTEGER, and this is the receipts:
- = 0: successful exit ✨
- < 0: if INFO = -i, the i-th argument had an illegal value 🙈
- > 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 🛠️

FUDGE — REAL, default = 2. A "fudge factor" to widen the Gerschgorin 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. 💅

🫶 Contributors 🫶

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

👑 Author 👑

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

✨ subroutine slarre ✨

subroutine slarre (character RANGE, integer N, real VL, real VU, integer IL,
    integer IU, real, dimension( * ) D, real, dimension( * ) E, real,
    dimension( * ) E2, real RTOL1, real RTOL2, real SPLTOL, integer NSPLIT,
    integer, dimension( * ) ISPLIT, integer M, real, dimension( * ) W, real,
    dimension( * ) WERR, real, dimension( * ) WGAP, integer, dimension( * )
    IBLOCK, integer, dimension( * ) INDEXW, real, dimension( * ) GERS, real
    PIVMIN, real, dimension( * ) WORK, integer, dimension( * ) IWORK,
    integer INFO)

SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues. 🔥

Purpose:

To find the desired eigenvalues of a given real symmetric tridiagonal matrix T, SLARRE sets any "small" off-diagonal elements to zero, and for each unreduced block T_i, it finds:

(a) a suitable shift at one end of the block's spectrum,
(b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
(c) eigenvalues of each L_i D_i L_i^T.

The representations and eigenvalues found are then used by SSTEMR to compute the eigenvectors of T. The accuracy varies depending on whether bisection is used to find a few eigenvalues or the dqds algorithm (subroutine SLASQ2) to conpute all and then discard any unwanted one. As an added benefit, SLARRE also outputs the n Gerschgorin intervals for the matrices L_i D_i L_i^T. 🌈

Parameters

RANGE — RANGE is CHARACTER*1, and it's giving options:
- '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.
N — N is INTEGER. The order of the matrix. N > 0.
VL — VL is REAL. If RANGE='V', the lower bound for the eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. If RANGE='I' or ='A', SLARRE computes bounds on the desired part of the spectrum.
VU —
VU (REAL) — If RANGE='V', the upper bound for the eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. If RANGE='I' or ='A', SLARRE computes bounds on the desired part of the spectrum. 💖
IL (INTEGER) — If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N. ✨
IU (INTEGER) — If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N. ✨
D (REAL array, dimension (N)) — On entry, the N diagonal elements of the tridiagonal matrix T. On exit, the N diagonal elements of the diagonal matrices D_i. 🦋
E (REAL array, dimension (N)) — On entry, the first (N-1) entries contain the subdiagonal elements of the tridiagonal matrix T; E(N) need not be set. On exit, E contains the subdiagonal elements of the unit bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, contain the base points sigma_i on output. 🦋
E2 (REAL array, dimension (N)) — On entry, the first (N-1) entries contain the SQUARES of the subdiagonal elements of the tridiagonal matrix T; E2(N) need not be set. On exit, the entries E2( ISPLIT( I ) ), 1 <= I <= NSPLIT, have been set to zero. 🦋
RTOL1 (REAL) — It's just a REAL, bestie. No extra deets here. 💅
RTOL2 (REAL) — Parameters for bisection. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ). 💅
SPLTOL (REAL) — The threshold for splitting. ✂️
NSPLIT (INTEGER) — The number of blocks T splits into. 1 <= NSPLIT <= N. 🧩
ISPLIT (INTEGER array, dimension (N)) — The splitting points, at which T breaks up into blocks. The first block consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 🧩
M (INTEGER) — The total number of eigenvalues (of all L_i D_i L_i^T) found. 🔢
W (REAL array, dimension (N)) — The first M elements contain the eigenvalues. The eigenvalues of each of the blocks, L_i D_i L_i^T, are sorted in ascending order (SLARRE may use the remaining N-M elements as workspace). 📈
WERR (REAL array, dimension (N)) — The error bound on the corresponding eigenvalue in W. 📉
WGAP (REAL array, dimension (N)) — The separation from the right neighbor eigenvalue in W. The gap is only with respect to the eigenvalues of the same block as each block has its own representation tree. Exception: at the right end of a block we store the left gap. 🌳
IBLOCK (INTEGER array, dimension (N)) — The indices of the blocks (submatrices) associated with the corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to the first block from the top, =2 if W(i) belongs to the second block, etc. 🧱
INDEXW (INTEGER array, dimension (N)) — The indices of the eigenvalues within each block (submatrix); for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the i-th eigenvalue W(i) is the 10-th eigenvalue in block 2. 🔢
GERS (REAL array, dimension (2*N)) — The N Gerschgorin intervals (the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i))). 🎯
PIVMIN (REAL) — The minimum pivot in the Sturm sequence for T. 🏋️
WORK (REAL array, dimension (6*N)) — Workspace. 🛠️
IWORK (INTEGER array, dimension (5*N)) — Workspace. 🛠️
INFO (INTEGER) —
- = 0: successful exit
- > 0: A problem occurred in SLARRE.
- < 0: One of the called subroutines signaled an internal problem. Needs inspection of the corresponding parameter IINFO for further information.
- =-1: Problem in SLARRD.
- = 2: No base representation could be found in MAXTRY iterations. Increasing MAXTRY and recompilation might be a remedy.
- =-3: Problem in SLARRB when computing the refined root representation for SLASQ2.
- =-4: Problem in SLARRB when preforming bisection on the desired part of the spectrum.
- =-5: Problem in SLASQ2.
- =-6: Problem in SLASQ2.

Author 💖

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

Further Details: ✨

The base representations are required to suffer very little element growth and consequently define all their eigenvalues to high relative accuracy.

Contributors: 🫶

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

✨ `subroutine slarrf` ✨

subroutine slarrf (integer N, real, dimension( * ) D, real, dimension( * )
    L, real, dimension( * ) LD, integer CLSTRT, integer CLEND, real,
    dimension( * ) W, real, dimension( * ) WGAP, real, dimension( * ) WERR,
    real SPDIAM, real CLGAPL, real CLGAPR, real PIVMIN, real SIGMA, real,
    dimension( * ) DPLUS, real, dimension( * ) LPLUS, real, dimension( * )
    WORK, integer INFO)

SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.

Purpose: 🌟

Given the initial representation L D L^T and its cluster of close eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... W( CLEND ), SLARRF finds a new relatively robust representation L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the eigenvalues of L(+) D(+) L(+)^T is relatively isolated.

Parameters 🔍

N (INTEGER) — The order of the matrix (subblock, if the matrix split). 🧮
D (REAL array, dimension (N)) — The N diagonal elements of the diagonal matrix D. 🦋
L (REAL array, dimension (N-1)) — The (N-1) subdiagonal elements of the unit bidiagonal matrix L. 🦋
LD (REAL array, dimension (N-1)) — The (N-1) elements L(i)*D(i). 🦋
CLSTRT (INTEGER) — The index of the first eigenvalue in the cluster. 🎯
CLEND (INTEGER) — The index of the last eigenvalue in the cluster. 🎯
W (REAL array, dimension is >= (CLEND-CLSTRT+1)) — The eigenvalue APPROXIMATIONS of L D L^T in ascending order. W( CLSTRT ) through W( CLEND ) form the cluster of relatively close eigenalues. 📈
WGAP (REAL array, dimension is >= (CLEND-CLSTRT+1)) — The separation from the right neighbor eigenvalue in W. 🌳
WERR (REAL array, dimension is >= (CLEND-CLSTRT+1)) — WERR contain the semiwidth of the uncertainty interval of the corresponding eigenvalue APPROXIMATION in W. 📉
SPDIAM (REAL) — estimate of the spectral diameter obtained from the Gerschgorin intervals. 🎯
CLGAPL (REAL) — No extra deets, just a REAL. 💅
CLGAPR (REAL) — absolute gap on each end of the cluster. Set by the calling routine to protect against shifts too close to eigenvalues outside the cluster. 🛡️
PIVMIN (REAL) — The minimum pivot allowed in the Sturm sequence. 🏋️
SIGMA (REAL) — The shift used to form L(+) D(+) L(+)^T. 🌀
DPLUS (REAL array, dimension (N)) — The N diagonal elements of the diagonal matrix D(+). 🦋
LPLUS (REAL array, dimension (N-1)) — The first (N-1) elements of LPLUS contain the subdiagonal elements of the unit bidiagonal matrix L(+). 🦋
WORK (REAL array, dimension (2*N)) — Workspace. 🛠️
INFO (INTEGER) — Signals processing OK (=0) or failure (=1). 🚦

Author 💖

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

Contributors: 🫶

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

✨ `subroutine slarrj` ✨

subroutine slarrj (integer N, real, dimension( * ) D, real, dimension( * )
    E2, integer IFIRST, integer ILAST, real RTOL, integer OFFSET, real,
    dimension( * ) W, real, dimension( * ) WERR, real, dimension( * ) WORK,
    integer, dimension( * ) IWORK, real PIVMIN, real SPDIAM, integer INFO)

SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.

Purpose: 🌟

okay so given the initial eigenvalue approximations of T, SLARRJ does bisection to refine the eigenvalues of T, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial guesses for these eigenvalues are input in W, the corresponding estimate of the error in these guesses in WERR. During bisection, intervals [left, right] are maintained by storing their mid-points and semi-widths in the arrays W and WERR respectively. periodt 💅

✨ Parameters ✨

N — N is INTEGER. The order of the matrix. slay
D — D is REAL array, dimension (N). The N diagonal elements of T.
E2 — E2 is REAL array, dimension (N-1). The Squares of the (N-1) subdiagonal elements of T.
IFIRST — IFIRST is INTEGER. The index of the first eigenvalue to be computed.
ILAST — ILAST is INTEGER. The index of the last eigenvalue to be computed.
RTOL — RTOL is REAL. Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|).
OFFSET — OFFSET is INTEGER. Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are to be used.
W — W is REAL array, dimension (N). On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. On output, these estimates are refined. ✨
WERR — WERR is REAL array, dimension (N). On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined.
WORK — WORK is REAL array, dimension (2*N). Workspace.
IWORK — IWORK is INTEGER array, dimension (2*N). Workspace.
PIVMIN — PIVMIN is REAL. The minimum pivot in the Sturm sequence for T.
SPDIAM — SPDIAM is REAL. The spectral diameter of T.
INFO — INFO is INTEGER. Error flag.

✨ Author ✨

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

✨ Contributors ✨

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

✨ subroutine slarrk ✨

   subroutine slarrk (integer N, integer IW, real GL, real GU, real, dimension(
       * ) D, real, dimension( * ) E2, real PIVMIN, real RELTOL, real W, real
       WERR, integer INFO)

SLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. bestie, it's giving auxiliary energy ✨

✨ Purpose ✨

SLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. This is an auxiliary code to be called from SSTEMR.

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966.

✨ Parameters ✨

N — N is INTEGER. The order of the tridiagonal matrix T. N >= 0.
IW — IW is INTEGER. The index of the eigenvalues to be returned.
GL — GL is REAL.
GU — GU is REAL. An upper and a lower bound on the eigenvalue.
D — D is REAL array, dimension (N). The n diagonal elements of the tridiagonal matrix T.
E2 — E2 is REAL array, dimension (N-1). The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
PIVMIN — PIVMIN is REAL. The minimum pivot allowed in the Sturm sequence for T.
RELTOL — RELTOL is REAL. The minimum relative width of an interval. When an interval is narrower than RELTOL times the larger (in magnitude) endpoint, then it is considered to be sufficiently small, i.e., converged. Note: this should always be at least radix*machine epsilon.
W — W is REAL.
WERR — WERR is REAL. The error bound on the corresponding eigenvalue approximation in W.
INFO — INFO is INTEGER. = 0: Eigenvalue converged. = -1: Eigenvalue did NOT converge.

✨ Internal Parameters ✨

FUDGE — FUDGE is REAL, default = 2. A "fudge factor" to widen the Gershgorin intervals.

✨ Author ✨

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

✨ subroutine slarrr ✨

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

SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues. vibe check 🦋

✨ Purpose ✨

Perform tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues.

✨ Parameters ✨

N — N is INTEGER. The order of the matrix. N > 0.
D — D is REAL array, dimension (N). The N diagonal elements of the tridiagonal matrix T.
E — E is REAL array, dimension (N). On entry, the first (N-1) entries contain the subdiagonal elements of the tridiagonal matrix T; E(N) is set to ZERO.
INFO — INFO is INTEGER. INFO = 0(default) : the matrix warrants computations preserving relative accuracy. INFO = 1 : the matrix warrants computations guaranteeing only absolute accuracy.

✨ Author ✨

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

✨ Contributors ✨

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

✨ subroutine slartgp ✨

   subroutine slartgp (real F, real G, real CS, real SN, real R)

SLARTGP generates a plane rotation so that the diagonal is nonnegative. iconic ✨

✨ Purpose ✨

SLARTGP generates a plane rotation so that

   [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1.
   [ -SN  CS  ]     [ G ]     [ 0 ]

This is a slower, more accurate version of the Level 1 BLAS routine SROTG, with the following other differences:

F and G are unchanged on return.
If G=0, then CS=(+/-)1 and SN=0.
If F=0 and (G .ne. 0), then CS=0 and SN=(+/-)1.

The sign is chosen so that R >= 0.

✨ Parameters ✨

F — F is REAL. The first component of vector to be rotated.
G — G is REAL. The second component of vector to be rotated.
CS — CS is REAL. The cosine of the rotation.
SN — SN is REAL. The sine of the rotation.
R — R is REAL. The nonzero component of the rotated vector.

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

✨ Author ✨

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

✨ subroutine slaruv ✨

   subroutine slaruv (integer, dimension( 4 ) ISEED, integer N, real,
       dimension( n ) X)

SLARUV returns a vector of n random real numbers from a uniform distribution. random queen 🎀

✨ Purpose ✨

SLARUV returns a vector of n random real numbers from a uniform (0,1) distribution (n <= 128).

This is an auxiliary routine called by SLARNV and CLARNV.

✨ Parameters ✨

ISEED — ISEED is INTEGER array, dimension (4). On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. On exit, the seed is updated.
N — N is INTEGER. The number of random numbers to be generated. N <= 128.
X — X is REAL array, dimension (N). The generated random numbers.

✨ Author ✨

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

✨ Further Details ✨

This routine uses a multiplicative congruential method with modulus 2**48 and multiplier 33952834046453 (see G.S.Fishman, 'Multiplicative congruential random number generators with modulus 2**b: an exhaustive analysis for b = 32 and a partial analysis for b = 48', Math. Comp. 189, pp 331-344, 1990).

48-bit integers are stored in 4 integer array elements with 12 bits per element. Hence the routine is portable across machines with integers of 32 bits or more.

✨ subroutine slas2 ✨

   subroutine slas2 (real F, real G, real H, real SSMIN, real SSMAX)

SLAS2 computes singular values of a 2-by-2 triangular matrix. ate 💅

✨ Purpose ✨

SLAS2 computes the singular values of the 2-by-2 matrix

   [  F   G  ]
   [  0   H  ].

On return, SSMIN is the smaller singular value and SSMAX is the larger singular value.

✨ Parameters ✨

F — F is REAL. The (1,1) element of the 2-by-2 matrix.
G — G is REAL. The (1,2) element of the 2-by-2 matrix.
H — H is REAL. The (2,2) element of the 2-by-2 matrix.
SSMIN
SSMIN — SSMIN is REAL, and it's the smaller singular value. that's it, that's the whole vibe 💅
SSMAX — SSMAX is REAL, and it's the larger singular value. the big one. main character singular value energy 👑

✨ Author ✨

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🌈 Further Details 🌈

okay so bestie, here's the fine print and it's actually important:

Barring over/underflow, all output quantities are correct to within a few units in the last place (ulps), even in the absence of a guard digit in addition/subtraction. 🎀

In IEEE arithmetic, the code works correctly if one matrix element is infinite. no cap, it just handles it.

Underflow is harmless if underflow is gradual. Otherwise, results may correspond to a matrix modified by perturbations of size near the underflow threshold. gradual underflow = fine, abrupt underflow = your results are lowkey not your results anymore 😭

🌸 subroutine `slascl` 🌸

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

SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. it's giving "scale the whole thing, but do it safely" 🫶

Purpose

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

like, the whole point is that it doesn't blow up on you mid-computation. iconic behaviour honestly ⭐

Parameters

TYPE — TYPE is CHARACTER*1. it tells the routine how your matrix is stored, and honestly it's the most important letter you'll type today:
- = 'G': A is a full matrix.
- = 'L': A is a lower triangular matrix.
- = 'U': A is an upper triangular matrix.
- = 'H': A is an upper Hessenberg matrix.
- = 'B': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the lower half stored.
- = 'Q': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the upper half stored.
- = 'Z': A is a band matrix with lower bandwidth KL and upper bandwidth KU. See SGBTRF for storage details.
KL — KL is INTEGER. The lower bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'. so if you're not doing bands, don't stress about it 🙈
KU — KU is INTEGER. The upper bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'. same deal, band gang only 🎀
CFROM — CFROM is REAL. (the denominator of the glow-up, if you will)
CTO — CTO is REAL. The matrix A is multiplied by CTO/CFROM. A(I,J) is computed without over/underflow if the final result CTO*A(I,J)/CFROM can be represented without over/underflow. CFROM must be nonzero.

⚠️ bestie. CFROM must be nonzero. it is literally the denominator. do not do that to yourself 💀

M — M is INTEGER. The number of rows of the matrix A. M >= 0.
N — N is INTEGER. The number of columns of the matrix A. N >= 0.
A — A is REAL array, dimension (LDA,N). The matrix to be multiplied by CTO/CFROM. See TYPE for the storage type.
LDA — LDA is INTEGER. The leading dimension of the array A. here's the vibe check, and yes it's a whole conditional situation:
- If TYPE = 'G', 'L', 'U', 'H', LDA >= max(1,M);
- TYPE = 'B', LDA >= KL+1;
- TYPE = 'Q', LDA >= KU+1;
- TYPE = 'Z', LDA >= 2*KL+KU+1.
INFO — INFO is INTEGER. the receipts 📋
- 0 — successful exit. slay
- <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 `slasd0` 🦋

subroutine slasd0 (integer N, integer SQRE, real, dimension( * ) D, real,
    dimension( * ) E, real, dimension( ldu, * ) U, integer LDU, real,
    dimension( ldvt, * ) VT, integer LDVT, integer SMLSIZ, integer,
    dimension( * ) IWORK, real, dimension( * ) WORK, integer INFO)

SLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc. it's the behind-the-scenes bestie, iykyk 💜

Purpose

Using a divide and conquer approach, SLASD0 computes the singular value decomposition (SVD) of a real upper bidiagonal N-by-M matrix B with diagonal D and offdiagonal E, where M = N + SQRE. The algorithm computes orthogonal matrices U and VT such that B = U * S * VT. The singular values S are overwritten on D.

A related subroutine, SLASDA, computes only the singular values, and optionally, the singular vectors in compact form. that's her low-effort cousin, still valid ✨

Parameters

N — N is INTEGER. On entry, the row dimension of the upper bidiagonal matrix. This is also the dimension of the main diagonal array D.
SQRE — SQRE is INTEGER. Specifies the column dimension of the bidiagonal matrix.
- = 0: The bidiagonal matrix has column dimension M = N;
- = 1: The bidiagonal matrix has column dimension M = N+1;
D — D is REAL array, dimension (N). On entry D contains the main diagonal of the bidiagonal matrix. On exit D, if INFO = 0, contains its singular values. it gets a whole glow-up in place 💅
E — E is REAL array, dimension (M-1). Contains the subdiagonal entries of the bidiagonal matrix. On exit, E has been destroyed. gone. reduced to atoms. do not go looking for it after 💀
U — U is REAL array, dimension (LDU, N). On exit, U contains the left singular vectors.
LDU — LDU is INTEGER. On entry, leading dimension of U.
VT — VT is REAL array, dimension (LDVT, M). On exit, VT**T contains the right singular vectors.
LDVT — LDVT is INTEGER. On entry, leading dimension of VT.
SMLSIZ — SMLSIZ is INTEGER. On entry, maximum size of the subproblems at the bottom of the computation tree. tiny subproblems, tiny tree leaves, big results 🍓
IWORK — IWORK is INTEGER array, dimension (8*N).
WORK — WORK is REAL array, dimension (3*M**2+2*M).
INFO — INFO is INTEGER.
- = 0: successful exit. we love to see it
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, a singular value did not converge. not the vibe. something didn't converge and you gotta deal with it 😭

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 `slasd1` ⭐

subroutine slasd1 (integer NL, integer NR, integer SQRE, real, dimension( *
    ) D, real ALPHA, real BETA, real, dimension( ldu, * ) U, integer LDU,
    real, dimension( ldvt, * ) VT, integer LDVT, integer, dimension( * )
    IDXQ, integer, dimension( * ) IWORK, real, dimension( * ) WORK, integer
    INFO)

SLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc. another one from the sbdsdc cinematic universe 🎬

Purpose

SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0.

A related subroutine SLASD7 handles the case in which the singular values (and the singular vectors in factored form) are desired.

SLASD1 computes the SVD as follows:

              ( D1(in)    0    0       0 )
  B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
              (   0       0   D2(in)   0 )

    = U(out) * ( D(out) 0) * VT(out)

where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros elsewhere; and the entry b is empty if SQRE = 0.

The left singular vectors of the original matrix are stored in U, and the transpose of the right singular vectors are stored in VT, and the singular values are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem when there are multiple singular values or when there are zeros 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 SLASD2.
The second stage consists of calculating the updated singular values. This is done by finding the square roots of the roots of the secular equation via the routine SLASD4 (as called by SLASD3). This routine also calculates the singular vectors of the current problem.
The final stage consists of computing the updated singular vectors directly using the updated singular values. The singular vectors for the current problem are multiplied with the singular vectors from the overall problem.

three stages, one pipeline, zero shortcuts. it's giving methodical. 💖

Parameters

NL — NL is INTEGER. The row dimension of the upper block. NL >= 1.
NR — NR is INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — SQRE is INTEGER.
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has row dimension N = NL + NR + 1, and column dimension M = N + SQRE.

D — D is REAL array, dimension (NL+NR+1). N = NL+NR+1. On entry D(1:NL,1:NL) contains the singular values of the upper block; and D(NL+2:N) contains the singular values of the lower block. On exit D(1:N) contains the singular values of the modified matrix.
ALPHA — ALPHA is REAL. Contains the diagonal element associated with the added row.
BETA — BETA is REAL. Contains the off-diagonal element associated with the added row.
U — U is REAL array, dimension (LDU,N). On entry U(1:NL, 1:NL) contains the left singular vectors of the upper block; U(NL+2:N, NL+2:N) contains the left singular vectors of the lower block. On exit U contains the left singular vectors of the bidiagonal matrix.
LDU
LDU — INTEGER, the leading dimension of the array U. Gotta be LDU >= max( 1, N ). no exceptions bestie
VT — REAL array, dimension (LDVT,M) where M = N + SQRE. On entry, VT(1:NL+1, 1:NL+1)**T is holding the right singular vectors of the upper block, and VT(NL+2:M, NL+2:M)**T is holding the right singular vectors of the lower block. On exit? VT**T contains the right singular vectors of the bidiagonal matrix. glow up complete ✨
LDVT — INTEGER, leading dimension of the array VT. LDVT >= max( 1, M ). respect the dimension or perish
IDXQ — INTEGER array, dimension (N). This one's carrying the permutation that reintegrates the subproblem we just solved back into sorted order — so like, D( IDXQ( I = 1, N ) ) will come out in ascending order. it's giving organized
IWORK — INTEGER array, dimension (4*N). just vibing as workspace
WORK — REAL array, dimension (3*M**2+2*M). big mood, big allocation
INFO — INTEGER
- = 0: successful exit. we ate 💅
- < 0: if INFO = -i, the i-th argument had an illegal value. bestie that argument is NOT valid
- > 0: if INFO = 1, a singular value did not converge. it literally would not converge, embarrassing for it

👑 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 `slasd2`

subroutine slasd2 (integer NL, integer NR, integer SQRE, integer K, real,
    dimension( * ) D, real, dimension( * ) Z, real ALPHA, real BETA, real,
    dimension( ldu, * ) U, integer LDU, real, dimension( ldvt, * ) VT,
    integer LDVT, real, dimension( * ) DSIGMA, real, dimension( ldu2, * )
    U2, integer LDU2, real, dimension( ldvt2, * ) VT2, integer LDVT2,
    integer, dimension( * ) IDXP, integer, dimension( * ) IDX, integer,
    dimension( * ) IDXC, integer, dimension( * ) IDXQ, integer, dimension( *
    ) COLTYP, integer INFO)

SLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc. a lil merge-and-deflate moment

✨ Purpose:

SLASD2 merges the two sets of singular values 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 singular values 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 era

SLASD2 is called from SLASD1.

🎀 Parameters

NL — INTEGER. The row dimension of the upper block. NL >= 1. mandatory, bestie
NR — INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — INTEGER
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns. math is mathing

K — INTEGER. Contains the dimension of the non-deflated matrix. This is the order of the related secular equation. 1 <= K <=N. watch that range fr
D — REAL array, dimension (N). On entry D contains the singular values of the two submatrices to be combined. On exit D contains the trailing (N-K) updated singular values (those which were deflated) sorted into increasing order. deflated but still sorted, we love a comeback
Z — REAL array, dimension (N). On exit Z contains the updating row vector in the secular equation.
ALPHA — REAL. Contains the diagonal element associated with the added row.
BETA — REAL. Contains the off-diagonal element associated with the added row.
U — REAL array, dimension (LDU,N). On entry U contains the left singular vectors of two submatrices in the two square blocks with corners at (1,1), (NL, NL), and (NL+2, NL+2), (N,N). On exit U contains the trailing (N-K) updated left singular vectors (those which were deflated) in its last N-K columns.
LDU — INTEGER. The leading dimension of the array U. LDU >= N.
VT — REAL array, dimension (LDVT,M). On entry VT**T contains the right singular vectors of two submatrices in the two square blocks with corners at (1,1), (NL+1, NL+1), and (NL+2, NL+2), (M,M). On exit VT**T contains the trailing (N-K) updated right singular vectors (those which were deflated) in its last N-K columns. In case SQRE =1, the last row of VT spans the right null space. spooky null space behavior ✨
LDVT — INTEGER. The leading dimension of the array VT. LDVT >= M.
DSIGMA — REAL array, dimension (N). Contains a copy of the diagonal elements (K-1 singular values and one zero) in the secular equation.
U2 — REAL array, dimension (LDU2,N). Contains a copy of the first K-1 left singular vectors which will be used by SLASD3 in a matrix multiply (SGEMM) to solve for the new left singular vectors. U2 is arranged into four blocks. The first block contains a column with 1 at NL+1 and zero everywhere else; the second block contains non-zero entries only at and above NL; the third contains non-zero entries only below NL+1; and the fourth is dense. four blocks, one slay
LDU2 — INTEGER. The leading dimension of the array U2. LDU2 >= N.
VT2 — REAL array, dimension (LDVT2,N). VT2**T contains a copy of the first K right singular vectors which will be used by SLASD3 in a matrix multiply (SGEMM) to solve for the new right singular vectors. VT2 is arranged into three blocks. The first block contains a row that corresponds to the special 0 diagonal element in SIGMA; the second block contains non-zeros only at and before NL +1; the third block contains non-zeros only at and after NL +2. three blocks this time, iconic
LDVT2 — INTEGER. The leading dimension of the array VT2. LDVT2 >= M.
IDXP — INTEGER array, dimension (N). This will contain the permutation used to place deflated values of D at the end of the array. On output IDXP(2:K) points to the nondeflated D-values and IDXP(K+1:N) points to the deflated singular values.
IDX — INTEGER array, dimension (N). This will contain the permutation used to sort the contents of D into ascending order.
IDXC — INTEGER array, dimension (N). This will contain the permutation used to arrange the columns of the deflated U matrix into three groups: the first group contains non-zero entries only at and above NL, the second contains non-zero entries only below NL+2, and the third is dense.
IDXQ — INTEGER array, dimension (N). This contains the permutation which separately sorts the two sub-problems in D into ascending order. Note that entries in the first half of this permutation must first be moved one position backward; and entries in the second half must first have NL+1 added to their values. heads up, both halves need a lil pre-edit
COLTYP — INTEGER array, dimension (N). As workspace, this will contain a label which will indicate which of the following types a column in the U2 matrix or a row in the VT2 matrix is:
- 1 : non-zero in the upper half only
- 2 : non-zero in the lower half only
- 3 : dense
- 4 : deflated

On exit, it is an array of dimension 4, with COLTYP(I) being the dimension of the I-th type columns. labels on labels

INFO — INTEGER
- = 0: successful exit. clean 💖
- < 0: if INFO = -i, the i-th argument had an illegal value. girl what argument was that

👑 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 `slasd3`

subroutine slasd3 (integer NL, integer NR, integer SQRE, integer K, real,
    dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, real,
    dimension( * ) DSIGMA, real, dimension( ldu, * ) U, integer LDU, real,
    dimension( ldu2, * ) U2, integer LDU2, real, dimension( ldvt, * ) VT,
    integer LDVT, real, dimension( ldvt2, * ) VT2, integer LDVT2, integer,
    dimension( * ) IDXC, integer, dimension( * ) CTOT, real, dimension( * )
    Z, integer INFO)

SLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc. roots of roots, very meta of it

✨ Purpose:

SLASD3 finds all the square roots of the roots of the secular equation, as defined by the values in D and Z. It makes the appropriate calls to SLASD4 and then updates the singular vectors by matrix multiplication.

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 XMP, Cray YMP, Cray C 90, or Cray 2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. conceivably fail?? on hex/decimal machines with no guard digit. dramatic.

SLASD3 is called from SLASD1.

🎀 Parameters

NL — INTEGER. The row dimension of the upper block. NL >= 1.
NR — INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — INTEGER
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns. same shape lore as before

K
K — K is INTEGER. The size of the secular equation, 1 =< K = < N. Yes bestie, that's a real double inequality, that's just how it is 💅
D — D is REAL array, dimension(K). On exit? It holds the square roots of the roots of the secular equation, in ascending order. Glow up from roots to square roots, we love an ascending arc ✨
Q — Q is REAL array, dimension (LDQ,K). Just sitting there, being an array, no notes.
LDQ — LDQ is INTEGER. The leading dimension of the array Q. LDQ >= K. It's the boundary, respect it 🫶
DSIGMA — DSIGMA is 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. The drama of it all, honestly iconic.
U — U is REAL array, dimension (LDU, N). The last N - K columns of this matrix contain the deflated left singular vectors.
LDU — LDU is INTEGER. The leading dimension of the array U. LDU >= N.
U2 — U2 is REAL array, dimension (LDU2, N). The first K columns of this matrix contain the non-deflated left singular vectors for the split problem. Not deflated. Hydrated, even 💧
LDU2 — LDU2 is INTEGER. The leading dimension of the array U2. LDU2 >= N.
VT — VT is REAL array, dimension (LDVT, M). The last M - K columns of VT**T contain the deflated right singular vectors. (Yes, VT is holding VT**T's stuff. It's giving transposed energy 🔄)
LDVT — LDVT is INTEGER. The leading dimension of the array VT. LDVT >= N.
VT2 — VT2 is REAL array, dimension (LDVT2, N). The first K columns of VT2**T contain the non-deflated right singular vectors for the split problem.
LDVT2 — LDVT2 is INTEGER. The leading dimension of the array VT2. LDVT2 >= N.
IDXC — IDXC is INTEGER array, dimension (N). The permutation used to arrange the columns of U (and rows of VT) into three groups: the first group contains non-zero entries only at and above (or before) NL +1; the second contains non-zero entries only at and below (or after) NL+2; and the third is dense. The first column of U and the row of VT are treated separately, however. Three groups, one iconique little exception.

The rows of the singular vectors found by SLASD4 must be likewise permuted before the matrix multiplies can take place. No cap — permute first, multiply later, that's the law 📋

CTOT — CTOT is INTEGER array, dimension (4). A count of the total number of the various types of columns in U (or rows in VT), as described in IDXC. The fourth column type is any column which has been deflated. A little census of the vibes 📊
Z — Z is REAL array, dimension (K). The first K elements of this array contain the components of the deflation-adjusted updating row vector. Deflation-adjusted = ✨self-care✨
INFO — INFO is INTEGER. The receipts:
- = 0: successful exit. 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value. Rude argument alert.
- > 0: if INFO = 1, a singular value did not converge. She did NOT converge, bestie. It happens.

👑 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 slasd4 🦋

subroutine slasd4 (integer N, integer I, real, dimension( * ) D, real,
    dimension( * ) Z, real, dimension( * ) DELTA, real RHO, real SIGMA,
    real, dimension( * ) WORK, integer INFO)

SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by sbdsdc. Very specific, very slay.

🌸 Purpose:

This subroutine computes the square root of the I-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix whose entries are given as the squares of the corresponding entries in the array d, and that

       0 <= 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 ) * 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 approximations, the audacity, the craft 👑

🎀 Parameters

N — N is INTEGER. The length of all arrays. All of them. Every single one 🫡
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, 0 <= D(I) < D(J) for I < J. Sorted, as they should be 💅
Z — Z is REAL array, dimension ( N ). The components of the updating vector.
DELTA — DELTA is REAL array, dimension ( N ). If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th component. If N = 1, then DELTA(1) = 1. The vector DELTA contains the information necessary to construct the (singular) eigenvectors. It's giving... the blueprint 📐
RHO — RHO is REAL. The scalar in the symmetric updating formula.
SIGMA — SIGMA is REAL. The computed sigma_I, the I-th updated eigenvalue. She's computed. She's sigma. 💖
WORK — WORK is REAL array, dimension ( N ). If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th component. If N = 1, then WORK( 1 ) = 1. Note the plus. DELTA was minus. They're a duo 🤝
INFO — INFO is INTEGER.
- = 0: successful exit ✅
- > 0: if INFO = 1, the updating process failed. She tried. She failed. We move on. 😭

🔍 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. Three poles. Say it with me 🔥

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 slasd5 🌈

subroutine slasd5 (integer I, real, dimension( 2 ) D, real, dimension( 2 )
    Z, real, dimension( 2 ) DELTA, real RHO, real DSIGMA, real, dimension( 2
    ) WORK)

SLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. Used by sbdsdc. Two by two. Tiny. Adorable. Powerful.

🌸 Purpose:

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

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

The diagonal entries in the array D are assumed to satisfy

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

We also assume RHO > 0 and that the Euclidean norm of the vector Z is one. Norm of one. Always norm of one. It's the house rule 🏠✨

🎀 Parameters

I — I is INTEGER. The index of the eigenvalue to be computed. I = 1 or I = 2. That's the whole menu, bestie 🍓
D — D is REAL array, dimension (2). The original eigenvalues. We assume 0 <= D(1) < D(2).
Z — Z is REAL array, dimension (2). The components of the updating vector.
DELTA — DELTA is REAL array, dimension (2). Contains (D(j) - sigma_I) in its j-th component. The vector DELTA contains the information necessary to construct the eigenvectors.
RHO — RHO is REAL. The scalar in the symmetric updating formula.
DSIGMA — DSIGMA is REAL. The computed sigma_I, the I-th updated eigenvalue. Look at her, all computed and everything 💖
WORK — WORK is REAL array, dimension (2). WORK contains (D(j) + sigma_I) in its j-th component.

👑 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 slasd6 💜

subroutine slasd6 (integer ICOMPQ, integer NL, integer NR, integer SQRE,
    real, dimension( * ) D, real, dimension( * ) VF, real, dimension( * )
    VL, real ALPHA, real BETA, integer, dimension( * ) IDXQ, integer,
    dimension( * ) PERM, integer GIVPTR, integer, dimension( ldgcol, * )
    GIVCOL, integer LDGCOL, real, dimension( ldgnum, * ) GIVNUM, integer
    LDGNUM, real, dimension( ldgnum, * ) POLES, real, dimension( * ) DIFL,
    real, dimension( * ) DIFR, real, dimension( * ) Z, integer K, real C,
    real S, real, dimension( * ) WORK, integer, dimension( * ) IWORK,
    integer INFO)

SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc. Merging. Appending. SVD. The full package 🎁

🌸 Purpose:

SLASD6 computes the SVD of an updated upper bidiagonal matrix B obtained by merging two smaller ones by appending a row. This routine is used only for the problem which requires all singular values and optionally singular vector matrices in factored form. B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. A related subroutine, SLASD1, handles the case in which all singular values and singular vectors of the bidiagonal matrix are desired. (Shoutout SLASD1, different routine, different vibe, same family 🫶)

SLASD6 computes the SVD as follows:

              ( D1(in)    0    0       0 )
  B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
              (   0       0   D2(in)   0 )

    = U(out) * ( D(out) 0) * VT(out)

The singular values of B can be computed using D1, D2, the first components of all the right singular vectors of the lower block, and the last components of all the right singular vectors of the upper block. These components are stored and updated in VF and VL, respectively, in SLASD6. Hence U and VT are not explicitly referenced. Not even referenced. They're just vibing offstage, unbothered 🙈

The singular values are stored in D. The algorithm consists of two stages:

okay so stage one is literally all about deflating the problem size when you've got multiple singular values or a zero just vibing in the Z vector. Every time that happens, the dimension of the secular equation problem gets reduced by one. This whole stage is handled by the routine SLASD7. 💅

then stage two is where we calculate the updated singular values. We do that by finding the roots of the secular equation using SLASD4 (which is called by SLASD8). That routine also updates VF and VL and computes the distances between the new singular values and the old ones. slay ✨

SLASD6 gets called from SLASDA. periodt 🦋

✨ Parameters ✨

okay bestie, here's the full parameter breakdown:

ICOMPQ — INTEGER. Tells us if singular vectors should be computed in factored form:
- = 0: compute singular values only.
- = 1: compute singular vectors in factored form as well.
NL — INTEGER. The row dimension of the upper block. NL >= 1.
NR — INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — INTEGER.
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has row dimension N = NL + NR + 1, and column dimension M = N + SQRE.

D — REAL array, dimension (NL+NR+1). On entry, D(1:NL,1:NL) contains the singular values of the upper block, and D(NL+2:N) contains the singular values of the lower block. On exit, D(1:N) contains the singular values of the modified matrix.
VF — REAL array, dimension (M). On entry, VF(1:NL+1) contains the first components of all right singular vectors of the upper block; and VF(NL+2:M) contains the first components of all right singular vectors of the lower block. On exit, VF contains the first components of all right singular vectors of the bidiagonal matrix.
VL — REAL array, dimension (M). On entry, VL(1:NL+1) contains the last components of all right singular vectors of the upper block; and VL(NL+2:M) contains the last components of all right singular vectors of the lower block. On exit, VL contains the last components of all right singular vectors of the bidiagonal matrix.
ALPHA — REAL. Contains the diagonal element associated with the added row.
BETA — REAL. Contains the off-diagonal element associated with the added row.
IDXQ — INTEGER array, dimension (N). This contains the permutation which will reintegrate the subproblem just solved back into sorted order, i.e. D( IDXQ( I = 1, N ) ) will be in ascending order.
PERM — INTEGER array, dimension (N). The permutations (from deflation and sorting) to be applied to each block. Not referenced if ICOMPQ = 0.
GIVPTR — INTEGER. The number of Givens rotations which took place in this subproblem. Not referenced if ICOMPQ = 0.
GIVCOL — INTEGER array, dimension ( LDGCOL, 2 ). Each pair of numbers indicates a pair of columns to take place in a Givens rotation. Not referenced if ICOMPQ = 0.
LDGCOL — INTEGER. Leading dimension of GIVCOL, must be at least N.
GIVNUM — REAL array, dimension ( LDGNUM, 2 ). Each number indicates the C or S value to be used in the corresponding Givens rotation. Not referenced if ICOMPQ = 0.
LDGNUM — INTEGER. The leading dimension of GIVNUM and POLES, must be at least N.
POLES — REAL array, dimension ( LDGNUM, 2 ). On exit, POLES(1,*) is an array containing the new singular values obtained from solving the secular equation, and POLES(2,*) is an array containing the poles in the secular equation. Not referenced if ICOMPQ = 0.
DIFL — REAL array, dimension (N). On exit, DIFL(I) is the distance between the I-th updated (undeflated) singular value and the I-th (undeflated) old singular value.
DIFR — REAL array, dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and dimension ( K ) if ICOMPQ = 0. On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not defined and will not be referenced.

If ICOMPQ = 1, DIFR(1:K,2) is an array containing the normalizing factors for the right singular vector matrix.

See SLASD8 for details on DIFL and DIFR.

Z — REAL array, dimension (M). The first elements of this array contain the components of the deflation-adjusted updating row vector.
K — INTEGER. Contains the dimension of the non-deflated matrix. This is the order of the related secular equation. 1 <= K <=N.
C — REAL. C contains garbage if SQRE =0 and the C-value of a Givens rotation related to the right null space if SQRE = 1. garbage means don't touch it, bestie 💀
S — REAL. S contains garbage if SQRE =0 and the S-value of a Givens rotation related to the right null space if SQRE = 1. same energy 💀
WORK — REAL array, dimension ( 4 * M ).
IWORK — INTEGER array, dimension ( 3 * N ).
INFO — INTEGER.
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, a singular value did not converge.

✨ Author ✨

Shoutout to the legends:

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

✨ Contributors ✨

Big ups to:

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

✨ subroutine slasd7 ✨

   subroutine slasd7 (integer ICOMPQ, integer NL, integer NR, integer SQRE,
       integer K, real, dimension( * ) D, real, dimension( * ) Z, real,
       dimension( * ) ZW, real, dimension( * ) VF, real, dimension( * ) VFW,
       real, dimension( * ) VL, real, dimension( * ) VLW, real ALPHA, real
       BETA, real, dimension( * ) DSIGMA, integer, dimension( * ) IDX, integer,
       dimension( * ) IDXP, integer, dimension( * ) IDXQ, integer, dimension( *
       ) PERM, integer GIVPTR, integer, dimension( ldgcol, * ) GIVCOL, integer
       LDGCOL, real, dimension( ldgnum, * ) GIVNUM, integer LDGNUM, real C,
       real S, integer INFO)

SLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc. 🔥

✨ Purpose ✨

SLASD7 merges the two sets of singular values 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 singular values 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.

SLASD7 is called from SLASD6. it's giving teamwork 🤝

✨ Parameters ✨

ICOMPQ — INTEGER. Specifies whether singular vectors are to be computed in compact form, as follows:
- = 0: Compute singular values only.
- = 1: Compute singular vectors of upper bidiagonal matrix in compact form.
NL — INTEGER. The row dimension of the upper block. NL >= 1.
NR — INTEGER. The row dimension of the lower block. NR >= 1.
SQRE — INTEGER.
- = 0: the lower block is an NR-by-NR square matrix.
- = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns.

K — INTEGER. Contains the dimension of the non-deflated matrix, this is the order of the related secular equation. 1 <= K <=N.
D — REAL array, dimension (N). On entry D contains the singular values of the two submatrices to be combined. On exit D contains the trailing (N-K) updated singular values (those which were deflated) sorted into increasing order.
Z — REAL array, dimension (M). On exit Z contains the updating row vector in the secular equation.
ZW — REAL array, dimension (M). Workspace for Z.
VF — REAL array, dimension (M). On entry, VF(1:NL+1) contains the first components of all right singular vectors of the upper block; and VF(NL+2:M) contains the first components of all right singular vectors of the lower block. On exit, VF contains the first components of all right singular vectors of the bidiagonal matrix.
VFW — REAL array, dimension (M). Workspace for VF.
VL — REAL array, dimension (M). On entry, VL(1:NL+1) contains the last components of all right singular vectors of the upper block; and VL(NL+2:M) contains the last components of all right singular vectors of the lower block. On exit, VL contains the last components of all right singular vectors of the bidiagonal matrix.
VLW — REAL array, dimension (M). Workspace for VL.
ALPHA — REAL. Contains the diagonal element associated with the added row.
BETA —
BETA — it's REAL, bestie. Holds the off-diagonal element that comes with the added row. that's it, that's the whole job 💅
DSIGMA — a REAL array, dimension ( N ). Contains a copy of the diagonal elements — that's K-1 singular values plus one (1) zero — living inside the secular equation. iconic behavior honestly.
IDX — INTEGER array, dimension ( N ). This one's gonna hold the permutation used to sort what's inside D into ascending order. sorting is self-care ✨
IDXP — INTEGER array, dimension ( N ). Contains the permutation used to plop the deflated values of D at the end of the array. On output, IDXP(2:K) points to the nondeflated D-values and IDXP(K+1:N) points to the deflated singular values. two neighborhoods, one array, we love to see it 🦋
IDXQ — INTEGER array, dimension ( N ). Contains the permutation that separately sorts the two sub-problems in D into ascending order. ⚠️ heads up tho: entries in the first half of this permutation must first be moved one position backward, and entries in the second half must first have NL+1 added to their values. do not skip that step, it's giving "I forgot the offset" energy and that's a whole mess 😭
PERM — INTEGER array, dimension ( N ). The permutations (from deflation and sorting, yes both) to be applied to each singular block. Not referenced if ICOMPQ = 0.
GIVPTR — INTEGER. The number of Givens rotations which took place in this subproblem. Not referenced if ICOMPQ = 0. it's literally just counting the rotations, no cap 🌀
GIVCOL — INTEGER array, dimension ( LDGCOL, 2 ). Each pair of numbers indicates a pair of columns to take place in a Givens rotation. Not referenced if ICOMPQ = 0.
LDGCOL — INTEGER. The leading dimension of GIVCOL, and it must be at least N. no exceptions bestie.
GIVNUM — REAL array, dimension ( LDGNUM, 2 ). Each number indicates the C or S value to be used in the corresponding Givens rotation. Not referenced if ICOMPQ = 0.
LDGNUM — INTEGER. The leading dimension of GIVNUM, must be at least N. same rules, same energy 🎀
C — REAL. C contains garbage if SQRE = 0, and the C-value of a Givens rotation related to the right null space if SQRE = 1. yes, literal garbage, don't go reading it 💀
S — REAL. S contains garbage if SQRE = 0, and the S-value of a Givens rotation related to the right null space if SQRE = 1. twin of C, same warning 🙈
INFO — INTEGER. the vibe check on how it all went:
- = 0: successful exit. slay.
- < 0: if INFO = -i, the i-th argument had an illegal value. somebody passed something illegal, rude.

🌸 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 slasd8 🧚

subroutine slasd8 (integer ICOMPQ, integer K, real, dimension( * ) D, real,
    dimension( * ) Z, real, dimension( * ) VF, real, dimension( * ) VL,
    real, dimension( * ) DIFL, real, dimension( lddifr, * ) DIFR, integer
    LDDIFR, real, dimension( * ) DSIGMA, real, dimension( * ) WORK, integer
    INFO)

SLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc. it's giving "I do one very specific thing and I do it flawlessly" 👑

Purpose:

SLASD8 finds the square roots of the roots of the secular equation, as defined by the values in DSIGMA and Z. It makes the appropriate calls to SLASD4, and stores, for each element in D, the distance to its two nearest poles (elements in DSIGMA). It also updates the arrays VF and VL, the first and last components of all the right singular vectors of the original bidiagonal matrix.

SLASD8 is called from SLASD6. family tree moment 🌳

Parameters

ICOMPQ — INTEGER. Specifies whether singular vectors are to be computed in factored form in the calling routine:
- = 0: Compute singular values only.
- = 1: Compute singular vectors in factored form as well. the full package, periodt ✨
K — INTEGER. The number of terms in the rational function to be solved by SLASD4. Constraint: K >= 1.
D — REAL array, dimension ( K ). On output, D contains the updated singular values. glow-up complete 💖
Z — REAL array, dimension ( K ). On entry, the first K elements of this array contain the components of the deflation-adjusted updating row vector. On exit, Z is updated.
VF — REAL array, dimension ( K ). On entry, VF contains information passed through DBEDE8. On exit, VF contains the first K components of the first components of all right singular vectors of the bidiagonal matrix.
VL — REAL array, dimension ( K ). On entry, VL contains information passed through DBEDE8. On exit, VL contains the first K components of the last components of all right singular vectors of the bidiagonal matrix. VF's sibling, obviously 🦋
DIFL — REAL array, dimension ( K ). On exit, DIFL(I) = D(I) - DSIGMA(I). math is math, bestie.
DIFR — REAL array, dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and dimension ( K ) if ICOMPQ = 0. On exit, DIFR(I,1) = D(I) - DSIGMA(I+1). DIFR(K,1) is not defined and will not be referenced — so don't come for it, it simply does not exist 💀

If ICOMPQ = 1, DIFR(1:K,2) is an array containing the normalizing factors for the right singular vector matrix.

LDDIFR — INTEGER. The leading dimension of DIFR, must be at least K.
DSIGMA — REAL array, dimension ( K ). On entry, the first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation. On exit, the elements of DSIGMA may be very slightly altered in value. very slightly. like a whisper of change 🌸
WORK — REAL array, dimension (3*K)
INFO — INTEGER:
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, a singular value did not converge. 😭 it tried, it really did.

🌸 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 slasda 🎀

subroutine slasda (integer ICOMPQ, integer SMLSIZ, integer N, integer SQRE,
    real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldu, *
    ) U, integer LDU, real, dimension( ldu, * ) VT, integer, dimension( * )
    K, real, dimension( ldu, * ) DIFL, real, dimension( ldu, * ) DIFR, real,
    dimension( ldu, * ) Z, real, dimension( ldu, * ) POLES, integer,
    dimension( * ) GIVPTR, integer, dimension( ldgcol, * ) GIVCOL, integer
    LDGCOL, integer, dimension( ldgcol, * ) PERM, real, dimension( ldu, * )
    GIVNUM, real, dimension( * ) C, real, dimension( * ) S, real, dimension(
    * ) WORK, integer, dimension( * ) IWORK, integer INFO)

SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. main character algorithm energy fr fr 👑

Purpose:

Using a divide and conquer approach, SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal N-by-M matrix B with diagonal D and offdiagonal E, where M = N + SQRE. The algorithm computes the singular values in the SVD B = U * S * VT. The orthogonal matrices U and VT are optionally computed in compact form.

A related subroutine, SLASD0, computes the singular values and the singular vectors in explicit form. cousins ✨

Parameters

ICOMPQ — INTEGER. Specifies whether singular vectors are to be computed in compact form, as follows:
- = 0: Compute singular values only.
- = 1: Compute singular vectors of upper bidiagonal matrix in compact form.
SMLSIZ — INTEGER. The maximum size of the subproblems at the bottom of the computation tree. the little leaves of the tree 🍓
N — INTEGER. The row dimension of the upper bidiagonal matrix. This is also the dimension of the main diagonal array D.
SQRE — INTEGER. Specifies the column dimension of the bidiagonal matrix.
- = 0: The bidiagonal matrix has column dimension M = N;
- = 1: The bidiagonal matrix has column dimension M = N + 1.
D — REAL array, dimension ( N ). On entry D contains the main diagonal of the bidiagonal matrix. On exit D, if INFO = 0, contains its singular values.
E — REAL array, dimension ( M-1 ). Contains the subdiagonal entries of the bidiagonal matrix. On exit, E has been destroyed. 💀 it served its purpose and then it was gone, no cap.
U — REAL array, dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left singular vector matrices of all subproblems at the bottom level.
LDU — INTEGER, LDU = > N. The leading dimension of arrays U, VT, DIFL, DIFR, POLES, GIVNUM, and Z. one dimension to rule them all 💅
VT — REAL array, dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT**T contains the right singular vector matrices of all subproblems at the bottom level. note the transpose, bestie, don't let it sneak past you 🦋
K — INTEGER array, dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th secular equation on the computation tree.
DIFL — REAL array, dimension ( LDU, NLVL ), where NLVL = floor(log_2 (N/SMLSIZ))).
DIFR
DIFR — a REAL array, dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and dimension ( N ) if ICOMPQ = 0. If ICOMPQ = 1, then on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1) record the distances between the singular values on the I-th level and the singular values on the (I -1)-th level, and DIFR(1:N, 2 * I ) contains the normalizing factors for the right singular vector matrix. See SLASD8 for details, bestie. 📖
Z — a REAL array, dimension ( LDU, NLVL ) if ICOMPQ = 1 and dimension ( N ) if ICOMPQ = 0. The first K elements of Z(1, I) contain the components of the deflation-adjusted updating row vector for subproblems on the I-th level. ✨
POLES — a REAL array, dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0 (she is simply not invited 💅). If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and POLES(1, 2*I) contain the new and old singular values involved in the secular equations on the I-th level.
GIVPTR — an INTEGER array, dimension ( N ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records the number of Givens rotations performed on the I-th problem on the computation tree. 🔄
GIVCOL — an INTEGER array, dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I, GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations of Givens rotations performed on the I-th level on the computation tree.
LDGCOL — an INTEGER, and literally: LDGCOL = > N. The leading dimension of arrays GIVCOL and PERM. (Yeah, that's how it's written in the source, no cap — we're just the messengers 🫶)
PERM — an INTEGER array, dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records permutations done on the I-th level of the computation tree.
GIVNUM — a REAL array, dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I, GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-values of Givens rotations performed on the I-th level on the computation tree.
C — a REAL array, dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1 and the I-th subproblem is not square, on exit, C( I ) contains the C-value of a Givens rotation related to the right null space of the I-th subproblem.
S — a REAL array, dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1 and the I-th subproblem is not square, on exit, S( I ) contains the S-value of a Givens rotation related to the right null space of the I-th subproblem.
WORK — a REAL array, dimension (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)). That's the workspace, and yes, those parentheses are doing important work. 🧮
IWORK — an INTEGER array, dimension (7*N).
INFO — an INTEGER, and it's the vibe check at the end:
- = 0: successful exit. 🎉 slay.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, a singular value did not converge. 😭

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

subroutine slasdq (character UPLO, integer SQRE, 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)

SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc.

Purpose:

SLASDQ computes the singular value decomposition (SVD) of a real (upper or lower) bidiagonal matrix with diagonal D and offdiagonal E, accumulating the transformations if desired. Letting B denote the input bidiagonal matrix, the algorithm computes orthogonal matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose of P). The singular values S are overwritten on D. 👑

The input matrix U is changed to U * Q if desired. The input matrix VT is changed to P**T * VT if desired. The input matrix C is changed to Q**T * C if desired.

Parameters

UPLO — UPLO is CHARACTER*1. On entry, UPLO specifies whether the input bidiagonal matrix is upper or lower bidiagonal, and whether it is square or not.
- UPLO = 'U' or 'u' — B is upper bidiagonal.
- UPLO = 'L' or 'l' — B is lower bidiagonal.
SQRE — SQRE is INTEGER.
- = 0: then the input matrix is N-by-N.
- = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and (N+1)-by-N if UPLU = 'L'. (spelled UPLU right there in the original — iconic typo behavior, we keep it 💀)

The bidiagonal matrix has N = NL + NR + 1 rows and M = N + SQRE >= N columns.

N — N is INTEGER. On entry, N specifies the number of rows and columns in the matrix. N must be at least 0.
NCVT — NCVT is INTEGER. On entry, NCVT specifies the number of columns of the matrix VT. NCVT must be at least 0.
NRU — NRU is INTEGER. On entry, NRU specifies the number of rows of the matrix U. NRU must be at least 0.
NCC — NCC is INTEGER. On entry, NCC specifies the number of columns of the matrix C. NCC must be at least 0.
D — D is a REAL array, dimension (N). On entry, D contains the diagonal entries of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in ascending order. 📈
E — E is a REAL array. Dimension is (N-1) if SQRE = 0 and N if SQRE = 1. On entry, the entries of E contain the offdiagonal entries of the bidiagonal matrix whose SVD is desired. On normal exit, E will contain 0. If the algorithm does not converge, D and E will contain the diagonal and superdiagonal entries of a bidiagonal matrix orthogonally equivalent to the one given as input. 😭
VT — VT is a REAL array, dimension (LDVT, NCVT). On entry, contains a matrix which on exit has been premultiplied by P**T, dimension N-by-NCVT if SQRE = 0 and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
LDVT — LDVT is INTEGER. On entry, LDVT specifies the leading dimension of VT as declared in the calling (sub) program. LDVT must be at least 1. If NCVT is nonzero LDVT must also be at least N.
U — U is a REAL array, dimension (LDU, N). On entry, contains a matrix which on exit has been postmultiplied by Q, dimension NRU-by-N if SQRE = 0 and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
LDU — LDU is INTEGER. On entry, LDU specifies the leading dimension of U as declared in the calling (sub) program. LDU must be at least max( 1, NRU ).
C — C is a REAL array, dimension (LDC, NCC). On entry, contains an N-by-NCC matrix which on exit has been premultiplied by Q**T, dimension N-by-NCC if SQRE = 0 and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
LDC — LDC is INTEGER. On entry, LDC specifies the leading dimension of C as declared in the calling (sub) program. LDC must be at least 1. If NCC is nonzero, LDC must also be at least N.
WORK — WORK is a REAL array, dimension (4*N). Workspace. Only referenced if one of NCVT, NRU, or NCC is nonzero, and if N is at least 2.
INFO — INFO is INTEGER. On exit, a value of 0 indicates a successful exit. If INFO < 0, argument number -INFO is illegal. If INFO > 0, the algorithm did not converge, and INFO specifies how many superdiagonals did not converge. 🫠

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 `slasdt` 🌈

subroutine slasdt (integer N, integer LVL, integer ND, integer, dimension( *
    ) INODE, integer, dimension( * ) NDIML, integer, dimension( * ) NDIMR,
    integer MSUB)

SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.

Purpose:

SLASDT creates a tree of subproblems for bidiagonal divide and conquer. 🌳

Parameters

N — N is INTEGER. On entry, the number of diagonal elements of the bidiagonal matrix.
LVL — LVL is INTEGER. On exit, the number of levels on the computation tree.
ND — ND is INTEGER. On exit, the number of nodes on the tree.
INODE — INODE is an INTEGER array, dimension ( N ). On exit, centers of subproblems. 🎯
NDIML — NDIML is an INTEGER array, dimension ( N ). On exit, row dimensions of left children.
NDIMR — NDIMR is an INTEGER array, dimension ( N ). On exit, row dimensions of right children.
MSUB
MSUB is INTEGER — on entry, this is the maximum row dimension that each subproblem at the bottom of the tree is allowed to be 🌳 it's giving "size cap for the lil baby subproblems at the leaves" and honestly? iconic.

🌸 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 slaset

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

Bestie, SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to values you pick. It's the "fill it in for me" queen 👑

✨ Purpose

SLASET initializes an m-by-n matrix A to BETA on the diagonal and
ALPHA on the offdiagonals.

🎀 Parameters

UPLO is CHARACTER*1 — tells you which chunk of matrix A is getting set 🦋
- = 'U': the upper triangular part is set; the strictly lower triangular part of A is not changed (hands off 💅)
- = 'L': the lower triangular part is set; the strictly upper triangular part of A is not changed
- Otherwise: all of the matrix A is set — every last crumb
M is INTEGER — the number of rows of the matrix A. M >= 0.
N is INTEGER — the number of columns of the matrix A. N >= 0.
ALPHA is REAL — the constant that the offdiagonal elements get set to.
BETA is REAL — the constant that the diagonal elements get set to.
A is REAL array, dimension (LDA,N) — on exit, the leading m-by-n submatrix of A is set as follows:

if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n,
if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n,
otherwise,     A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j,

and, for all UPLO:

A(i,i) = BETA, 1<=i<=min(m,n).

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

🌈 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🧚 subroutine slasr

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

SLASR applies a sequence of plane rotations to a general rectangular matrix 🔄 it's literally rotating, bestie.

✨ Purpose

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

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

A := P*A

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

A := A*P**T

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

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

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

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

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

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

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

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

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

where R(k) appears as a rank-2 modification to the identity matrix in rows and columns k and k+1.

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

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

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

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

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

where R(k) appears in rows and columns k and z. The rotations are performed without ever forming P(k) explicitly — efficient queen behavior, no cap 👑

🎀 Parameters

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

(M-1) if SIDE = 'L'
(N-1) if SIDE = 'R'

The cosines c(k) of the plane rotations.

S is REAL array, dimension

(M-1) if SIDE = 'L'
(N-1) if SIDE = 'R'

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

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

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

🌈 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

💫 subroutine slasv2

subroutine slasv2 (real F, real G, real H, real SSMIN, real SSMAX, real SNR, real CSR, real SNL, real CSL)

SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix 🔺 teeny tiny matrix, absolutely huge energy.

✨ Purpose

SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix

   [  F   G  ]
   [  0   H  ].

   [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
   [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].

🎀 Parameters

F is REAL — the (1,1) element of the 2-by-2 matrix.
G is REAL — the (1,2) element of the 2-by-2 matrix.
H is REAL — the (2,2) element of the 2-by-2 matrix.
SSMIN is REAL — abs(SSMIN) is the smaller singular value.
SSMAX is REAL — abs(SSMAX) is the larger singular value.
SNL is REAL
CSL is REAL — the vector (CSL, SNL) is a unit left singular vector for the singular value abs(SSMAX).
SNR is REAL
CSR is REAL — the vector (CSR, SNR) is a unit right singular vector for the singular value abs(SSMAX).

🌈 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

🔍 Further Details:

Any input parameter may be aliased with any output parameter. (chaotic? yes. allowed? also yes 💅)
Barring over/underflow and assuming a guard digit in subtraction, all output quantities are correct to within a few units in the last place (ulps).
In IEEE arithmetic, the code works correctly if one matrix element is infinite.
Overflow will not occur unless the largest singular value itself overflows or is within a few ulps of overflow. (On machines with partial overflow, like the Cray, overflow may occur if the largest singular value is within a factor of 2 of overflow.)
Underflow is harmless if underflow is gradual. Otherwise, results may correspond to a matrix modified by perturbations of size near the underflow threshold. ⚠️

🦋 subroutine xerbla

subroutine xerbla (character*(*) SRNAME, integer INFO)

XERBLA

✨ Purpose

XERBLA is an error handler for the LAPACK routines.

It is called by an LAPACK routine if an input parameter has an invalid value. A message is printed and execution stops. 😭

Installers may consider modifying the STOP statement in order to call system-specific exception-handling facilities.

🎀 Parameters

SRNAME

SRNAME is CHARACTER*(*) — the name of the routine which called XERBLA. 💅
INFO is INTEGER — the position of the invalid parameter in the parameter list of the calling routine. 🔥

Author ✨

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

subroutine xerbla_array (character(1), dimension(srname_len) SRNAME_ARRAY,
    integer SRNAME_LEN, integer INFO)

XERBLA_ARRAY 🌈

Purpose:

XERBLA_ARRAY assists other languages in calling XERBLA, the LAPACK and BLAS error handler. No cap, it's the ultimate translator bestie. Rather than taking a Fortran string argument as the function's name, XERBLA_ARRAY takes an array of single characters along with the array's length. XERBLA_ARRAY then copies up to 32 characters of that array into a Fortran string and passes that to XERBLA. If called with a non-positive SRNAME_LEN, XERBLA_ARRAY will call XERBLA with a string of all blank characters. 💀

Say some macro or other device makes XERBLA_ARRAY available to C99 by a name lapack_xerbla and with a common Fortran calling convention. Then a C99 program could invoke XERBLA via:

{
  int flen = strlen(__func__);
  lapack_xerbla(__func__, &flen, &info);
}

Providing XERBLA_ARRAY is not necessary for intercepting LAPACK errors. XERBLA_ARRAY calls XERBLA. 🫶

Parameters 🎀

SRNAME_ARRAY is CHARACTER(1) array, dimension (SRNAME_LEN) — the name of the routine which called XERBLA_ARRAY. 🦋
SRNAME_LEN is INTEGER — the length of the name in SRNAME_ARRAY. 💜
INFO is INTEGER — the position of the invalid parameter in the parameter list of the calling routine. 🌸

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.

And that's a wrap! 🎬 This page is OTHERauxiliary(3), version 3.10.0, dated Wed Jan 12 2022. Stay slaying, bestie! 💅✨

📦 where this page comes from

ubuntu jammy: liblapack-doc 3.10.0-2ubuntu1 · 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

OTHERauxiliary(3)

🌸 NAME 🌸

✨ SYNOPSIS ✨

🎀 Modules 🎀

✨ Functions ✨

✨ Detailed Description ✨

🌸 Function Documentation 🌸

💖 subroutine dcombssq (double precision, dimension( 2 ) V1, double precision, dimension( 2 ) V2)

Purpose:

Parameters

Author

🦋 logical function disnan (double precision, intent(in) DIN)

Purpose:

Parameters

Author

🌈 subroutine dlabad (double precision SMALL, double precision LARGE)

Purpose:

Parameters

Author

🌟 subroutine dlacpy (character UPLO, integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( ldb, * ) B, integer LDB)

Purpose:

Parameters

Author

💜 subroutine dlae2 (double precision A, double precision B, double precision C, double precision RT1, double precision RT2)

Purpose:

Parameters

Author 🫶

Further Details: 💖

subroutine dlaebz

Purpose: 🎯

Parameters

Author 👑

Further Details: 🌈

💅 subroutine dlaev2

🌈 subroutine dlagts

✨ logical function dlaisnan

🦋 integer function dlaneg

🌸 Author 🌸

🎀 Contributors 🎀

👑 dlanst — the norm queen

💜 dlapy2 — pythagoras but make it safe

🦋 dlapy3 — pythagoras but in 3D, main character energy

✨ dlarnv — the random number generator, no cap

🧚 dlarra — she finds the splitting points

💖 dlarrb — limited bisection, maximum accuracy

👑 Author

🫶 Contributors:

🦋 subroutine dlarrc

Purpose:

Parameters

Author

Contributors:

💜 subroutine dlarrd

Purpose:

Parameters

W 💖

WERR 🌸

WL 🦋

WU ✨

IBLOCK 💅

INDEXW 🎀

WORK 🧚

IWORK 🌈

INFO 💀

✨ Internal Parameters ✨

FUDGE 🍓

🌟 Contributors 🌟

👑 Author 👑

🦋 subroutine dlarre 🦋

🎯 Purpose 🎯

💅 Parameters 💅

RANGE

N

VL

VU

IL

IU

D

E

E2

💖 `subroutine dcombssq (double precision, dimension( 2 ) V1, double precision, dimension( 2 ) V2)`

🦋 `logical function disnan (double precision, intent(in) DIN)`

🌈 `subroutine dlabad (double precision SMALL, double precision LARGE)`

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

💜 `subroutine dlae2 (double precision A, double precision B, double precision C, double precision RT1, double precision RT2)`

💅 `subroutine dlaev2`

🌈 `subroutine dlagts`

✨ `logical function dlaisnan`

🦋 `integer function dlaneg`

👑 `dlanst` — the norm queen

💜 `dlapy2` — pythagoras but make it safe

🦋 `dlapy3` — pythagoras but in 3D, main character energy

✨ `dlarnv` — the random number generator, no cap

🧚 `dlarra` — she finds the splitting points

💖 `dlarrb` — limited bisection, maximum accuracy

`W` 💖

`WERR` 🌸

`WL` 🦋

`WU` ✨

`IBLOCK` 💅

`INDEXW` 🎀

`WORK` 🧚

`IWORK` 🌈

`INFO` 💀

`FUDGE` 🍓

🦋 `subroutine dlarre` 🦋

`RANGE`

`N`

`VL`

`VU`

`IL`

`IU`

`D`

`E`

`E2`

`RTOL1`

`RTOL2`

`SPLTOL`

`NSPLIT`

`ISPLIT`

`M`

`W`

`WERR`

`WGAP`

🧚 subroutine `dlarrf` 🧚

🦋 subroutine `dlarrj` 🦋

💫 `dlarrk`

🦋 `dlarrr`

🌸 `dlartgp`

💜 `dlaruv`

🧚 `dlas2`

👑 `dlascl`

🦋 `subroutine dlasd0`

🦋 `subroutine dlasd1`

🦋 `dlasd2` — the merge-and-deflate bestie 🦋

🌈 `dlasd3` — square roots of the secular equation, but make it fashion 🌈