realOTHEReigen(3)
💖 realOTHEReigen(3) 💖
okay bestie, welcome to LAPACK's realOTHEReigen(3) page — this is the "real Other Eigenvalue routines" corner of the library, and yes it's a LOT, so buckle up 🎀 (this is part 1 of 12, we're just getting started, don't you dare close the tab)
✨ NAME ✨
realOTHEReigen — real 💅
🦋 SYNOPSIS 🦋
🌸 Functions 🌸
here's the whole squad, fr fr. every signature, exactly as written, no cap:
SBDSVDX
SGGGLM
SSBEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSBGV
SSBGVD
SSBGVX
SSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSPGV
SSPGVD
SSPGVX
SSTEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSTEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
SSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
---
# 💜 Detailed Description 💜
This is the group of <pink>real Other Eigenvalue routines</pink> 🌸
periodt, that's the whole description, she's concise like that.
---
# 👑 Function Documentation 👑
## 🌈 subroutine sbdsvdx 🌈
SBDSVDX
⭐ Purpose ⭐
SBDSVDX computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT, where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices of left and right singular vectors, respectively.
Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ] and superdiagonal E = [ e_1 e_2 ... e_N-1 ], SBDSVDX computes the singular value decompositon of B through the eigenvalues and eigenvectors of the N*2-by-N*2 tridiagonal matrix
| 0 d_1 |
| d_1 0 e_1 |
TGK = | e_1 0 d_2 |
| d_2 . . |
| . . . |If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) / sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
Given a TGK matrix, one can either a) compute -s,-v and change signs so that the singular values (and corresponding vectors) are already in descending order (as in SGESVD/SGESDD) or b) compute s,v and reorder the values (and corresponding vectors). SBDSVDX implements a) by calling SSTEVX (bisection plus inverse iteration, to be replaced with a version of the Multiple Relative Robust Representation algorithm. (See P. Willems and B. Lang, A framework for the MR^3 algorithm: theory and implementation, SIAM J. Sci. Comput., 35:740-766, 2013.)
🎀 Parameters 🎀
UPLO
UPLO is CHARACTER*1
= 'U':B is upper bidiagonal;= 'L':B is lower bidiagonal.
JOBZ
JOBZ is CHARACTER*1
= 'N':Compute singular values only;= 'V':Compute singular values and singular vectors.
RANGE
RANGE is CHARACTER*1
= 'A':all singular values will be found.= 'V':all singular values in the half-open interval[VL,VU)will be found.= 'I':theIL-th throughIU-th singular values will be found.
N
N is INTEGER The order of the bidiagonal matrix. N >= 0.
D
D is REAL array, dimension (N) The n diagonal elements of the bidiagonal matrix B.
E
E is REAL array, dimension (max(1,N-1)) The (n-1) superdiagonal elements of the bidiagonal matrix B in elements 1 to N-1.
VL
VL is REAL If RANGE='V', the lower bound of the interval to be searched for singular values. VU > VL. Not referenced if RANGE = 'A' or 'I'.
VU
VU is REAL If RANGE='V', the upper bound of the interval to be searched for singular values. VU > VL. Not referenced if RANGE = 'A' or 'I'.
IL
IL is INTEGER If RANGE='I', the index of the smallest singular value to be returned. 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. Not referenced if RANGE = 'A' or 'V'.
IU
IU is INTEGER If RANGE='I', the index of the largest singular value to be returned. 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. Not referenced if RANGE = 'A' or 'V'.
NS
NS is INTEGER The total number of singular values found. 0 <= NS <= N. If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
S
S is REAL array, dimension (N) The first NS elements contain the selected singular values in ascending order.
Z
Z is REAL array, dimension (2*N,K) If JOBZ = 'V', then if INFO = 0 the first NS columns of Z contain the singular vectors of the matrix B corresponding to the selected singular values, with U in rows 1 to N and V in rows N+1 to N*2, i.e.
Z = [ U ]
[ V ]If JOBZ = 'N', then Z is not referenced.
Note: The user must ensure that at leastK = NS+1columns are supplied in the array Z; ifRANGE = 'V', the exact value of NS is not known in advance and an upper bound must be used. 😭 bestie this is the part where you double-check your array size, I'm begging
LDZ
LDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', LDZ >= max(2,N*2).
WORK
WORK is REAL array, dimension (14*N)
IWORK
IWORK is INTEGER array, dimension (12*N) If JOBZ = 'V', then if INFO = 0, the first NS elements of IWORK are zero. If INFO > 0, then IWORK contains the indices of the eigenvectors that failed to converge in DSTEVX.
INFO
INFO is INTEGER
= 0:successful exit< 0:ifINFO = -i, the i-th argument had an illegal value> 0:ifINFO = i, then i eigenvectors failed to converge in SSTEVX. The indices of the eigenvectors (as returned by SSTEVX) are stored in the array IWORK.- if
INFO = N*2 + 1, an internal error occurred.
💅 Author 💅
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
🦋 subroutine sggglm 🦋
SGGGLM
### ✨ Purpose ✨
SGGGLM solves a general Gauss-Markov linear model (GLM) problem:
minimize || y ||_2 subject to d = A*x + B*y x
where A is an N-by-M matrix, B is an N-by-P matrix, and d is a given N-vector. It is assumed that `M <= N <= M+P`, and
rank(A) = M and rank( A B ) = N.
🌸 (and that's where we're leaving it for this part, bestie — more coming, don't go anywhere 💖)
Under these assumptions, the constrained equation is **always consistent**, and there is a unique solution `x` and a minimal 2-norm solution `y` — and bestie, you get them via a generalized QR factorization of the matrices `(A, B)` given by:
A = Q*(R), B = Q*T*Z. (0)
In particular, if matrix `B` is square nonsingular, then the problem GLM is equivalent to the following weighted linear least squares problem
minimize || inv(B)*(d-A*x) ||_2 x
where `inv(B)` denotes the inverse of `B`. <pink>(yes, it's the inverse, we said what we said 💅)</pink>
### Parameters
- **`N`** — `N` is `INTEGER`. The number of rows of the matrices `A` and `B`. `N >= 0`. <cyan>no negative rows, that's illegal bestie</cyan>
- **`M`** — `M` is `INTEGER`. The number of columns of the matrix `A`. `0 <= M <= N`.
- **`P`** — `P` is `INTEGER`. The number of columns of the matrix `B`. `P >= N-M`.
- **`A`** — `A` is `REAL` array, dimension `(LDA,M)`.
- On entry, the `N`-by-`M` matrix `A`.
- On exit, the upper triangular part of the array `A` contains the `M`-by-`M` upper triangular matrix `R`. <purple>she's been transformed, she's a different girl now ✨</purple>
- **`LDA`** — `LDA` is `INTEGER`. The leading dimension of the array `A`. `LDA >= max(1,N)`.
- **`B`** — `B` is `REAL` array, dimension `(LDB,P)`.
- On entry, the `N`-by-`P` matrix `B`.
- On exit, if `N <= P`, the upper triangle of the subarray `B(1:N,P-N+1:P)` contains the `N`-by-`N` upper triangular matrix `T`; if `N > P`, the elements on and above the `(N-P)`th subdiagonal contain the `N`-by-`P` upper trapezoidal matrix `T`.
- **`LDB`** — `LDB` is `INTEGER`. The leading dimension of the array `B`. `LDB >= max(1,N)`.
- **`D`** — `D` is `REAL` array, dimension `(N)`.
- On entry, `D` is the left hand side of the GLM equation.
- On exit, `D` is **destroyed**. <red>gone. reduced to atoms. do not go looking for her afterwards 😭</red>
- **`X`** — `X` is `REAL` array, dimension `(M)`.
- **`Y`** — `Y` is `REAL` array, dimension `(P)`.
- On exit, `X` and `Y` are the solutions of the GLM problem. <rainbow>the girls ATE 🫶</rainbow>
- **`WORK`** — `WORK` is `REAL` array, dimension `(MAX(1,LWORK))`.
- On exit, if `INFO = 0`, `WORK(1)` returns the optimal `LWORK`.
- **`LWORK`** — `LWORK` is `INTEGER`. The dimension of the array `WORK`. `LWORK >= max(1,N+M+P)`.
- For optimum performance, `LWORK >= M+min(N,P)+max(N,P)*NB`, where `NB` is an upper bound for the optimal blocksizes for `SGEQRF`, `SGERQF`, `SORMQR` and `SORMRQ`.
- If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal size of the `WORK` array, returns this value as the first entry of the `WORK` array, and no error message related to `LWORK` is issued by `XERBLA`. <mint>a little sneak peek, no commitment, we love a soft launch</mint>
- **`INFO`** — `INFO` is `INTEGER`.
- `= 0`: successful exit. <green>slay</green>
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value.
- `= 1`: the upper triangular factor `R` associated with `A` in the generalized QR factorization of the pair `(A, B)` is singular, so that `rank(A) < M`; the least squares solution could not be computed.
- `= 2`: the bottom `(N-M)` by `(N-M)` part of the upper trapezoidal factor `T` associated with `B` in the generalized QR factorization of the pair `(A, B)` is singular, so that `rank( A B ) < N`; the least squares solution could not be computed.
### Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
## subroutine `ssbev`
subroutine ssbev (character JOBZ, character UPLO, integer N, integer KD, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)
<hotpink>**SSBEV**</hotpink> computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices. <purple>OTHER matrices?? yes bestie, OTHER. that's the vibe, don't question it 🦋</purple>
### Purpose
`SSBEV` computes **all** the eigenvalues and, optionally, eigenvectors of a real symmetric band matrix `A`. <sparkle>every single one, we're thorough here</sparkle>
### Parameters
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: Upper triangle of `A` is stored;
- `= 'L'`: Lower triangle of `A` is stored.
- **`N`** — `N` is `INTEGER`. The order of the matrix `A`. `N >= 0`.
- **`KD`** — `KD` is `INTEGER`. The number of superdiagonals of the matrix `A` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KD >= 0`.
- **`AB`** — `AB` is `REAL` array, dimension `(LDAB, N)`.
- On entry, the upper or lower triangle of the symmetric band matrix `A`, stored in the first `KD+1` rows of the array. The j-th column of `A` is stored in the j-th column of the array `AB` as follows:
- if `UPLO = 'U'`, `AB(kd+1+i-j,j) = A(i,j)` for `max(1,j-kd)<=i<=j`;
- if `UPLO = 'L'`, `AB(1+i-j,j) = A(i,j)` for `j<=i<=min(n,j+kd).`
- On exit, `AB` is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the first superdiagonal and the diagonal of the tridiagonal matrix `T` are returned in rows `KD` and `KD+1` of `AB`, and if `UPLO = 'L'`, the diagonal and first subdiagonal of `T` are returned in the first two rows of `AB`. <pink>she got a glow-up mid-routine, iconic 💖</pink>
- **`LDAB`** — `LDAB` is `INTEGER`. The leading dimension of the array `AB`. `LDAB >= KD + 1`.
- **`W`** — `W` is `REAL` array, dimension `(N)`. If `INFO = 0`, the eigenvalues in ascending order. <cyan>sorted, organized, main character energy 📋</cyan>
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, N)`.
- If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the orthonormal eigenvectors of the matrix `A`, with the i-th column of `Z` holding the eigenvector associated with `W(i)`.
- If `JOBZ = 'N'`, then `Z` is not referenced. <lavender>she's not even invited to the function, no hard feelings</lavender>
- **`LDZ`** — `LDZ` is `INTEGER`. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL` array, dimension `(max(1,3*N-2))`
- **`INFO`** — `INFO` is `INTEGER`
- `= 0`: successful exit
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value
- `> 0`: if `INFO = i`, the algorithm failed to converge; `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero. <red>it just didn't converge bestie, we tried, we really did 💀</red>
### Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
## subroutine `ssbev_2stage`
subroutine ssbev_2stage (character JOBZ, character UPLO, integer N, integer KD, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer INFO)
<hotpink>**SSBEV_2STAGE**</hotpink> computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices. <purple>the sequel. the two-stage remix. she came back with a whole new technique 🎀</purple>
### Purpose
`SSBEV_2STAGE` computes all the eigenvalues and, optionally, eigenvectors of a real symmetric band matrix `A` using the **2stage technique** for the reduction to tridiagonal. <sparkle>two stages, twice the drama, same slay</sparkle>
### Parameters
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- <yellow>**Not available in this release.**</yellow> <red>ma'am we do not have that yet 🙈</red>
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: Upper triangle of `A` is stored;
- `= 'L'`: Lower triangle of `A` is stored.
- **`N`** — `N` is `INTEGER`. The order of the matrix `A`. `N >= 0`.
- **`KD`** — `KD` is `INTEGER`. The number of superdiagonals of the matrix `A` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KD >= 0`.
- **`AB`** — `AB` is `REAL` array, dimension `(LDAB, N)`.
- On entry, the upper or lower triangle of the symmetric band matrix `A`, stored in the first `KD+1` rows of the array. The j-th column of `A` is stored in the j-th column of the array `AB` as follows:
- if `UPLO = 'U'`, `AB(kd+1+i-j,j) = A(i,j)` for `max(1,j-kd)<=i<=j`;
- if `UPLO = 'L'`, `AB(1+i-j,j) = A(i,j)` for `j<=i<=min(n,j+kd).`
- On exit, `AB` is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the first superdiagonal and the diagonal of the tridiagonal matrix `T` are returned in rows `KD` and `KD+1` of `AB`, and if `UPLO = 'L'`, the diagonal and first subdiagonal of `T` are returned in the first two rows of `AB`. <pink>same glow-up, different technique, consistency is key ✨</pink>
- **`LDAB`** — `LDAB` is `INTEGER`. The leading dimension of the array `AB`. `LDAB >= KD + 1`.
- **`W`** — `W` is `REAL` array, dimension `(N)`. If `INFO = 0`, the eigenvalues in ascending order.
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, N)`.
- If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the orthonormal eigenvectors of the matrix `A`, with the i-th column of `Z` holding the eigenvector associated with `W(i)`.
- If `JOBZ = 'N'`, then `Z` is not referenced.
- **`LDZ`** — `LDZ` is `INTEGER`. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL` array, dimension `LWORK`.
- On exit, if `INFO = 0`, `WORK(1)` returns the optimal `LWORK`.
- **`LWORK`** — `LWORK` is `INTEGER`. The length of the array `WORK`. `LWORK >= 1`, when `N <= 1`; otherwise
- If `JOBZ = 'N'` and `N > 1`, `LWORK` must be queried.
- `LWORK = MAX(1, dimension)` where
- `dimension = (2KD+1)*N + KD*NTHREADS + N`
- where `KD` is the size of the band.
- `NTHREADS` is the number of threads used when openMP compilation is enabled, otherwise `=1`.
- If `JOBZ = 'V'` and `N > 1`, `LWORK` must be queried. **Not yet available.** <lavender>coming soon to a release near you, we promise nothing yet bestie 🦋</lavender>
If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal size of the `WORK` array, returns this value as the first entry of the `WORK` array, and no error message related to `LWORK` is issued by `XERBLA`. <rainbow>we love a lil sneak peek at the workspace 💅</rainbow>
- **`INFO`** is `INTEGER` — the status report bestie, read it carefully:
- `= 0`: successful exit ✨
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value 🙈
- `> 0`: if `INFO = i`, the algorithm failed to converge; `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero 💀
**Author**
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
**Further Details:**
All details about the 2stage techniques are available in:
- Azzam Haidar, Hatem Ltaief, and Jack Dongarra. *Parallel reduction to condensed forms for symmetric eigenvalue problems using aggregated fine-grained and memory-aware kernels.* In Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis (SC '11), New York, NY, USA, Article 8 , 11 pages. `http://doi.acm.org/10.1145/2063384.2063394`
- A. Haidar, J. Kurzak, P. Luszczek, 2013. *An improved parallel singular value algorithm and its implementation for multicore hardware*, In Proceedings of 2013 International Conference for High Performance Computing, Networking, Storage and Analysis (SC '13). Denver, Colorado, USA, 2013. Article 90, 12 pages. `http://doi.acm.org/10.1145/2503210.2503292`
- A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra. *A novel hybrid CPU-GPU generalized eigensolver for electronic structure calculations based on fine-grained memory aware tasks.* International Journal of High Performance Computing Applications. Volume 28 Issue 2, Pages 196-209, May 2014. `http://hpc.sagepub.com/content/28/2/196`
---
### ✨ subroutine `ssbevd` ✨
subroutine ssbevd (character JOBZ, character UPLO, integer N, integer KD, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
<cyan>SSBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</cyan> 🦋
**Purpose:**
Okay so `SSBEVD` computes <pink>all</pink> the eigenvalues and, if you want them, the eigenvectors of a real symmetric band matrix `A`. If eigenvectors are desired, it uses a divide and conquer algorithm. No cap. 💅
The divide and conquer algorithm makes <purple>very mild assumptions</purple> about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. <yellow>not one single one bestie</yellow> 😭
**Parameters**
- **`JOBZ`** is `CHARACTER*1` — the whole personality of this call:
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`UPLO`** is `CHARACTER*1` — which triangle you're giving me:
- `= 'U'`: Upper triangle of A is stored;
- `= 'L'`: Lower triangle of A is stored.
- **`N`** is `INTEGER`
- The order of the matrix A. `N >= 0`.
- **`KD`** is `INTEGER`
- The number of superdiagonals of the matrix A if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KD >= 0`.
- **`AB`** is `REAL` array, dimension `(LDAB, N)`
- On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first `KD+1` rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:
- if `UPLO = 'U'`, `AB(kd+1+i-j,j) = A(i,j)` for `max(1,j-kd)<=i<=j`;
- if `UPLO = 'L'`, `AB(1+i-j,j) = A(i,j)` for `j<=i<=min(n,j+kd)`.
- On exit, AB is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the first superdiagonal and the diagonal of the tridiagonal matrix T are returned in rows `KD` and `KD+1` of AB, and if `UPLO = 'L'`, the diagonal and first subdiagonal of T are returned in the first two rows of AB. <hotpink>she does not come back the way you left her 💀</hotpink>
- **`LDAB`** is `INTEGER`
- The leading dimension of the array AB. `LDAB >= KD + 1`.
- **`W`** is `REAL` array, dimension `(N)`
- If `INFO = 0`, the eigenvalues in ascending order. 🍓
- **`Z`** is `REAL` array, dimension `(LDZ, N)`
- If `JOBZ = 'V'`, then if `INFO = 0`, Z contains the orthonormal eigenvectors of the matrix A, with the i-th column of Z holding the eigenvector associated with `W(i)`.
- If `JOBZ = 'N'`, then Z is not referenced. <lavender>she's just there, vibing, untouched</lavender>
- **`LDZ`** is `INTEGER`
- The leading dimension of the array Z. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** is `REAL` array, dimension `(LWORK)`
- On exit, if `INFO = 0`, `WORK(1)` returns the optimal `LWORK`.
- **`LWORK`** is `INTEGER`
- The dimension of the array `WORK`.
- If `N <= 1`, `LWORK` must be at least `1`.
- If `JOBZ = 'N'` and `N > 2`, `LWORK` must be at least `2*N`.
- If `JOBZ = 'V'` and `N > 2`, `LWORK` must be at least `( 1 + 5*N + 2*N**2 )`.
- If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`.
- **`IWORK`** is `INTEGER` array, dimension `(MAX(1,LIWORK))`
- On exit, if `INFO = 0`, `IWORK(1)` returns the optimal `LIWORK`.
- **`LIWORK`** is `INTEGER`
- The dimension of the array `IWORK`.
- If `JOBZ = 'N'` or `N <= 1`, `LIWORK` must be at least `1`.
- If `JOBZ = 'V'` and `N > 2`, `LIWORK` must be at least `3 + 5*N`.
- If `LIWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`.
- **`INFO`** is `INTEGER`
- `= 0`: successful exit ✨
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value 🙈
- `> 0`: if `INFO = i`, the algorithm failed to converge; `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero 💀
**Author**
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
### ✨ subroutine `ssbevd_2stage` ✨
subroutine ssbevd_2stage (character JOBZ, character UPLO, integer N, integer KD, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
<cyan>SSBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</cyan> 🧚
**Purpose:**
`SSBEVD_2STAGE` computes all the eigenvalues and, optionally, eigenvectors of a real symmetric band matrix `A` <pink>using the 2stage technique</pink> for the reduction to tridiagonal. If eigenvectors are desired, it uses a divide and conquer algorithm. <rainbow>two stages?? we stan a multi-step queen 👑</rainbow>
The divide and conquer algorithm makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none.
**Parameters**
- **`JOBZ`** is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- <red>Not available in this release.</red> 😭 she said not yet bestie
- **`UPLO`** is `CHARACTER*1`
- `= 'U'`: Upper triangle of A is stored;
- `= 'L'`: Lower triangle of A is stored.
- **`N`** is `INTEGER`
- The order of the matrix A. `N >= 0`.
- **`KD`** is `INTEGER`
- The number of superdiagonals of the matrix A if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KD >= 0`.
- **`AB`** is `REAL` array, dimension `(LDAB, N)`
- On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first `KD+1` rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:
- if `UPLO = 'U'`, `AB(kd+1+i-j,j) = A(i,j)` for `max(1,j-kd)<=i<=j`;
- if `UPLO = 'L'`, `AB(1+i-j,j) = A(i,j)` for `j<=i<=min(n,j+kd)`.
- On exit, AB is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the first superdiagonal and the diagonal of the tridiagonal matrix T are returned in rows `KD` and `KD+1` of AB, and if `UPLO = 'L'`, the diagonal and first subdiagonal of T are returned in the first two rows of AB.
- **`LDAB`** is `INTEGER`
- The leading dimension of the array AB. `LDAB >= KD + 1`.
- **`W`** is `REAL` array, dimension `(N)`
- If `INFO = 0`, the eigenvalues in ascending order. 🍓
- **`Z`** is `REAL` array, dimension `(LDZ, N)`
- If `JOBZ = 'V'`, then if `INFO = 0`, Z contains the orthonormal eigenvectors of the matrix A, with the i-th column of Z holding the eigenvector associated with `W(i)`.
- If `JOBZ = 'N'`, then Z is not referenced.
- **`LDZ`** is `INTEGER`
- The leading dimension of the array Z. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** —
- **`WORK`** is `REAL` array, dimension `LWORK`. On exit, if `INFO = 0`, `WORK(1)` hands you back the optimal `LWORK`. 💖
- **`LWORK`** — `LWORK` is `INTEGER`. It's literally just how long the `WORK` array is. Bestie, the rules here are kinda ✨conditional✨:
- `LWORK >= 1`, when `N <= 1`; otherwise…
- If `JOBZ = 'N'` and `N > 1`, `LWORK` **must be queried**. 🚨
- `LWORK = MAX(1, dimension)` where
- `dimension = (2KD+1)*N + KD*NTHREADS + N`
- where `KD` is the size of the band.
- `NTHREADS` is the number of threads used when openMP compilation is enabled, otherwise `=1`.
- If `JOBZ = 'V'` and `N > 1`, `LWORK` must be queried. **Not yet available.** 💀 (fr fr, it's just not there yet)
- If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`. So `-1` is the "bestie just tell me how much space I need" move. 🫶
- **`IWORK`** — `IWORK` is `INTEGER` array, dimension `(MAX(1,LIWORK))`. On exit, if `INFO = 0`, `IWORK(1)` returns the optimal `LIWORK`. 👑
- **`LIWORK`** — `LIWORK` is `INTEGER`. The dimension of the array `IWORK`.
- If `JOBZ = 'N'` or `N <= 1`, `LIWORK` must be at least `1`.
- If `JOBZ = 'V'` and `N > 2`, `LIWORK` must be at least `3 + 5*N`.
- If `LIWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`. Same energy as above, just for the int array. ✨
- **`INFO`** — `INFO` is `INTEGER`:
- `= 0`: successful exit. slay. 🎀
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value. not me passing garbage 😭
- `> 0`: if `INFO = i`, the algorithm failed to converge; `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero. it tried, it really did. 💀
---
## 👑 Author 👑
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
The whole academic friend group, no cap. 🫶
---
## 🔬 Further Details: 🔬
All the tea about the 2stage techniques lives in these, verbatim:
Azzam Haidar, Hatem Ltaief, and Jack Dongarra. Parallel reduction to condensed forms for symmetric eigenvalue problems using aggregated fine-grained and memory-aware kernels. In Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis (SC '11), New York, NY, USA, Article 8 , 11 pages. http://doi.acm.org/10.1145/2063384.2063394
A. Haidar, J. Kurzak, P. Luszczek, 2013. An improved parallel singular value algorithm and its implementation for multicore hardware, In Proceedings of 2013 International Conference for High Performance Computing, Networking, Storage and Analysis (SC '13). Denver, Colorado, USA, 2013. Article 90, 12 pages. http://doi.acm.org/10.1145/2503210.2503292
A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra. A novel hybrid CPU-GPU generalized eigensolver for electronic structure calculations based on fine-grained memory aware tasks. International Journal of High Performance Computing Applications. Volume 28 Issue 2, Pages 196-209, May 2014. http://hpc.sagepub.com/content/28/2/196
---
# 🌈 subroutine `ssbevx` 🌈
subroutine ssbevx (character JOBZ, character RANGE, character UPLO, integer N, integer KD, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldq, * ) Q, integer LDQ, real VL, real VU, integer IL, integer IU, real ABSTOL, integer M, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer, dimension( * ) IFAIL, integer INFO)
`SSBEVX` computes the eigenvalues and, optionally, the left and/or right eigenvectors for **OTHER** matrices. yes, that's the whole vibe of this one. 💅
## 🎯 Purpose: 🎯
`SSBEVX` computes selected eigenvalues and, optionally, eigenvectors of a **real symmetric band matrix** `A`. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues. So you get to *choose*, which is very main character of you. ✨
## 🌸 Parameters 🌸
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors. (the full package 💖)
- **`RANGE`** — `RANGE` is `CHARACTER*1`
- `= 'A'`: all eigenvalues will be found;
- `= 'V'`: all eigenvalues in the half-open interval `(VL,VU]` will be found;
- `= 'I'`: the `IL`-th through `IU`-th eigenvalues will be found.
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: Upper triangle of `A` is stored;
- `= 'L'`: Lower triangle of `A` is stored.
- **`N`** — `N` is `INTEGER`. The order of the matrix `A`. `N >= 0`.
- **`KD`** — `KD` is `INTEGER`. The number of superdiagonals of the matrix `A` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KD >= 0`.
- **`AB`** — `AB` is `REAL` array, dimension `(LDAB, N)`.
On entry, the upper or lower triangle of the symmetric band matrix `A`, stored in the first `KD+1` rows of the array. The `j`-th column of `A` is stored in the `j`-th column of the array `AB` as follows:
if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
On exit, `AB` is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the first superdiagonal and the diagonal of the tridiagonal matrix `T` are returned in rows `KD` and `KD+1` of `AB`, and if `UPLO = 'L'`, the diagonal and first subdiagonal of `T` are returned in the first two rows of `AB`. your input? gone. replaced. glow up. 💅
- **`LDAB`** — `LDAB` is `INTEGER`. The leading dimension of the array `AB`. `LDAB >= KD + 1`.
- **`Q`** — `Q` is `REAL` array, dimension `(LDQ, N)`.
- If `JOBZ = 'V'`, the `N`-by-`N` orthogonal matrix used in the reduction to tridiagonal form.
- If `JOBZ = 'N'`, the array `Q` is not referenced. ignored, as it should be. 🙈
- **`LDQ`** — `LDQ` is `INTEGER`. The leading dimension of the array `Q`. If `JOBZ = 'V'`, then `LDQ >= max(1,N)`.
- **`VL`** — `VL` is `REAL`. If `RANGE='V'`, the lower bound of the interval to be searched for eigenvalues. `VL < VU`. Not referenced if `RANGE = 'A'` or `'I'`.
- **`VU`** — `VU` is `REAL`. If `RANGE='V'`, the upper bound of the interval to be searched for eigenvalues. `VL < VU`. Not referenced if `RANGE = 'A'` or `'I'`.
- **`IL`** — `IL` is `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`** — `IU` is `INTEGER`. If `RANGE='I'`, the index of the largest eigenvalue to be returned. `1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`. Not referenced if `RANGE = 'A'` or `'V'`.
- **`ABSTOL`** — `ABSTOL` is `REAL`. The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval `[a,b]` of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where `EPS` is the machine precision. If `ABSTOL` is less than or equal to zero, then `EPS*|T|` will be used in its place, where `|T|` is the 1-norm of the tridiagonal matrix obtained by reducing `AB` to tridiagonal form.
Eigenvalues will be computed most accurately when `ABSTOL` is set to twice the underflow threshold `2*SLAMCH('S')`, **not zero**. If this routine returns with `INFO>0`, indicating that some eigenvectors did not converge, try setting `ABSTOL` to `2*SLAMCH('S')`. 👀
See *'Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy,'* by Demmel and Kahan, LAPACK Working Note #3.
- **`M`** — `M` is `INTEGER`. The total number of eigenvalues found. `0 <= M <= N`. If `RANGE = 'A'`, `M = N`, and if `RANGE = 'I'`, `M = IU-IL+1`.
- **`W`** — `W` is `REAL` array, dimension `(N)`. The first `M` elements contain the selected eigenvalues in ascending order. organized. sorted. iconic. ✨
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, max(1,M))`.
- If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` columns of `Z` contain the orthonormal eigenvectors of the matrix `A` corresponding to the selected eigenvalues, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`.
- If an eigenvector fails to converge, then that column of `Z` contains the latest approximation to the eigenvector, and the index of the eigenvector is returned in `IFAIL`. (close enough, bestie 😭)
- If `JOBZ = 'N'`, then `Z` is not referenced.
- **Note:** the user must ensure that at least `max(1,M)` columns are supplied in the array `Z`; if `RANGE = 'V'`, the exact value of `M` is not known in advance and an upper bound must be used.
- **`LDZ`** — `LDZ` is `INTEGER`. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL` array, dimension `(7*N)`.
- **`IWORK`** — `IWORK` is `INTEGER` array, dimension `(5*N)`.
- **`IFAIL`** — `IFAIL` is `INTEGER` array, dimension `(N)`.
- If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` elements of `IFAIL` are zero.
- If `INFO > 0`, then `IFAIL` contains the indices of the eigenvectors that failed to converge. the receipts, basically 🧾
- If `JOBZ = 'N'`, then `IFAIL` is not referenced.
- **`INFO`** — …
- **`INFO`** is `INTEGER`
- `= 0`: <green>successful exit</green> ✨ we love a clean run
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value 🙈 (so like, argument number *i* was a mess)
- `> 0`: if `INFO = i`, then i eigenvectors failed to converge 😭 and their indices are stored in array `IFAIL` — go check there bestie
---
### Author 🌸
<purple>Univ. of Tennessee</purple>
<purple>Univ. of California Berkeley</purple>
<purple>Univ. of Colorado Denver</purple>
<purple>NAG Ltd.</purple>
---
### Further Details: 📚
All the tea about the 2stage techniques lives in these papers, okay? Like, the full lore:
> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
> Parallel reduction to condensed forms for symmetric eigenvalue problems
> using aggregated fine-grained and memory-aware kernels. In Proceedings
> of 2011 International Conference for High Performance Computing,
> Networking, Storage and Analysis (SC '11), New York, NY, USA,
> Article 8 , 11 pages.
> http://doi.acm.org/10.1145/2063384.2063394
> A. Haidar, J. Kurzak, P. Luszczek, 2013.
> An improved parallel singular value algorithm and its implementation
> for multicore hardware, In Proceedings of 2013 International Conference
> for High Performance Computing, Networking, Storage and Analysis (SC '13).
> Denver, Colorado, USA, 2013.
> Article 90, 12 pages.
> http://doi.acm.org/10.1145/2503210.2503292
> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
> calculations based on fine-grained memory aware tasks.
> International Journal of High Performance Computing Applications.
> Volume 28 Issue 2, Pages 196-209, May 2014.
> http://hpc.sagepub.com/content/28/2/196
---
### subroutine ssbgv 🦋
Okay new function dropped, and she's giving *generalized* symmetric band problem energy:
subroutine ssbgv (character JOBZ, character UPLO, integer N, integer KA, integer KB, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldbb, * ) BB, integer LDBB, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)
<pink>**SSBGV**</pink>
#### Purpose: 💖
So here's the tea 🫖: `SSBGV` computes all the eigenvalues, and optionally the eigenvectors, of a real generalized symmetric-definite banded eigenproblem, of the form `A*x=(lambda)*B*x`. Both `A` and `B` are assumed to be symmetric and banded, and `B` is also **positive definite**. That last bit is not optional bestie, it's load-bearing 💅.
### ✨ Parameters ✨
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: compute eigenvalues only;
- `= 'V'`: compute eigenvalues *and* eigenvectors.
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: upper triangles of `A` and `B` are stored;
- `= 'L'`: lower triangles of `A` and `B` are stored.
- **`N`** — `N` is `INTEGER`. The order of the matrices `A` and `B`. `N >= 0`.
- **`KA`** — `KA` is `INTEGER`. The number of superdiagonals of the matrix `A` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KA >= 0`.
- **`KB`** — `KB` is `INTEGER`. The number of superdiagonals of the matrix `B` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KB >= 0`.
- **`AB`** — `AB` is `REAL` array, dimension `(LDAB, N)`. On entry, the upper or lower triangle of the symmetric band matrix `A`, stored in the first `ka+1` rows of the array. The `j`-th column of `A` is stored in the `j`-th column of the array `AB` as follows:
- if `UPLO = 'U'`, `AB(ka+1+i-j,j) = A(i,j)` for `max(1,j-ka)<=i<=j`;
- if `UPLO = 'L'`, `AB(1+i-j,j) = A(i,j)` for `j<=i<=min(n,j+ka)`.
On exit, the contents of `AB` are <red>**destroyed**</red> 💀 like actually gone, do not go looking for your original matrix here.
- **`LDAB`** — `LDAB` is `INTEGER`. The leading dimension of the array `AB`. `LDAB >= KA+1`.
- **`BB`** — `BB` is `REAL` array, dimension `(LDBB, N)`. On entry, the upper or lower triangle of the symmetric band matrix `B`, stored in the first `kb+1` rows of the array. The `j`-th column of `B` is stored in the `j`-th column of the array `BB` as follows:
- if `UPLO = 'U'`, `BB(kb+1+i-j,j) = B(i,j)` for `max(1,j-kb)<=i<=j`;
- if `UPLO = 'L'`, `BB(1+i-j,j) = B(i,j)` for `j<=i<=min(n,j+kb)`.
On exit, `BB` holds the factor `S` from the split Cholesky factorization `B = S**T*S`, as returned by `SPBSTF`. Glow up ✨.
- **`LDBB`** — `LDBB` is `INTEGER`. The leading dimension of the array `BB`. `LDBB >= KB+1`.
- **`W`** — `W` is `REAL` array, dimension `(N)`. If `INFO = 0`, the eigenvalues in ascending order. Sorted and iconic 👑.
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, N)`. If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the matrix `Z` of eigenvectors, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`. The eigenvectors are normalized so that `Z**T*B*Z = I`. If `JOBZ = 'N'`, then `Z` is **not referenced** — it's giving "you don't know her" 🙈.
- **`LDZ`** — `LDZ` is `INTEGER`. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= N`.
- **`WORK`** — `WORK` is `REAL` array, dimension `(3*N)`.
- **`INFO`** — `INFO` is `INTEGER`
- `= 0`: successful exit 🎉
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value
- `> 0`: if `INFO = i`, and `i` is:
- `<= N`: the algorithm failed to converge: `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero;
- `> N`: if `INFO = N + i`, for `1 <= i <= N`, then `SPBSTF` returned `INFO = i`: `B` is not positive definite. The factorization of `B` could not be completed and no eigenvalues or eigenvectors were computed. <red>Rough.</red> 😭
### ✨ Author ✨
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
subroutine ssbgvd (character JOBZ, character UPLO, integer N, integer KA, integer KB, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldbb, * ) BB, integer LDBB, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO) SSBGVD
### ✨ Purpose ✨
`SSBGVD` computes all the eigenvalues, and optionally, the eigenvectors of a real generalized symmetric-definite banded eigenproblem, of the form `A*x=(lambda)*B*x`. Here `A` and `B` are assumed to be symmetric and banded, and `B` is also positive definite. If eigenvectors are desired, it uses a **divide and conquer algorithm**. Yes bestie, it ate and left no crumbs 🔥.
The divide and conquer algorithm makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the `Cray X-MP`, `Cray Y-MP`, `Cray C-90`, or `Cray-2`. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. So like, vintage Cray energy only, it's fine ✨.
### ✨ Parameters ✨
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: compute eigenvalues only;
- `= 'V'`: compute eigenvalues *and* eigenvectors.
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: upper triangles of `A` and `B` are stored;
- `= 'L'`: lower triangles of `A` and `B` are stored.
- **`N`** — `N` is `INTEGER`. The order of the matrices `A` and `B`. `N >= 0`.
- **`KA`** — `KA` is `INTEGER`. The number of superdiagonals of the matrix `A` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KA >= 0`.
- **`KB`** — `KB` is `INTEGER`. The number of superdiagonals of the matrix `B` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. `KB >= 0`.
- **`AB`** — `AB` is `REAL` array, dimension `(LDAB, N)`. On entry, the upper or lower triangle of the symmetric band matrix `A`, stored in the first `ka+1` rows of the array. The `j`-th column of `A` is stored in the `j`-th column of the array `AB` as follows:
- if `UPLO = 'U'`, `AB(ka+1+i-j,j) = A(i,j)` for `max(1,j-ka)<=i<=j`;
- if `UPLO = 'L'`, `AB(1+i-j,j) = A(i,j)` for `j<=i<=min(n,j+ka)`.
On exit, the contents of `AB` are <red>**destroyed**</red> 💀.
- **`LDAB`** — `LDAB` is `INTEGER`. The leading dimension of the array `AB`. `LDAB >= KA+1`.
- **`BB`** — `BB` is `REAL` array, dimension `(LDBB, N)`. On entry, the upper or lower triangle of the symmetric band matrix `B`, stored in the first `kb+1` rows of the array. The `j`-th column of `B` is stored in the `j`-th column of the array `BB` as follows:
- if `UPLO = 'U'`, `BB(ka+1+i-j,j) = B(i,j)` for `max(1,j-kb)<=i<=j`;
- if `UPLO = 'L'`, `BB(1+i-j,j) = B(i,j)` for `j<=i<=min(n,j+kb)`.
On exit, `BB` holds the factor `S` from the split Cholesky factorization `B = S**T*S`, as returned by `SPBSTF`.
- **`LDBB`** — `LDBB` is `INTEGER`. The leading dimension of the array `BB`. `LDBB >= KB+1`.
- **`W`** — `W` is `REAL` array, dimension `(N)`. If `INFO = 0`, the eigenvalues in ascending order.
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, N)`. If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the matrix `Z` of eigenvectors, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`. The eigenvectors are normalized so `Z**T*B*Z = I`. If `JOBZ = 'N'`, then `Z` is not referenced.
- **`LDZ`** — `LDZ` is `INTEGER`. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL` array, dimension `(MAX(1,LWORK))`. On exit, if `INFO = 0`, `WORK(1)` returns the optimal `LWORK`.
- **`LWORK`** — `LWORK` is `INTEGER`. The dimension of the array `WORK`.
- If `N <= 1`, `LWORK >= 1`.
- If `JOBZ = 'N'` and `N > 1`, `LWORK >= 3*N`.
- If `JOBZ = 'V'` and `N > 1`, `LWORK >= 1 + 5*N + 2*N**2`.
If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`. It's giving "just asking, not committing" 🫶.
- **`IWORK`** — `IWORK` is `INTEGER` array, dimension `(MAX(1,LIWORK))`. On exit, if `LIWORK > 0`, `IWORK(1)` returns the optimal `LIWORK`.
- **`LIWORK`** — `LIWORK` is `INTEGER`. The dimension of the array `IWORK`.
- If `JOBZ = 'N'` or `N <= 1`, `LIWORK >= 1`.
- If `JOBZ = 'V'` and `N > 1`, `LIWORK >= 3 + 5*N`.
If `LIWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`.
- **`INFO`** — `INFO` is `INTEGER`
- `= 0`: successful exit 🎉
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value
- `> 0`: if `INFO = i`, and `i` is:
- `<= N`: the algorithm failed to converge: `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero;
- `> N`: if `INFO = N + i`, for `1 <= i <= N`, then `SPBSTF` returned `INFO = i`: `B` is not positive definite. The factorization of `B` could not be completed and no eigenvalues or eigenvectors were computed.
### ✨ Author ✨
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
### ✨ Contributors ✨
- Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
---
subroutine ssbgvx (character JOBZ, character RANGE, character UPLO, integer N, integer KA, integer KB, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldbb, * ) BB, integer LDBB, real, dimension( ldq, * ) Q, integer LDQ, real VL, real VU, integer IL, integer IU, real ABSTOL, integer M, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer, dimension( * ) IFAIL, integer INFO) SSBGVX
### ✨ Purpose
- **`JOBZ`** — <cyan>`CHARACTER*1`</cyan> 💖 okay so what do you actually want out of this routine, bestie?
- `'N'`: eigenvalues only, that's the whole order, nothing extra.
- `'V'`: eigenvalues **and** eigenvectors, we're going full main character 👑
- **`RANGE`** — <cyan>`CHARACTER*1`</cyan> this one decides *which* eigenvalues get pulled up ✨
- `'A'`: all eigenvalues will be found. every single one. no cap.
- `'V'`: all eigenvalues in the half-open interval `(VL,VU]` will be found.
- `'I'`: the `IL`-th through `IU`-th eigenvalues will be found.
- **`UPLO`** — <cyan>`CHARACTER*1`</cyan> which triangle of `A` and `B` did you actually store? 👀
- `'U'`: Upper triangles of `A` and `B` are stored;
- `'L'`: Lower triangles of `A` and `B` are stored.
- **`N`** — <cyan>`INTEGER`</cyan> the order of the matrices `A` and `B`. <hotpink>`N >= 0`</hotpink>.
- **`KA`** — <cyan>`INTEGER`</cyan> the number of superdiagonals of the matrix `A` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. <hotpink>`KA >= 0`</hotpink>.
- **`KB`** — <cyan>`INTEGER`</cyan> the number of superdiagonals of the matrix `B` if `UPLO = 'U'`, or the number of subdiagonals if `UPLO = 'L'`. <hotpink>`KB >= 0`</hotpink>.
- **`AB`** — <cyan>`REAL array, dimension (LDAB, N)`</cyan> 🦋 on entry, the upper or lower triangle of the symmetric band matrix `A`, stored in the first `ka+1` rows of the array. The `j`-th column of `A` is stored in the `j`-th column of the array `AB` as follows:
if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
On exit, the contents of `AB` are destroyed. 😭 she's gone, she served her purpose, let her rest.
- **`LDAB`** — <cyan>`INTEGER`</cyan> the leading dimension of the array `AB`. <hotpink>`LDAB >= KA+1`</hotpink>.
- **`BB`** — <cyan>`REAL array, dimension (LDBB, N)`</cyan> on entry, the upper or lower triangle of the symmetric band matrix `B`, stored in the first `kb+1` rows of the array. The `j`-th column of `B` is stored in the `j`-th column of the array `BB` as follows:
if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb).
On exit, the factor `S` from the split Cholesky factorization `B = S**T*S`, as returned by `SPBSTF`. 🫶
- **`LDBB`** — <cyan>`INTEGER`</cyan> the leading dimension of the array `BB`. <hotpink>`LDBB >= KB+1`</hotpink>.
- **`Q`** — <cyan>`REAL array, dimension (LDQ, N)`</cyan> if `JOBZ = 'V'`, the `n`-by-`n` matrix used in the reduction of `A*x = (lambda)*B*x` to standard form, i.e. `C*x = (lambda)*x`, and consequently `C` to tridiagonal form. If `JOBZ = 'N'`, the array `Q` is not referenced. she's off duty, don't @ her 🙈
- **`LDQ`** — <cyan>`INTEGER`</cyan> the leading dimension of the array `Q`. If `JOBZ = 'N'`, <hotpink>`LDQ >= 1`</hotpink>. If `JOBZ = 'V'`, <hotpink>`LDQ >= max(1,N)`</hotpink>.
- **`VL`** — <cyan>`REAL`</cyan> if `RANGE='V'`, the lower bound of the interval to be searched for eigenvalues. <hotpink>`VL < VU`</hotpink>. Not referenced if `RANGE = 'A'` or `'I'`.
- **`VU`** — <cyan>`REAL`</cyan> if `RANGE='V'`, the upper bound of the interval to be searched for eigenvalues. <hotpink>`VL < VU`</hotpink>. Not referenced if `RANGE = 'A'` or `'I'`.
- **`IL`** — <cyan>`INTEGER`</cyan> if `RANGE='I'`, the index of the smallest eigenvalue to be returned. <hotpink>`1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`</hotpink>. Not referenced if `RANGE = 'A'` or `'V'`.
- **`IU`** — <cyan>`INTEGER`</cyan> if `RANGE='I'`, the index of the largest eigenvalue to be returned. <hotpink>`1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`</hotpink>. Not referenced if `RANGE = 'A'` or `'V'`.
- **`ABSTOL`** — <cyan>`REAL`</cyan> the absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval `[a,b]` of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where `EPS` is the machine precision. If `ABSTOL` is less than or equal to zero, then `EPS*|T|` will be used in its place, where `|T|` is the 1-norm of the tridiagonal matrix obtained by reducing `A` to tridiagonal form.
Eigenvalues will be computed most accurately when `ABSTOL` is set to twice the underflow threshold `2*SLAMCH('S')`, not zero. If this routine returns with `INFO>0`, indicating that some eigenvectors did not converge, try setting `ABSTOL` to `2*SLAMCH('S')`.
> 🚨 bestie listen — set it to `2*SLAMCH('S')`, NOT zero. zero is not the vibe here, that's the accuracy you're leaving on the table. iykyk.
- **`M`** — <cyan>`INTEGER`</cyan> the total number of eigenvalues found. <hotpink>`0 <= M <= N`</hotpink>. If `RANGE = 'A'`, `M = N`, and if `RANGE = 'I'`, `M = IU-IL+1`.
- **`W`** — <cyan>`REAL array, dimension (N)`</cyan> if `INFO = 0`, the eigenvalues in ascending order. sorted, organized, iconic 💅
- **`Z`** — <cyan>`REAL array, dimension (LDZ, N)`</cyan> if `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the matrix `Z` of eigenvectors, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`. The eigenvectors are normalized so `Z**T*B*Z = I`. If `JOBZ = 'N'`, then `Z` is not referenced. ✨
- **`LDZ`** — <cyan>`INTEGER`</cyan> the leading dimension of the array `Z`. <hotpink>`LDZ >= 1`</hotpink>, and if `JOBZ = 'V'`, <hotpink>`LDZ >= max(1,N)`</hotpink>.
- **`WORK`** — <cyan>`REAL array, dimension (7*N)`</cyan> 🎀 the workspace. give her room, she needs `7*N`.
- **`IWORK`** — <cyan>`INTEGER array, dimension (5*N)`</cyan> 💜 integer workspace, `5*N`, no negotiation.
- **`IFAIL`** — <cyan>`INTEGER array, dimension (M)`</cyan> if `JOBZ = 'V'`, then if `INFO = 0`, the first `M` elements of `IFAIL` are zero. If `INFO > 0`, then `IFAIL` contains the indices of the eigenvalues that failed to converge. If `JOBZ = 'N'`, then `IFAIL` is not referenced.
- **`INFO`** — <cyan>`INTEGER`</cyan> the receipts 📋
- `= 0`: successful exit. slay.
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value. whoopsie.
- `<= N`: if `INFO = i`, then `i` eigenvectors failed to converge. Their indices are stored in `IFAIL`. 😭
- `> N`: `SPBSTF` returned an error code; i.e., if `INFO = N + i`, for `1 <= i <= N`, then the leading minor of order `i` of `B` is not positive definite. The factorization of `B` could not be completed and no eigenvalues or eigenvectors were computed.
> 💀 that last one is the "we got nothing done today" ending and i'm so serious about it — no eigenvalues, no eigenvectors, nada.
---
### 💖 Author 💖
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
### 🌟 Contributors 🌟
- Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
---
## 🦋 subroutine sspev 🦋
subroutine sspev (character JOBZ, character UPLO, integer N, real, dimension( * ) AP, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)
SSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices. 🫶 (yes, "OTHER matrices", that's literally what it says, i don't make the rules bestie)
### 🎯 Purpose 🎯
SSPEV computes all the eigenvalues and, optionally, eigenvectors of a real symmetric matrix `A` in packed storage. 🍓
### ✨ Parameters ✨
- **`JOBZ`** — <cyan>`CHARACTER*1`</cyan>
- `'N'`: Compute eigenvalues only;
- `'V'`: Compute eigenvalues and eigenvectors.
- **`UPLO`** — <cyan>`CHARACTER*1`</cyan>
- `'U'`: Upper triangle of `A` is stored;
- `'L'`: Lower triangle of `A` is stored.
- **`N`** — <cyan>`INTEGER`</cyan> the order of the matrix `A`. <hotpink>`N >= 0`</hotpink>.
- **`AP`** — <cyan>`REAL array, dimension (N*(N+1)/2)`</cyan> on entry, the upper or lower triangle of the symmetric matrix `A`, packed columnwise in a linear array. The `j`-th column of `A` is stored in the array `AP` as follows:
if UPLO = 'U', AP(i + (j-1)j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)(2*n-j)/2) = A(i,j) for j<=i<=n.
On exit, `AP` is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the diagonal and first superdiagonal of the tridiagonal matrix `T` overwrite the corresponding elements of `A`, and if `UPLO = 'L'`, the diagonal and first subdiagonal of `T` overwrite the corresponding elements of `A`. 🌀
- **`W`** — <cyan>`REAL array, dimension (N)`</cyan> if `INFO = 0`, the eigenvalues in ascending order. ✨
- **`Z`** — <cyan>`REAL array, dimension (LDZ, N)`</cyan> if `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the orthonormal eigenvectors of the matrix `A`, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`. If `JOBZ = 'N'`, then `Z` is not referenced. 🙈
- **`LDZ`** — <cyan>`INTEGER`</cyan> the leading dimension of the array `Z`. <hotpink>`LDZ >= 1`</hotpink>, and if `JOBZ = 'V'`, <hotpink>`LDZ >= max(1,N)`</hotpink>.
- **`WORK`** — <cyan>`REAL array, dimension (3*N)`</cyan> 🧚 the workspace, `3*N` of it.
- **`INFO`** — <cyan>`INTEGER`</cyan> the receipts 📋
- `= 0`: successful exit.
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value.
- `> 0`: if `INFO = i`, the algorithm failed to converge; `i` off-diagonal elements of an intermediate tridiagonal form did not converge to zero.
### 💜 Author 💜
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
## 🧚 subroutine sspevd 🧚
subroutine sspevd (character JOBZ, character UPLO, integer N, real, dimension( * ) AP, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices.
### 🎀 Purpose 🎀
> 🌈 annnnd that's exactly where this part runs out — the rest of SSPEVD's purpose + all its parameters are rolling in with the next chunk, keep scrolling bestie 💖
<cyan>`SSPEVD`</cyan> computes all the eigenvalues and, optionally, the eigenvectors of a real symmetric matrix <purple>A</purple> in packed storage. If eigenvectors are desired, it uses a divide and conquer algorithm. Bestie, it eats. 🔥
The divide and conquer algorithm makes <yellow>very mild</yellow> assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. (Not the Crays catching strays 💀 but okay, noted.)
### Parameters
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: Upper triangle of A is stored;
- `= 'L'`: Lower triangle of A is stored.
- **`N`** — `N` is `INTEGER`
- The order of the matrix A. `N >= 0`.
- **`AP`** — `AP` is `REAL` array, dimension `(N*(N+1)/2)`
- On entry, the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array `AP` as follows:
- if `UPLO = 'U'`, `AP(i + (j-1)*j/2) = A(i,j)` for `1<=i<=j`;
- if `UPLO = 'L'`, `AP(i + (j-1)*(2*n-j)/2) = A(i,j)` for `j<=i<=n`.
- On exit, `AP` is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the diagonal and first superdiagonal of the tridiagonal matrix T overwrite the corresponding elements of A, and if `UPLO = 'L'`, the diagonal and first subdiagonal of T overwrite the corresponding elements of A. <peach>(yes bestie your input is gone, it's giving sacrifice-for-science 🫶)</peach>
- **`W`** — `W` is `REAL` array, dimension `(N)`
- If `INFO = 0`, the eigenvalues in ascending order.
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, N)`
- If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the orthonormal eigenvectors of the matrix A, with the i-th column of `Z` holding the eigenvector associated with `W(i)`.
- If `JOBZ = 'N'`, then `Z` is not referenced.
- **`LDZ`** — `LDZ` is `INTEGER`
- The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL` array, dimension `(MAX(1,LWORK))`
- On exit, if `INFO = 0`, `WORK(1)` returns the required `LWORK`.
- **`LWORK`** — `LWORK` is `INTEGER`
- The dimension of the array `WORK`.
- If `N <= 1`, `LWORK` must be at least `1`.
- If `JOBZ = 'N'` and `N > 1`, `LWORK` must be at least `2*N`.
- If `JOBZ = 'V'` and `N > 1`, `LWORK` must be at least `1 + 6*N + N**2`.
- If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the required sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`. <mint>(the little peek-behind-the-curtain move, iconic ✨)</mint>
- **`IWORK`** — `IWORK` is `INTEGER` array, dimension `(MAX(1,LIWORK))`
- On exit, if `INFO = 0`, `IWORK(1)` returns the required `LIWORK`.
- **`LIWORK`** — `LIWORK` is `INTEGER`
- The dimension of the array `IWORK`.
- If `JOBZ = 'N'` or `N <= 1`, `LIWORK` must be at least `1`.
- If `JOBZ = 'V'` and `N > 1`, `LIWORK` must be at least `3 + 5*N`.
- If `LIWORK = -1`, then a workspace query is assumed; the routine only calculates the required sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`.
- **`INFO`** — `INFO` is `INTEGER`
- `= 0`: successful exit
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value.
- `> 0`: if `INFO = i`, the algorithm failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero. <red>(sad girl hours, it tried 😭)</red>
### 👑 Author 👑
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
## 🦋 subroutine `sspevx` 🦋
subroutine sspevx (character JOBZ, character RANGE, character UPLO, integer N, real, dimension( * ) AP, real VL, real VU, integer IL, integer IU, real ABSTOL, integer M, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer, dimension( * ) IFAIL, integer INFO)
`SSPEVX` computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
**Purpose:**
<cyan>`SSPEVX`</cyan> computes <yellow>selected</yellow> eigenvalues and, optionally, eigenvectors of a real symmetric matrix A in packed storage. Eigenvalues/vectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues. <lavender>(picky queen energy, and honestly? valid 💅)</lavender>
### Parameters
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`RANGE`** — `RANGE` is `CHARACTER*1`
- `= 'A'`: all eigenvalues will be found;
- `= 'V'`: all eigenvalues in the half-open interval `(VL,VU]` will be found;
- `= 'I'`: the IL-th through IU-th eigenvalues will be found.
- **`UPLO`** — `UPLO` is `CHARACTER*1`
- `= 'U'`: Upper triangle of A is stored;
- `= 'L'`: Lower triangle of A is stored.
- **`N`** — `N` is `INTEGER`
- The order of the matrix A. `N >= 0`.
- **`AP`** — `AP` is `REAL` array, dimension `(N*(N+1)/2)`
- On entry, the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array `AP` as follows:
- if `UPLO = 'U'`, `AP(i + (j-1)*j/2) = A(i,j)` for `1<=i<=j`;
- if `UPLO = 'L'`, `AP(i + (j-1)*(2*n-j)/2) = A(i,j)` for `j<=i<=n`.
- On exit, `AP` is overwritten by values generated during the reduction to tridiagonal form. If `UPLO = 'U'`, the diagonal and first superdiagonal of the tridiagonal matrix T overwrite the corresponding elements of A, and if `UPLO = 'L'`, the diagonal and first subdiagonal of T overwrite the corresponding elements of A.
- **`VL`** — `VL` is `REAL`
- If `RANGE='V'`, the lower bound of the interval to be searched for eigenvalues. `VL < VU`.
- Not referenced if `RANGE = 'A'` or `'I'`.
- **`VU`** — `VU` is `REAL`
- If `RANGE='V'`, the upper bound of the interval to be searched for eigenvalues. `VL < VU`.
- Not referenced if `RANGE = 'A'` or `'I'`.
- **`IL`** — `IL` is `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`** — `IU` is `INTEGER`
- If `RANGE='I'`, the index of the largest eigenvalue to be returned.
- `1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`.
- Not referenced if `RANGE = 'A'` or `'V'`.
- **`ABSTOL`** — `ABSTOL` is `REAL`
- The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval `[a,b]` of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
- where `EPS` is the machine precision. If `ABSTOL` is less than or equal to zero, then `EPS*|T|` will be used in its place, where `|T|` is the 1-norm of the tridiagonal matrix obtained by reducing `AP` to tridiagonal form.
- Eigenvalues will be computed most accurately when `ABSTOL` is set to twice the underflow threshold `2*SLAMCH('S')`, <hotpink>not zero</hotpink>. If this routine returns with `INFO>0`, indicating that some eigenvectors did not converge, try setting `ABSTOL` to `2*SLAMCH('S')`. <peach>(the glow-up tip of the century fr fr 🎀)</peach>
- See 'Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy,' by Demmel and Kahan, LAPACK Working Note #3.
- **`M`** — `M` is `INTEGER`
- The total number of eigenvalues found. `0 <= M <= N`.
- If `RANGE = 'A'`, `M = N`, and if `RANGE = 'I'`, `M = IU-IL+1`.
- **`W`** — `W` is `REAL` array, dimension `(N)`
- If `INFO = 0`, the selected eigenvalues in ascending order.
- **`Z`** — `Z` is `REAL` array, dimension `(LDZ, max(1,M))`
- If `JOBZ = 'V'`, then if `INFO = 0`, the first M columns of `Z` contain the orthonormal eigenvectors of the matrix A corresponding to the selected eigenvalues, with the i-th column of `Z` holding the eigenvector associated with `W(i)`.
- If an eigenvector fails to converge, then that column of `Z` contains the latest approximation to the eigenvector, and the index of the eigenvector is returned in `IFAIL`.
- If `JOBZ = 'N'`, then `Z` is not referenced.
- Note: the user must ensure that at least `max(1,M)` columns are supplied in the array `Z`; if `RANGE = 'V'`, the exact value of M is not known in advance and an upper bound must be used. <yellow>(so like, overbook the venue bestie, no cap 🫶)</yellow>
- **`LDZ`** — `LDZ` is `INTEGER`
- The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL` array, dimension `(8*N)`
- **`IWORK`** — `IWORK` is `INTEGER` array, dimension `(5*N)`
- **`IFAIL`** — `IFAIL` is `INTEGER` array, dimension `(N)`
- If `JOBZ = 'V'`, then if `INFO = 0`, the first M elements of `IFAIL` are zero. If `INFO > 0`, then `IFAIL` contains the indices of the eigenvectors that failed to converge.
- If `JOBZ = 'N'`, then `IFAIL` is not referenced.
- **`INFO`**
- **`INFO`** — an INTEGER, the receipts:
= 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
0: if INFO = i, then i eigenvectors failed to converge.
Their indices are stored in array IFAIL.
So if `INFO` comes back **positive**, bestie, some eigenvectors just refused to converge 😭 and their little index numbers are chilling in `IFAIL` so you can go find them.
## 👑 Author 👑
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
---
### 🌸 subroutine `sspgv` 🌸
subroutine sspgv (integer ITYPE, character JOBZ, character UPLO, integer N, real, dimension( * ) AP, real, dimension( * ) BP, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)
**`SSPGV`** — say it with me 💅
#### ✨ Purpose ✨
<pink>`SSPGV` computes **all** the eigenvalues and, optionally, the eigenvectors of a real generalized symmetric-definite eigenproblem</pink>, of the form
A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
Here `A` and `B` are assumed to be symmetric, stored in **packed format**, and `B` is also **positive definite** 🫶
#### 🌈 Parameters 🌈
- **`ITYPE`** — INTEGER, tells us which flavor of problem we're solving:
= 1: A*x = (lambda)*B*x = 2: A*B*x = (lambda)*x = 3: B*A*x = (lambda)*x
- **`JOBZ`** — CHARACTER*1, the "do I want the eigenvectors or nah" switch:
= 'N': Compute eigenvalues only; = 'V': Compute eigenvalues and eigenvectors.
- **`UPLO`** — CHARACTER*1, which triangle of `A` and `B` you actually stored:
= 'U': Upper triangles of A and B are stored; = 'L': Lower triangles of A and B are stored.
- **`N`** — INTEGER, the order of the matrices `A` and `B`. `N >= 0`.
- **`AP`** — REAL array, dimension `(N*(N+1)/2)`. On entry, it holds the upper or lower triangle of the symmetric matrix `A`, packed columnwise in a linear array 🔥 The j-th column of `A` lives in `AP` like this:
if UPLO = 'U', AP(i + (j-1)j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)(2*n-j)/2) = A(i,j) for j<=i<=n.
> On exit, the contents of `AP` are **destroyed**. Gone. Not coming back 💀
- **`BP`** — REAL array, dimension `(N*(N+1)/2)`. On entry, same deal but for the symmetric matrix `B`, packed columnwise in a linear array:
if UPLO = 'U', BP(i + (j-1)j/2) = B(i,j) for 1<=i<=j; if UPLO = 'L', BP(i + (j-1)(2*n-j)/2) = B(i,j) for j<=i<=n.
On exit, `BP` holds the triangular factor `U` or `L` from the Cholesky factorization `B = U**T*U` or `B = L*L**T`, in the same storage format as `B` ✨
- **`W`** — REAL array, dimension `(N)`. If `INFO = 0`, this is your eigenvalues, in **ascending order** 👑
- **`Z`** — REAL array, dimension `(LDZ, N)`. If `JOBZ = 'V'`, then when `INFO = 0`, `Z` contains the matrix `Z` of eigenvectors. The eigenvectors are normalized as follows:
if ITYPE = 1 or 2, ZT*B*Z = I; if ITYPE = 3, ZT*inv(B)*Z = I.
If `JOBZ = 'N'`, then `Z` is **not referenced** — she's off duty 💅
- **`LDZ`** — INTEGER, the leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, then `LDZ >= max(1,N)`.
- **`WORK`** — REAL array, dimension `(3*N)`.
- **`INFO`** — INTEGER, the final vibe check:
= 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
0: SPPTRF or SSPEV returned an error code:
<= N: if INFO = i, SSPEV failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero.
N: if INFO = n + i, for 1 <= i <= n, then the leading
minor of order i of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed.
> Bestie, if you get `INFO > N` it's giving "`B` wasn't positive definite" and literally **no** eigenvalues or eigenvectors got computed 😭
## 👑 Author 👑
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
---
### 🦋 subroutine `sspgvd` 🦋
subroutine sspgvd (integer ITYPE, character JOBZ, character UPLO, integer N, real, dimension( * ) AP, real, dimension( * ) BP, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
**`SSPGVD`** ✨
#### ✨ Purpose ✨
<pink>`SSPGVD` computes all the eigenvalues, and optionally, the eigenvectors of a real generalized symmetric-definite eigenproblem</pink>, of the form
A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
Here `A` and `B` are assumed to be symmetric, stored in **packed format**, and `B` is also **positive definite**.
If eigenvectors are desired, it uses a **divide and conquer algorithm** 🧚
The divide and conquer algorithm makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none.
#### 🌈 Parameters 🌈
- **`ITYPE`** — INTEGER, which problem type are we solving:
= 1: A*x = (lambda)*B*x = 2: A*B*x = (lambda)*x = 3: B*A*x = (lambda)*x
- **`JOBZ`** — CHARACTER*1:
= 'N': Compute eigenvalues only; = 'V': Compute eigenvalues and eigenvectors.
- **`UPLO`** — CHARACTER*1:
= 'U': Upper triangles of A and B are stored; = 'L': Lower triangles of A and B are stored.
- **`N`** — INTEGER, the order of the matrices `A` and `B`. `N >= 0`.
- **`AP`** — REAL array, dimension `(N*(N+1)/2)`. On entry, the upper or lower triangle of the symmetric matrix `A`, packed columnwise in a linear array. The j-th column of `A` is stored in `AP` as follows:
if UPLO = 'U', AP(i + (j-1)j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)(2*n-j)/2) = A(i,j) for j<=i<=n.
> On exit, the contents of `AP` are **destroyed** 💀
- **`BP`** — REAL array, dimension `(N*(N+1)/2)`. On entry, the upper or lower triangle of the symmetric matrix `B`, packed columnwise in a linear array. The j-th column of `B` is stored in `BP` as follows:
if UPLO = 'U', BP(i + (j-1)j/2) = B(i,j) for 1<=i<=j; if UPLO = 'L', BP(i + (j-1)(2*n-j)/2) = B(i,j) for j<=i<=n.
On exit, the triangular factor `U` or `L` from the Cholesky factorization `B = U**T*U` or `B = L*L**T`, in the same storage format as `B` 🌸
- **`W`** — REAL array, dimension `(N)`. If `INFO = 0`, the eigenvalues in ascending order ⭐
- **`Z`** — REAL array, dimension `(LDZ, N)`. If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the matrix `Z` of eigenvectors. The eigenvectors are normalized as follows:
if ITYPE = 1 or 2, ZT*B*Z = I; if ITYPE = 3, ZT*inv(B)*Z = I.
If `JOBZ = 'N'`, then `Z` is not referenced.
- **`LDZ`** — INTEGER, the leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — REAL array, dimension `(MAX(1,LWORK))`. On exit, if `INFO = 0`, `WORK(1)` returns the required `LWORK`.
- **`LWORK`** — INTEGER, the dimension of the array `WORK`:
If N <= 1, LWORK >= 1. If JOBZ = 'N' and N > 1, LWORK >= 2*N. If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
If `LWORK = -1`, then a **workspace query** is assumed; the routine only calculates the required sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and **no error message related to `LWORK` or `LIWORK` is issued by XERBLA** 👀
- **`IWORK`** — INTEGER array, dimension `(MAX(1,LIWORK))`. On exit, if `INFO = 0`, `IWORK(1)` returns the required `LIWORK`.
- **`LIWORK`** — INTEGER, the dimension of the array `IWORK`:
If JOBZ = 'N' or N <= 1, LIWORK >= 1. If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
If `LIWORK = -1`, then a workspace query is assumed; the routine only calculates the required sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`.
- **`INFO`** — INTEGER:
= 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value
0: SPPTRF or SSPEVD returned an error code:
<= N: if INFO = i, SSPEVD failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero;
N: if INFO = N + i, for 1 <= i <= N, then the leading
minor of order i of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed.
## 👑 Author 👑
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
## 🫶 Contributors 🫶
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
---
### 💜 subroutine `sspgvx` 💜
subroutine sspgvx (integer ITYPE, character JOBZ, character RANGE, character UPLO, integer N, real, dimension( * ) AP, real, dimension( * ) BP, real VL, real VU, integer IL, integer IU, real ABSTOL, integer M, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer, dimension(
- ) IFAIL, integer INFO)
**`SSPGVX`** 🔥
#### ✨ Purpose ✨
<pink>okay so `SSPGVX`</pink> is the girl who hunts down <hotpink>*selected*</hotpink> eigenvalues — and, if you ask nicely, the eigenvectors too — for a real generalized symmetric-definite eigenproblem, and the three fits she serves are:
A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x
A and B are assumed to be <cyan>symmetric</cyan>, stored in <cyan>packed storage</cyan>, and B is <hotpink>also positive definite</hotpink> (non-negotiable, bestie 💅). You get to pick which eigenvalues show up by giving either a <lavender>range of values</lavender> or a <lavender>range of indices</lavender> for the ones you want. Choose your fighter ✨
---
## 💖 Parameters 💖
- **`ITYPE`** — <cyan>INTEGER</cyan>. The ✨vibe selector✨ that says which problem type we're solving:
- `= 1`: `A*x = (lambda)*B*x`
- `= 2`: `A*B*x = (lambda)*x`
- `= 3`: `B*A*x = (lambda)*x`
- **`JOBZ`** — <cyan>CHARACTER*1</cyan>. Do you want just the numbers or the whole fit? 👑
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues <hotpink>and</hotpink> eigenvectors.
- **`RANGE`** — <cyan>CHARACTER*1</cyan>. Which ones are we pulling up?
- `= 'A'`: all eigenvalues will be found.
- `= 'V'`: all eigenvalues in the half-open interval `(VL,VU]` will be found.
- `= 'I'`: the `IL`-th through `IU`-th eigenvalues will be found.
- **`UPLO`** — <cyan>CHARACTER*1</cyan>. Which half of A and B are we even looking at? 🙈
- `= 'U'`: Upper triangle of A and B are stored;
- `= 'L'`: Lower triangle of A and B are stored.
- **`N`** — <cyan>INTEGER</cyan>. The order of the matrix pencil `(A,B)`. `N >= 0`.
- **`AP`** — <cyan>REAL array, dimension (N*(N+1)/2)</cyan>. On entry, the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The `j`-th column of A is stored in the array `AP` as follows:
- if `UPLO = 'U'`, `AP(i + (j-1)*j/2) = A(i,j)` for `1<=i<=j`;
- if `UPLO = 'L'`, `AP(i + (j-1)*(2*n-j)/2) = A(i,j)` for `j<=i<=n`.
On exit, <red>the contents of `AP` are destroyed</red>. 💀 she does not survive the trip, sorry.
- **`BP`** — <cyan>REAL array, dimension (N*(N+1)/2)</cyan>. On entry, the upper or lower triangle of the symmetric matrix B, packed columnwise in a linear array. The `j`-th column of B is stored in the array `BP` as follows:
- if `UPLO = 'U'`, `BP(i + (j-1)*j/2) = B(i,j)` for `1<=i<=j`;
- if `UPLO = 'L'`, `BP(i + (j-1)*(2*n-j)/2) = B(i,j)` for `j<=i<=n`.
On exit, the triangular factor `U` or `L` from the Cholesky factorization `B = U**T*U` or `B = L*L**T`, in the same storage format as B. ✨glow up✨
- **`VL`** — <cyan>REAL</cyan>. If `RANGE='V'`, the lower bound of the interval to be searched for eigenvalues. `VL < VU`. <lavender>Not referenced if `RANGE = 'A'` or `'I'`.</lavender>
- **`VU`** — <cyan>REAL</cyan>. If `RANGE='V'`, the upper bound of the interval to be searched for eigenvalues. `VL < VU`. <lavender>Not referenced if `RANGE = 'A'` or `'I'`.</lavender>
- **`IL`** — <cyan>INTEGER</cyan>. 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`. <lavender>Not referenced if `RANGE = 'A'` or `'V'`.</lavender>
- **`IU`** — <cyan>INTEGER</cyan>. 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`. <lavender>Not referenced if `RANGE = 'A'` or `'V'`.</lavender>
- **`ABSTOL`** — <cyan>REAL</cyan>. The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval `[a,b]` of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where `EPS` is the machine precision. If `ABSTOL` is less than or equal to zero, then `EPS*|T|` will be used in its place, where `|T|` is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when `ABSTOL` is set to twice the underflow threshold `2*SLAMCH('S')`, <hotpink>not zero</hotpink>. If this routine returns with `INFO>0`, indicating that some eigenvectors did not converge, try setting `ABSTOL` to `2*SLAMCH('S')`. 🫶 the fix is right there, take it.
- **`M`** — <cyan>INTEGER</cyan>. The total number of eigenvalues found. `0 <= M <= N`. If `RANGE = 'A'`, `M = N`, and if `RANGE = 'I'`, `M = IU-IL+1`.
- **`W`** — <cyan>REAL array, dimension (N)</cyan>. On normal exit, the first `M` elements contain the selected eigenvalues in ascending order. ✨sorted, organized, thriving✨
- **`Z`** — <cyan>REAL array, dimension (LDZ, max(1,M))</cyan>. If `JOBZ = 'N'`, then `Z` is not referenced. If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` columns of `Z` contain the orthonormal eigenvectors of the matrix A corresponding to the selected eigenvalues, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`. The eigenvectors are normalized as follows:
- if `ITYPE = 1` or `2`, `Z**T*B*Z = I`;
- if `ITYPE = 3`, `Z**T*inv(B)*Z = I`.
If an eigenvector fails to converge, then that column of `Z` contains the latest approximation to the eigenvector, and the index of the eigenvector is returned in `IFAIL`. 😭
> Note: the user must ensure that at least `max(1,M)` columns are supplied in the array `Z`; if `RANGE = 'V'`, the exact value of `M` is not known in advance and an <hotpink>upper bound must be used</hotpink>. no cap, size it big.
- **`LDZ`** — <cyan>INTEGER</cyan>. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — <cyan>REAL array, dimension (8*N)</cyan>
- **`IWORK`** — <cyan>INTEGER array, dimension (5*N)</cyan>
- **`IFAIL`** — <cyan>INTEGER array, dimension (N)</cyan>. If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` elements of `IFAIL` are zero. If `INFO > 0`, then `IFAIL` contains the indices of the eigenvectors that failed to converge. If `JOBZ = 'N'`, then `IFAIL` is not referenced.
- **`INFO`** — <cyan>INTEGER</cyan>. The receipts 📜:
- `= 0`: successful exit
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value
- `> 0`: `SPPTRF` or `SSPEVX` returned an error code:
- `<= N`: if `INFO = i`, `SSPEVX` failed to converge; `i` eigenvectors failed to converge. Their indices are stored in array `IFAIL`.
- `> N`: if `INFO = N + i`, for `1 <= i <= N`, then the leading minor of order `i` of B is not positive definite. The factorization of B could not be completed and <red>no eigenvalues or eigenvectors were computed</red>. 💀
---
## 🌈 Author 🌈
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
## 🍓 Contributors 🍓
- Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
---
## 🦋 subroutine sstev 🦋
subroutine sstev (character JOBZ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)
<pink>SSTEV</pink> computes the eigenvalues and, optionally, the left and/or right eigenvectors for <hotpink>OTHER</hotpink> matrices. 💅
### ✨ Purpose: ✨
`SSTEV` computes <cyan>all</cyan> eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix A.
### 🎀 Parameters 🎀
- **`JOBZ`** — <cyan>CHARACTER*1</cyan>. Pick your energy:
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`N`** — <cyan>INTEGER</cyan>. The order of the matrix. `N >= 0`.
- **`D`** — <cyan>REAL array, dimension (N)</cyan>. On entry, the `n` diagonal elements of the tridiagonal matrix A. On exit, if `INFO = 0`, the eigenvalues in ascending order. ✨
- **`E`** — <cyan>REAL array, dimension (N-1)</cyan>. On entry, the `(n-1)` subdiagonal elements of the tridiagonal matrix A, stored in elements `1` to `N-1` of `E`. On exit, <red>the contents of `E` are destroyed</red>. 💀 gone, reduced to atoms.
- **`Z`** — <cyan>REAL array, dimension (LDZ, N)</cyan>. If `JOBZ = 'V'`, then if `INFO = 0`, `Z` contains the orthonormal eigenvectors of the matrix A, with the `i`-th column of `Z` holding the eigenvector associated with `D(i)`. If `JOBZ = 'N'`, then `Z` is not referenced.
- **`LDZ`** — <cyan>INTEGER</cyan>. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — <cyan>REAL array, dimension (max(1,2*N-2))</cyan>. If `JOBZ = 'N'`, `WORK` is not referenced.
- **`INFO`** — <cyan>INTEGER</cyan>. The tea ☕:
- `= 0`: successful exit
- `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value
- `> 0`: if `INFO = i`, the algorithm failed to converge; `i` off-diagonal elements of `E` did not converge to zero. 😭
### 💜 Author 💜
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
---
## 🌸 subroutine sstevd 🌸
subroutine sstevd (character JOBZ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
<pink>SSTEVD</pink> computes the eigenvalues and, optionally, the left and/or right eigenvectors for <hotpink>OTHER</hotpink> matrices. 💅
### ✨ Purpose: ✨
`SSTEVD` computes <cyan>all</cyan> eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix. If eigenvectors are desired, it uses a <hotpink>divide and conquer algorithm</hotpink>. 👑
The divide and conquer algorithm makes <lavender>very mild assumptions</lavender> about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. 🦋
### 🎀 Parameters 🎀
- **`JOBZ`** — <cyan>CHARACTER*1</cyan>. The eternal question:
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`N`** — <cyan>INTEGER</cyan>. The order of the matrix. `N >= 0`.
- **`D`**
- **`D`** — REAL array, dimension `(N)`. On the way *in*, bestie, `D` is holding the `n` diagonal elements of your tridiagonal matrix `A`. On the way *out*, if `INFO = 0`, those babies are your eigenvalues in **ascending order**. 💖
- **`E`** — REAL array, dimension `(N-1)`. Coming in, it's the `(n-1)` subdiagonal elements of the tridiagonal matrix `A`, stored in elements `1` to `N-1` of `E`. Going out? The contents of `E` are **destroyed**. Gone. Poof. 😭 So don't come running back for them, they don't live here anymore.
- **`Z`** — REAL array, dimension `(LDZ, N)`. If `JOBZ = 'V'`, then when `INFO = 0`, `Z` contains the orthonormal eigenvectors of the matrix `A`, with the `i`-th column of `Z` holding the eigenvector that belongs to `D(i)`. If `JOBZ = 'N'`, then `Z` is simply **not referenced** — she's off duty, unbothered, thriving elsewhere. 💅
- **`LDZ`** — INTEGER. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, then `LDZ >= max(1,N)`.
- **`WORK`** — REAL array, dimension `(LWORK)`. On exit, if `INFO = 0`, `WORK(1)` hands you back the **optimal `LWORK`**. 🎀
- **`LWORK`** — INTEGER. The dimension of the array `WORK`.
- If `JOBZ = 'N'` or `N <= 1` then `LWORK` must be at least `1`.
- If `JOBZ = 'V'` and `N > 1` then `LWORK` must be at least `( 1 + 4*N + N**2 )`.
If `LWORK = -1`, then a **workspace query** is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and **no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`**. So you can ask "hey how much room do I need" without getting yelled at. Iconic. ✨
- **`IWORK`** — INTEGER array, dimension `(MAX(1,LIWORK))`. On exit, if `INFO = 0`, `IWORK(1)` returns the optimal `LIWORK`.
- **`LIWORK`** — INTEGER. The dimension of the array `IWORK`.
- If `JOBZ = 'N'` or `N <= 1` then `LIWORK` must be at least `1`.
- If `JOBZ = 'V'` and `N > 1` then `LIWORK` must be at least `3+5*N`.
If `LIWORK = -1`, then a **workspace query** is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and **no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`**. Same deal, different array, still no drama. 🫶
- **`INFO`** — INTEGER.
- `= 0`: successful exit 🎉
- `< 0`: if `INFO = -i`, the `i`-th argument had an **illegal value** 💀
- `> 0`: if `INFO = i`, the algorithm **failed to converge**; `i` off-diagonal elements of `E` did not converge to zero 😭
---
## 👑 Author 👑
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
---
### 🦋 subroutine `sstevr` 🦋
subroutine sstevr (character JOBZ, character RANGE, integer N, real, dimension( * ) D, real, dimension( * ) E, real VL, real VU, integer IL, integer IU, real ABSTOL, integer M, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
SSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for **OTHER** matrices. Yes, `OTHER`. That's literally what it says. 💅
### 🌸 Purpose 🌸
SSTEVR computes **selected** eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix `T`. And here's the fun part — you get to *choose* which ones! Eigenvalues and eigenvectors can be selected by specifying either a **range of values** or a **range of indices** for the eigenvalues you're craving. 🍓
Whenever possible, SSTEVR calls `SSTEMR` to compute the eigenspectrum using **Relatively Robust Representations**. `SSTEMR` computes eigenvalues by the **dqds algorithm**, while orthogonal eigenvectors are computed from various 'good' `L D L^T` representations (also known as Relatively Robust Representations). Gram-Schmidt orthogonalization is avoided as far as possible. More specifically, the various steps of the algorithm are as follows. For the `i`-th unreduced block of `T`,
- **(a)** Compute `T - sigma_i = L_i D_i L_i^T`, such that `L_i D_i L_i^T` is a relatively robust representation,
- **(b)** Compute the eigenvalues, `lambda_j`, of `L_i D_i L_i^T` to high relative accuracy by the dqds algorithm,
- **(c)** If there is a cluster of close eigenvalues, 'choose' `sigma_i` close to the cluster, and go to step (a),
- **(d)** Given the approximate eigenvalue `lambda_j` of `L_i D_i L_i^T`, compute the corresponding eigenvector by forming a rank-revealing twisted factorization.
The desired accuracy of the output can be specified by the input parameter `ABSTOL`.
For more details, see *'A new O(n^2) algorithm for the symmetric tridiagonal eigenvalue/eigenvector problem'*, by Inderjit Dhillon, Computer Science Division Technical Report No. `UCB//CSD-97-971`, UC Berkeley, May 1997.
> **Note 1** : SSTEVR calls `SSTEMR` when the **full spectrum** is requested on machines which conform to the ieee-754 floating point standard. SSTEVR calls `SSTEBZ` and `SSTEIN` on non-ieee machines and when **partial spectrum** requests are made.
> Normal execution of `SSTEMR` may create **NaNs and infinities** and hence may **abort** due to a floating point exception in environments which do not handle NaNs and infinities in the ieee standard default manner. So like… heads up. 🚨
### 🎀 Parameters 🎀
- **`JOBZ`** — CHARACTER*1.
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues **and** eigenvectors.
- **`RANGE`** — CHARACTER*1.
- `= 'A'`: all eigenvalues will be found.
- `= 'V'`: all eigenvalues in the half-open interval `(VL,VU]` will be found.
- `= 'I'`: the `IL`-th through `IU`-th eigenvalues will be found.
- For `RANGE = 'V'` or `'I'` and `IU - IL < N - 1`, `SSTEBZ` and `SSTEIN` are called
- **`N`** — INTEGER. The order of the matrix. `N >= 0`.
- **`D`** — REAL array, dimension `(N)`. On entry, the `n` diagonal elements of the tridiagonal matrix `A`. On exit, `D` **may be multiplied by a constant factor** chosen to avoid over/underflow in computing the eigenvalues. So don't be shocked when the numbers look scaled, it's on purpose bestie. 🙈
- **`E`** — REAL array, dimension `(max(1,N-1))`. On entry, the `(n-1)` subdiagonal elements of the tridiagonal matrix `A` in elements `1` to `N-1` of `E`. On exit, `E` **may be multiplied by a constant factor** chosen to avoid over/underflow in computing the eigenvalues. Same vibe, same disclaimer. 😌
- **`VL`** — REAL. If `RANGE='V'`, the lower bound of the interval to be searched for eigenvalues. `VL < VU`. Not referenced if `RANGE = 'A'` or `'I'`.
- **`VU`** — REAL. If `RANGE='V'`, the upper bound of the interval to be searched for eigenvalues. `VL < VU`. Not referenced if `RANGE = 'A'` or `'I'`.
- **`IL`** — INTEGER. If `RANGE='I'`, the index of the smallest eigenvalue to be returned. `1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`. Not referenced if `RANGE = 'A'` or `'V'`.
- **`IU`** — INTEGER. If `RANGE='I'`, the index of the largest eigenvalue to be returned. `1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`. Not referenced if `RANGE = 'A'` or `'V'`.
- **`ABSTOL`** — REAL. The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval `[a,b]` of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where `EPS` is the machine precision. If `ABSTOL` is less than or equal to zero, then `EPS*|T|` will be used in its place, where `|T|` is the 1-norm of the tridiagonal matrix obtained by reducing `A` to tridiagonal form.
See *'Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy,'* by Demmel and Kahan, LAPACK Working Note #3.
If high relative accuracy is important, set `ABSTOL` to `SLAMCH( 'Safe minimum' )`. Doing so will guarantee that eigenvalues are computed to high relative accuracy when possible **in future releases**. The current code does **not** make any guarantees about high relative accuracy, but future releases will. See J. Barlow and J. Demmel, *'Computing Accurate Eigensystems of Scaled Diagonally Dominant Matrices'*, LAPACK Working Note #7, for a discussion of which matrices define their eigenvalues to high relative accuracy. No cap, that's what the docs say. 🔥
- **`M`** — INTEGER. The total number of eigenvalues found. `0 <= M <= N`. If `RANGE = 'A'`, `M = N`, and if `RANGE = 'I'`, `M = IU-IL+1`.
- **`W`** — REAL array, dimension `(N)`. The first `M` elements contain the selected eigenvalues in **ascending order**. ⭐
- **`Z`** — REAL array, dimension `(LDZ, max(1,M) )`. If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` columns of `Z` contain the orthonormal eigenvectors of the matrix `A` corresponding to the selected eigenvalues, with the `i`-th column of `Z` holding the eigenvector associated with `W(i)`.
> **Note**: the user must ensure that at least `max(1,M)` columns are supplied in the array `Z`; if `RANGE = 'V'`, the exact value of `M` is **not known in advance** and an **upper bound must be used**. So like… allocate generously, bestie. She's unpredictable. 😭
- **`LDZ`** — INTEGER. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`ISUPPZ`** — INTEGER array, dimension `( 2*max(1,M) )`. The support of the eigenvectors in `Z`, i.e., the indices indicating the nonzero elements in `Z`. The `i`-th eigenvector is nonzero only in elements `ISUPPZ( 2*i-1 )` through `ISUPPZ( 2*i )`. Implemented only for `RANGE = 'A'` or `'I'` and `IU - IL = N - 1`
- **`WORK`** — *(and yes bestie, the `WORK` deets are right around the corner 💅)*
- **`WORK`** — <pink>REAL array</pink>, dimension `(MAX(1,LWORK))`. On exit, if `INFO = 0`, `WORK(1)` returns the optimal (and minimal) `LWORK`. <rainbow>It's giving "workspace knows its own size".</rainbow> 💅
- **`LWORK`** — <pink>INTEGER</pink>. The dimension of the array `WORK`. `LWORK >= 20*N`. If `LWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`. <cyan>So if you're just peeking at sizes, set it to `-1` and vibe.</cyan>
- **`IWORK`** — <pink>INTEGER array</pink>, dimension `(MAX(1,LIWORK))`. On exit, if `INFO = 0`, `IWORK(1)` returns the optimal (and minimal) `LIWORK`. <rainbow>Same energy, different array.</rainbow>
- **`LIWORK`** — <pink>INTEGER</pink>. The dimension of the array `IWORK`. `LIWORK >= 10*N`. If `LIWORK = -1`, then a workspace query is assumed; the routine only calculates the optimal sizes of the `WORK` and `IWORK` arrays, returns these values as the first entries of the `WORK` and `IWORK` arrays, and no error message related to `LWORK` or `LIWORK` is issued by `XERBLA`. <cyan>Yep, `-1` is the magic "just tell me the size" number.</cyan>
- **`INFO`** — <pink>INTEGER</pink>.
- `= 0`: successful exit <rainbow>slay</rainbow>
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value <red>uh oh, bad arg</red>
- `> 0`: Internal error <red>yikes, internal oopsie</red>
## ✨ Author ✨
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
## ✨ Contributors: ✨
- Inderjit Dhillon, IBM Almaden, USA
- Osni Marques, LBNL/NERSC, USA
- Ken Stanley, Computer Science Division, University of California at Berkeley, USA
- Jason Riedy, Computer Science Division, University of California at Berkeley, USA
<rainbow>And now, a whole new subroutine enters the chat:</rainbow> <hotpink>`sstevx`</hotpink>
subroutine sstevx (character JOBZ, character RANGE, integer N, real, dimension( * ) D, real, dimension( * ) E, real VL, real VU, integer IL, integer IU, real ABSTOL, integer M, real, dimension( * ) W, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer, dimension( * ) IFAIL, integer INFO)
**`SSTEVX`** computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices.
## ✨ Purpose ✨
`SSTEVX` computes selected eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix `A`. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
## ✨ Parameters ✨
- **`JOBZ`** — `JOBZ` is `CHARACTER*1`
- `= 'N'`: Compute eigenvalues only;
- `= 'V'`: Compute eigenvalues and eigenvectors.
- **`RANGE`** — `RANGE` is `CHARACTER*1`
- `= 'A'`: all eigenvalues will be found.
- `= 'V'`: all eigenvalues in the half-open interval `(VL,VU]` will be found.
- `= 'I'`: the `IL`-th through `IU`-th eigenvalues will be found.
- **`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 `A`. On exit, `D` may be multiplied by a constant factor chosen to avoid over/underflow in computing the eigenvalues. <cyan>So it can get scaled, no cap.</cyan>
- **`E`** — `E` is `REAL array`, dimension `(max(1,N-1))`. On entry, the (n-1) subdiagonal elements of the tridiagonal matrix `A` in elements `1` to `N-1` of `E`. On exit, `E` may be multiplied by a constant factor chosen to avoid over/underflow in computing the eigenvalues.
- **`VL`** — `VL` is `REAL`. If `RANGE='V'`, the lower bound of the interval to be searched for eigenvalues. `VL < VU`. Not referenced if `RANGE = 'A'` or `'I'`.
- **`VU`** — `VU` is `REAL`. If `RANGE='V'`, the upper bound of the interval to be searched for eigenvalues. `VL < VU`. Not referenced if `RANGE = 'A'` or `'I'`.
- **`IL`** — `IL` is `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`** — `IU` is `INTEGER`. If `RANGE='I'`, the index of the largest eigenvalue to be returned. `1 <= IL <= IU <= N`, if `N > 0`; `IL = 1` and `IU = 0` if `N = 0`. Not referenced if `RANGE = 'A'` or `'V'`.
- **`ABSTOL`** — `ABSTOL` is `REAL`. The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval `[a,b]` of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where `EPS` is the machine precision. If `ABSTOL` is less than or equal to zero, then `EPS*|T|` will be used in its place, where `|T|` is the 1-norm of the tridiagonal matrix.
Eigenvalues will be computed most accurately when `ABSTOL` is set to twice the underflow threshold `2*SLAMCH('S')`, not zero. If this routine returns with `INFO>0`, indicating that some eigenvectors did not converge, try setting `ABSTOL` to `2*SLAMCH('S')`.
See 'Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy,' by Demmel and Kahan, LAPACK Working Note #3.
- **`M`** — `M` is `INTEGER`. The total number of eigenvalues found. `0 <= M <= N`. If `RANGE = 'A'`, `M = N`, and if `RANGE = 'I'`, `M = IU-IL+1`.
- **`W`** — `W` is `REAL array`, dimension `(N)`. The first `M` elements contain the selected eigenvalues in ascending order.
- **`Z`** — `Z` is `REAL array`, dimension `(LDZ, max(1,M) )`. If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` columns of `Z` contain the orthonormal eigenvectors of the matrix `A` corresponding to the selected eigenvalues, with the i-th column of `Z` holding the eigenvector associated with `W(i)`. If an eigenvector fails to converge (`INFO > 0`), then that column of `Z` contains the latest approximation to the eigenvector, and the index of the eigenvector is returned in `IFAIL`. If `JOBZ = 'N'`, then `Z` is not referenced. Note: the user must ensure that at least `max(1,M)` columns are supplied in the array `Z`; if `RANGE = 'V'`, the exact value of `M` is not known in advance and an upper bound must be used.
- **`LDZ`** — `LDZ` is `INTEGER`. The leading dimension of the array `Z`. `LDZ >= 1`, and if `JOBZ = 'V'`, `LDZ >= max(1,N)`.
- **`WORK`** — `WORK` is `REAL array`, dimension `(5*N)`.
- **`IWORK`** — `IWORK` is `INTEGER array`, dimension `(5*N)`.
- **`IFAIL`** — `IFAIL` is `INTEGER array`, dimension `(N)`. If `JOBZ = 'V'`, then if `INFO = 0`, the first `M` elements of `IFAIL` are zero. If `INFO > 0`, then `IFAIL` contains the indices of the eigenvectors that failed to converge. If `JOBZ = 'N'`, then `IFAIL` is not referenced.
- **`INFO`** — `INFO` is `INTEGER`
- `= 0`: successful exit
- `< 0`: if `INFO = -i`, the i-th argument had an illegal value
- `> 0`: if `INFO = i`, then i eigenvectors failed to converge. Their indices are stored in array `IFAIL`.
## ✨ 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.
---
<rainbow>And that's a wrap, bestie!</rainbow> This was LAPACK version 3.11.0, dated Sun Nov 27 2022, for `realOTHEReigen(3)`. <pink>Stay slaying those eigenvalues. 💖✨</pink>
📦 where this page comes from
- debian bookworm: liblapack-doc 3.11.0-2 · BSD-3-Clause
- 💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on : a language model's rewrite of the original, under its terms