realOTHERauxiliary(3)
โจ realOTHERauxiliary(3) โจ
okay so bestie, welcome to the <strong>LAPACK</strong> cinematic universe ๐ </purgle>
๐ NAME ๐
realOTHERauxiliary โ real ๐ฆ
that's it. that's the whole name section. we are SO back, and this is part 1 of 21, so buckle up babe, the lore is only just starting ๐
โจ SYNOPSIS โจ
๐ Functions ๐
integer function ilaslc (M, N, A, LDA)โ ILASLC scans a matrix for its last non-zero column. ๐integer function ilaslr (M, N, A, LDA)โ ILASLR scans a matrix for its last non-zero row. ๐ฆsubroutine slabrd (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)โ SLABRD reduces the firstnbrows and columns of a general matrix to a bidiagonal form. She's giving factorization energy. โจsubroutine slacn2 (N, V, X, ISGN, EST, KASE, ISAVE)โ SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. Lowkey math genius. ๐งฎsubroutine slacon (N, V, X, ISGN, EST, KASE)โ SLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. Same vibe, different name. ๐ซsubroutine sladiv (A, B, C, D, P, Q)โ SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. She's careful with her numbers, periodt. ๐subroutine sladiv1 (A, B, C, D, P, Q)โ no description given, just the signature. She's mysterious. ๐real function sladiv2 (A, B, C, D, R, T)โ no description given either. Iconic silence. ๐subroutine slaein (RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO)โ SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. She's finding eigenvectors like it's nothing. ๐subroutine slaexc (WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, INFO)โ SLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation. Swapping blocks, no big deal. ๐ฅsubroutine slag2 (A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)โ SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow. 2-by-2? She ate. ๐subroutine slags2 (UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)โ SLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such that the rows of the transformed A and B are parallel. Parallel rows, slay. ๐subroutine slagtm (TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)โ SLAGTM performs a matrix-matrix product of the form C = ฮฑAB+ฮฒC, where A is a tridiagonal matrix, B and C are rectangular matrices, and ฮฑ and ฮฒ are scalars, which may be 0, 1, or -1. That's a whole mouthful but she did that. โจsubroutine slagv2 (A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)โ SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. Generalized Schur? She's sophisticated. ๐subroutine slahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)โ SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. QR algorithm but make it fashion. ๐subroutine slahr2 (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)โ SLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. That's a lot of words but she's precise. ๐ฆsubroutine slaic1 (JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)โ SLAIC1 applies one step of incremental condition estimation. One step at a time, bestie. ๐ธsubroutine slaln2 (LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)โ SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. Tiny systems, big brain. ๐งreal function slangt (NORM, N, DL, D, DU)โ SLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. Norm queen. ๐real function slanhs (NORM, N, A, LDA, WORK)โ SLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix. Another norm, another slay. ๐ซreal function slansb (NORM, UPLO, N, K, AB, LDAB, WORK)โ SLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix. Band matrix? She's in her band era. ๐ธreal function slansp (NORM, UPLO, N, AP, WORK)โ SLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form. Packed form, packed with talent. โจreal function slantb (NORM, UPLO, DIAG, N, K, AB, LDAB, WORK)โ SLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular band matrix. Triangular band, iconic. ๐real function slantp (NORM, UPLO, DIAG, N, AP, WORK)โ SLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. Packed triangular, she's flexible. ๐ฆreal function slantr (NORM, UPLO, DIAG, M, N, A, LDA, WORK)โ SLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix. Trapezoidal or triangular, she does both. ๐subroutine slanv2 (A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)โ SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form. 2-by-2 nonsymmetric? She's not symmetric, she's that girl. ๐subroutine slapll (N, X, INCX, Y, INCY, SSMIN)โ SLAPLL measures the linear dependence of two vectors. Measuring the vibes between two vectors. ๐subroutine slapmr (FORWRD, M, N, X, LDX, K)โ SLAPMR rearranges rows of a matrix as specified by a permutation vector. Rearranging rows, rearranging lives. ๐ซsubroutine slapmt (FORWRD, M, N, X, LDX, K)โ SLAPMT performs a forward or backward permutation of the columns of a matrix. Columns getting the permutation treatment. ๐subroutine slaqp2 (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)โ SLAQP2 computes a QR factorization with column pivoting of the matrix block. QR with pivoting, she's pivoting. ๐subroutine slaqps (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)โ SLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. BLAS level 3? She's using the good stuff. ๐ฅsubroutine slaqr0 (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)โ SLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. Eigenvalues and Schur, optional but iconic. โจsubroutine slaqr1 (N, H, LDH, SR1, SI1, SR2, SI2, V)โ SLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. Setting scalar multiples, periodt. ๐subroutine slaqr2 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)โ SLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). Aggressive early deflation? She's aggressive and early. ๐subroutine slaqr3 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)โ SLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). Same energy, different number. ๐ฆsubroutine slaqr4 (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)โ SLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. She's back with eigenvalues. ๐subroutine slaqr5 (WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)โ SLAQR5 performs a single small-bulge multi-shift QR sweep. Small-bulge multi-shift, say that five times fast. ๐ธsubroutine slaqsb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)โ SLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. Scaling with factors fromspbequ. ๐subroutine slaqsp (UPLO, N, AP, S, SCOND, AMAX, EQUED)โ SLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. Packed storage, scaling factors fromsppequ. โจsubroutine slaqtr (LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)โ SLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of special form, in real arithmetic. Real arithmetic, complex special form, she's versatile. ๐subroutine slar1v (N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)โ SLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - ฮปI. That's a whole sentence but she's precise. ๐ฆsubroutine slar2v (N, X, Y, Z, INCX, C, S, INCC)โ SLAR2V applies a vector of plane rotations with real cosines and real sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. Plane rotations, cosines, sines, both sides. She's doing the most. ๐ซsubroutine slarf (SIDE, M, N, V, INCV, TAU, C, LDC, WORK)โ SLARF applies an elementary reflector to a general rectangular matrix. Elementary reflector, but make it fashion. ๐subroutine slarfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)โ SLARFB applies a block reflector or its transpose to a general rectangular matrix. Block reflector or transpose, she's flexible. ๐subroutine slarfg (N, ALPHA, X, INCX, TAU)โ SLARFG generates an elementary reflector (Householder matrix). Householder matrix, iconic. ๐subroutine slarfgp (N, ALPHA, X, INCX, TAU)โ SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. Non-negative beta, she's positive. โจsubroutine slarft (DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)โ SLARFT forms the triangular factor T of a block reflector H = I - vtvH. Triangular factor T, she's factoring. ๐subroutine slarfx (SIDE, M, N, V, TAU, C, LDC, WORK)โ SLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order โค 10. Loop unrolling when order โค 10, she's optimized. ๐ฅsubroutine slarfy (UPLO, N, V, INCV, TAU, C, LDC, WORK)โ SLARFY (that's the whole description, bestie). ๐subroutine slargv (N, X, INCX, Y, INCY, C, INCC)โ SLARGV generates a vector of plane rotations with real cosines and real sines. Generating plane rotations, she's generating. ๐ซsubroutine slarrv (N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, WORK, IWORK, INFO)โ SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. Eigenvectors of tridiagonal, she's computing. ๐subroutine slartv (N, X, INCX, Y, INCY, C, S, INCC)โ SLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair of vectors. Applying rotations to a pair of vectors, she's pairing. ๐subroutine slaswp (N, A, LDA, K1, K2, IPIV, INCX)โ SLASWP performs a series of row interchanges on a general rectangular matrix. Row interchanges, she's swapping. ๐subroutine slatbs (UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)โ SLATBS solves a triangular banded system of equations. Triangular banded system, she's solving. ๐งฎsubroutine slatdf (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)โ SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. LU factorization fromsgetc2, reciprocal Dif-estimate, she's estimating. โจsubroutine slatps (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO)โ SLATPS solves a triangular system of equations with the matrix held in packed storage. Packed storage, triangular system, she's solving again. ๐subroutine slatrs (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)โ no description given, just the signature. She's mysterious. ๐- SLATRS solves a triangular system of equations with the scale factor set to prevent overflow. ๐
subroutine slauu2 (UPLO, N, A, LDA, INFO)โ SLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm). ๐ฆsubroutine slauum (UPLO, N, A, LDA, INFO)โ SLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). ๐ฆsubroutine srscl (N, SA, SX, INCX)โ SRSCL multiplies a vector by the reciprocal of a real scalar. ๐ธsubroutine stprfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK)โ STPRFB applies a real 'triangular-pentagonal' block reflector to a real matrix, which is composed of two blocks. ๐
โจ Detailed Description โจ
This is the group of real other auxiliary routines ๐
๐ Function Documentation ๐
integer function ilaslc (integer M, integer N, real, dimension( lda, * ) A, integer LDA)
ILASLC scans a matrix for its last non-zero column. ๐
Purpose
ILASLC scans A for its last non-zero column. โจ
Parameters
MโMis INTEGER. The number of rows of the matrix A. ๐ธNโNis INTEGER. The number of columns of the matrix A. ๐ธAโAis REAL array, dimension (LDA,N). The m by n matrix A. ๐ฆLDAโLDAis INTEGER. The leading dimension of the array A.LDA >= max(1,M). ๐
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
integer function ilaslr (integer M, integer N, real, dimension( lda, * ) A, integer LDA)
ILASLR scans a matrix for its last non-zero row. ๐
Purpose
ILASLR scans A for its last non-zero row. โจ
Parameters
MโMis INTEGER. The number of rows of the matrix A. ๐ธNโNis INTEGER. The number of columns of the matrix A. ๐ธAโAis REAL array, dimension (LDA,N). The m by n matrix A. ๐ฆLDAโLDAis INTEGER. The leading dimension of the array A.LDA >= max(1,M). ๐
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slabrd (integer M, integer N, integer NB, real, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) TAUQ, real, dimension( * ) TAUP, real, dimension( ldx, * ) X, integer LDX, real, dimension( ldy, * ) Y, integer LDY)
SLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form. ๐ฅ
Purpose
SLABRD reduces the first NB rows and columns of a real general m by n matrix A to upper or lower bidiagonal form by an orthogonal transformation Q**T * A * P, and returns the matrices X and Y which are needed to apply the transformation to the unreduced part of A. ๐
If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower bidiagonal form. ๐
This is an auxiliary routine called by SGEBRD ๐
Parameters
MโMis INTEGER. The number of rows in the matrix A. ๐ธNโNis INTEGER. The number of columns in the matrix A. ๐ธNBโNBis INTEGER. The number of leading rows and columns of A to be reduced. ๐ฆAโAis REAL array, dimension (LDA,N). On entry, the m by n general matrix to be reduced. On exit, the first NB rows and columns of the matrix are overwritten; the rest of the array is unchanged. If m >= n, elements on and below the diagonal in the first NB columns, with the array TAUQ, represent the orthogonal matrix Q as a product of elementary reflectors; and elements above the diagonal in the first NB rows, with the array TAUP, represent the orthogonal matrix P as a product of elementary reflectors. If m < n, elements below the diagonal in the first NB columns, with the array TAUQ, represent the orthogonal matrix Q as a product of elementary reflectors, and elements on and above the diagonal in the first NB rows, with the array TAUP, represent the orthogonal matrix P as a product of elementary reflectors. See Further Details. ๐LDAโLDAis INTEGER. The leading dimension of the array A.LDA >= max(1,M). ๐DโDis REAL array, dimension (NB). The diagonal elements of the first NB rows and columns of the reduced matrix.D(i) = A(i,i). ๐EโEis REAL array, dimension (NB). The off-diagonal elements of the first NB rows and columns of the reduced matrix. ๐TAUQโTAUQis REAL array, dimension (NB). The scalar factors of the elementary reflectors which represent the orthogonal matrix Q. See Further Details. ๐ซTAUPโTAUPis REAL array, dimension (NB). The scalar factors of the elementary reflectors which represent the orthogonal matrix P. See Further Details. ๐ซXโXis REAL array, dimension (LDX,NB). The m-by-nb matrix X required to update the unreduced part of A. ๐ฆLDXโLDXis INTEGER. The leading dimension of the array X.LDX >= max(1,M). ๐YโYis REAL array, dimension (LDY,NB). The n-by-nb matrix Y required to update the unreduced part of A. ๐ฆLDYโLDYis INTEGER. The leading dimension of the array Y.LDY >= max(1,N). ๐
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Further Details ๐
The matrices Q and P are represented as products of elementary reflectors:
Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb)
Each H(i) and G(i) has the form:
H(i) = I - tauq * v * v**T and G(i) = I - taup * u * u**T
where tauq and taup are real scalars, and v and u are real vectors.
If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
The elements of the vectors v and u together form the m-by-nb matrix V and the nb-by-n matrix U**T which are needed, with X and Y, to apply the transformation to the unreduced part of the matrix, using a block update of the form: A := A - V*Y**T - X*U**T.
The contents of A on exit are illustrated by the following examples with nb = 2:
m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 )
( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 )
( v1 v2 a a a ) ( v1 1 a a a a )
( v1 v2 a a a ) ( v1 v2 a a a a )
( v1 v2 a a a ) ( v1 v2 a a a a )
( v1 v2 a a a )
where a denotes an element of the original matrix which is unchanged, vi denotes an element of the vector defining H(i), and ui an element of the vector defining G(i). โจ
subroutine slacn2 (integer N, real, dimension( * ) V, real, dimension( * ) X, integer, dimension( * ) ISGN, real EST, integer KASE, integer, dimension( 3 ) ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. ๐
Purpose
SLACN2 estimates the 1-norm of a square, real matrix A. Reverse communication is used for evaluating matrix-vector products. โจ
Parameters
NโNis INTEGER. The order of the matrix.N >= 1. ๐ธVโVis REAL array, dimension (N). On the final return,V = A*W, whereEST = norm(V)/norm(W)(W is not returned). ๐ฆXโXis REAL array, dimension (N). On an intermediate return, X should be overwritten by
A * X, if KASE=1,
A**T * X, if KASE=2,
and SLACN2 must be re-called with all the other parameters unchanged. ๐
ISGNโISGNis INTEGER array, dimension (N). ๐ESTโESTis REAL. On entry withKASE = 1or2andISAVE(1) = 3, EST should be unchanged from the previous call to SLACN2. On exit, EST is an estimate (a lower bound) fornorm(A). ๐KASEโKASEis INTEGER. On the initial call to SLACN2, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether X should be overwritten byA * XorA**T * X. On the final return from SLACN2, KASE will again be 0. ๐ฅISAVEโISAVEis INTEGER array, dimension (3). ISAVE is used to save variables between calls to SLACN2. ๐
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Further Details ๐
Originally named SONEST, dated March 16, 1988.
This is a thread safe version of SLACON, which uses the array ISAVE in place of a SAVE statement, as follows:
SLACON SLACN2
JUMP ISAVE(1)
J ISAVE(2)
ITER ISAVE(3)
Contributors
Nick Higham, University of Manchester
References
N.J. Higham, 'FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. ๐
slacon
subroutine slacon (integer N, real, dimension( * ) V, real, dimension( * ) X, integer, dimension( * ) ISGN, real EST, integer KASE)
SLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. bestie, it's giving matrix norm energy ๐
Purpose:
SLACON estimates the 1-norm of a square, real matrix A. Reverse communication is used for evaluating matrix-vector products.
Parameters
NโNisINTEGER. The order of the matrix.N >= 1.VโVisREALarray, dimension(N). On the final return,V = A*W, whereEST = norm(V)/norm(W)(Wis not returned). W is a ghost, she's not coming back ๐ปXโXisREALarray, dimension(N). On an intermediate return,Xshould be overwritten by
A * X, if KASE=1,
A**T * X, if KASE=2,
and SLACON must be re-called with all the other parameters unchanged.
ISGNโISGNisINTEGERarray, dimension(N).ESTโESTisREAL. On entry withKASE = 1or2andJUMP = 3,ESTshould be unchanged from the previous call toSLACON. On exit,ESTis an estimate (a lower bound) fornorm(A).KASEโKASEisINTEGER. On the initial call toSLACON,KASEshould be0. On an intermediate return,KASEwill be1or2, indicating whetherXshould be overwritten byA * XorA**T * X. On the final return fromSLACON,KASEwill again be0.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Nick Higham, University of Manchester. Originally named SONEST, dated March 16, 1988.
References:
N.J. Higham, 'FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
sladiv
subroutine sladiv (real A, real B, real C, real D, real P, real Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. no overflow, no problem โจ
Purpose:
SLADIV performs complex division in real arithmetic
a + i*b
p + i*q = ---------
c + i*d
The algorithm is due to Michael Baudin and Robert L. Smith and can be found in the paper 'A Robust Complex Division in Scilab'
Parameters
AโAisREALBโBisREALCโCisREALDโDisREAL. The scalars a, b, c, and d in the above expression.PโPisREALQโQisREAL. The scalars p and q in the above expression.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
slaein
subroutine slaein (logical RIGHTV, logical NOINIT, integer N, real, dimension( ldh, * ) H, integer LDH, real WR, real WI, real, dimension( * ) VR, real, dimension( * ) VI, real, dimension( ldb, * ) B, integer LDB, real, dimension( * ) WORK, real EPS3, real SMLNUM, real BIGNUM, integer INFO)
SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. eigenvector era ๐
Purpose:
SLAEIN uses inverse iteration to find a right or left eigenvector corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg matrix H.
Parameters
RIGHTVโRIGHTVisLOGICAL= .TRUE.: compute right eigenvector;= .FALSE.: compute left eigenvector.
NOINITโNOINITisLOGICAL= .TRUE.: no initial vector supplied in(VR,VI).= .FALSE.: initial vector supplied in(VR,VI).
NโNisINTEGER. The order of the matrixH.N >= 0.HโHisREALarray, dimension(LDH,N). The upper Hessenberg matrixH.LDHโLDHisINTEGER. The leading dimension of the arrayH.LDH >= max(1,N).WRโWRisREALWIโWIisREAL. The real and imaginary parts of the eigenvalue ofHwhose corresponding right or left eigenvector is to be computed.VRโVRisREALarray, dimension(N)VIโVIisREALarray, dimension(N). On entry, ifNOINIT = .FALSE.andWI = 0.0,VRmust contain a real starting vector for inverse iteration using the real eigenvalueWR; ifNOINIT = .FALSE.andWI.ne.0.0,VRandVImust contain the real and imaginary parts of a complex starting vector for inverse iteration using the complex eigenvalue(WR,WI); otherwiseVRandVIneed not be set. On exit, ifWI = 0.0(real eigenvalue),VRcontains the computed real eigenvector; ifWI.ne.0.0(complex eigenvalue),VRandVIcontain the real and imaginary parts of the computed complex eigenvector. The eigenvector is normalized so that the component of largest magnitude has magnitude1; here the magnitude of a complex number(x,y)is taken to be|x| + |y|.VIis not referenced ifWI = 0.0.BโBisREALarray, dimension(LDB,N)LDBโLDBisINTEGER. The leading dimension of the arrayB.LDB >= N+1.WORKโWORKisREALarray, dimension(N)EPS3โEPS3isREAL. A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots.SMLNUMโSMLNUMisREAL. A machine-dependent value close to the underflow threshold.BIGNUMโBIGNUMisREAL. A machine-dependent value close to the overflow threshold.INFOโINFOisINTEGER= 0: successful exit= 1: inverse iteration did not converge;VRis set to the last iterate, and so isVIifWI.ne.0.0.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
slaexc
subroutine slaexc (logical WANTQ, integer N, real, dimension( ldt, * ) T, integer LDT, real, dimension( ldq, * ) Q, integer LDQ, integer J1, integer N1, integer N2, real, dimension( * ) WORK, integer INFO)
SLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation. swap swap swap ๐
Purpose:
SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in an upper quasi-triangular matrix T by an orthogonal similarity transformation.
T must be in Schur canonical form, that is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has its diagonal elements equal and its off-diagonal elements of opposite sign.
Parameters
WANTQโWANTQisLOGICAL= .TRUE.: accumulate the transformation in the matrixQ;= .FALSE.: do not accumulate the transformation.
NโNisINTEGER. The order of the matrixT.N >= 0.TโTisREALarray, dimension(LDT,N). On entry, the upper quasi-triangular matrixT, in Schur canonical form. On exit, the updated matrixT, again in Schur canonical form.LDTโLDTisINTEGER. The leading dimension of the arrayT.LDT >= max(1,N).QโQisREALarray, dimension(LDQ,N). On entry, ifWANTQis.TRUE., the orthogonal matrixQ. On exit, ifWANTQis.TRUE., the updated matrixQ. IfWANTQis.FALSE.,Qis not referenced.LDQโLDQisINTEGER. The leading dimension of the arrayQ.LDQ >= 1; and ifWANTQis.TRUE.,LDQ >= N.J1โJ1isINTEGER. The index of the first row of the first blockT11.N1โN1isINTEGER. The order of the first blockT11.N1 = 0,1or2.N2โN2isINTEGER. The order of the second blockT22.N2 = 0,1or2.WORKโWORKisREALarray, dimension(N)INFOโINFOisINTEGER= 0: successful exit= 1: the transformed matrixTwould be too far from Schur form; the blocks are not swapped andTandQare unchanged.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
slag2
subroutine slag2 (real, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, real SAFMIN, real SCALE1, real SCALE2, real WR1, real WR2, real WI)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow. scaling queen ๐
Purpose:
SLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue problem A - w B, with scaling as necessary to avoid over-/underflow.
The scaling factor 's' results in a modified eigenvalue equation
s A - w B
where s is a non-negative scaling factor chosen so that w, w B, and s A do not overflow and, if possible, do not underflow, either.
Parameters
AโAisREALarray, dimension(LDA, 2). On entry, the 2 x 2 matrixA. It is assumed that its 1-norm is less than1/SAFMIN. Entries less thansqrt(SAFMIN)*norm(A)are subject to being treated as zero.LDAโLDAisINTEGER. The leading dimension of the arrayA.LDA >= 2.BโBโBis a REAL array, dimension(LDB, 2). On entry, it's the 2 x 2 upper triangular matrix B. The vibe is that the one-norm of B is assumed less than1/SAFMIN. The diagonals should be at leastsqrt(SAFMIN)times the largest element of B (in absolute value); if a diagonal is smaller than that, then+/- sqrt(SAFMIN)will be used instead of that diagonal. No cap. ๐LDBโLDBis INTEGER. The leading dimension of the array B.LDB >= 2.SAFMINโSAFMINis REAL. The smallest positive number s.t.1/SAFMINdoes not overflow. (This should always beSLAMCH('S')-- it is an argument in order to avoid having to callSLAMCHfrequently.) ๐ซถSCALE1โSCALE1is REAL. A scaling factor used to avoid over-/underflow in the eigenvalue equation which defines the first eigenvalue. If the eigenvalues are complex, then the eigenvalues are( WR1 +/- WI i ) / SCALE1(which may lie outside the exponent range of the machine),SCALE1=SCALE2, andSCALE1will always be positive. If the eigenvalues are real, then the first (real) eigenvalue isWR1 / SCALE1, but this may overflow or underflow, and in fact,SCALE1may be zero or less than the underflow threshold if the exact eigenvalue is sufficiently large. ๐ญSCALE2โSCALE2is REAL. A scaling factor used to avoid over-/underflow in the eigenvalue equation which defines the second eigenvalue. If the eigenvalues are complex, thenSCALE2=SCALE1. If the eigenvalues are real, then the second (real) eigenvalue isWR2 / SCALE2, but this may overflow or underflow, and in fact,SCALE2may be zero or less than the underflow threshold if the exact eigenvalue is sufficiently large. ๐WR1โWR1is REAL. If the eigenvalue is real, thenWR1isSCALE1times the eigenvalue closest to the(2,2)element ofA B**(-1). If the eigenvalue is complex, thenWR1=WR2isSCALE1times the real part of the eigenvalues. ๐WR2โWR2is REAL. If the eigenvalue is real, thenWR2isSCALE2times the other eigenvalue. If the eigenvalue is complex, thenWR1=WR2isSCALE1times the real part of the eigenvalues. ๐WIโWIis REAL. If the eigenvalue is real, thenWIis zero. If the eigenvalue is complex, thenWIisSCALE1times the imaginary part of the eigenvalues.WIwill always be non-negative. ๐
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slags2 (logical UPPER, real A1, real A2, real A3, real B1, real B2, real B3, real CSU, real SNU, real CSV, real SNV, real CSQ, real SNQ)
SLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such that the rows of the transformed A and B are parallel. ๐ฅ
Purpose โจ
SLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such that if ( UPPER ) then
U**T *A*Q = U**T *( A1 A2 )*Q = ( x 0 )
( 0 A3 ) ( x x )
and
V**T*B*Q = V**T *( B1 B2 )*Q = ( x 0 )
( 0 B3 ) ( x x )
or if ( .NOT.UPPER ) then
U**T *A*Q = U**T *( A1 0 )*Q = ( x x )
( A2 A3 ) ( 0 x )
and
V**T*B*Q = V**T*( B1 0 )*Q = ( x x )
( B2 B3 ) ( 0 x )
The rows of the transformed A and B are parallel, where
U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ )
( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ )
Z**T denotes the transpose of Z. ๐
Parameters ๐
UPPERโUPPERis LOGICAL.= .TRUE.: the input matrices A and B are upper triangular.= .FALSE.: the input matrices A and B are lower triangular.A1โA1is REAL.A2โA2is REAL.A3โA3is REAL. On entry,A1,A2andA3are elements of the input 2-by-2 upper (lower) triangular matrix A.B1โB1is REAL.B2โB2is REAL.B3โB3is REAL. On entry,B1,B2andB3are elements of the input 2-by-2 upper (lower) triangular matrix B.CSUโCSUis REAL.SNUโSNUis REAL. The desired orthogonal matrix U.CSVโCSVis REAL.SNVโSNVis REAL. The desired orthogonal matrix V.CSQโCSQis REAL.SNQโSNQis REAL. The desired orthogonal matrix Q.
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slagtm (character TRANS, integer N, integer NRHS, real ALPHA, real, dimension( * ) DL, real, dimension( * ) D, real, dimension( * ) DU, real, dimension( ldx, * ) X, integer LDX, real BETA, real, dimension( ldb, * ) B, integer LDB)
SLAGTM performs a matrix-matrix product of the form C = ฮฑAB+ฮฒC, where A is a tridiagonal matrix, B and C are rectangular matrices, and ฮฑ and ฮฒ are scalars, which may be 0, 1, or -1. โจ
Purpose ๐
SLAGTM performs a matrix-vector product of the form
B := alpha * A * X + beta * B
where A is a tridiagonal matrix of order N, B and X are N by NRHS matrices, and alpha and beta are real scalars, each of which may be 0., 1., or -1. ๐
Parameters ๐
TRANSโTRANSis CHARACTER*1. Specifies the operation applied to A.= 'N':No transpose,B := alpha * A * X + beta * B= 'T':Transpose,B := alpha * A'* X + beta * B= 'C':Conjugate transpose = Transpose
NโNis INTEGER. The order of the matrix A.N >= 0.NRHSโNRHSis INTEGER. The number of right hand sides, i.e., the number of columns of the matrices X and B.ALPHAโALPHAis REAL. The scalar alpha.ALPHAmust be 0., 1., or -1.; otherwise, it is assumed to be 0.DLโDLis REAL array, dimension(N-1). The (n-1) sub-diagonal elements of T.DโDis REAL array, dimension(N). The diagonal elements of T.DUโDUis REAL array, dimension(N-1). The (n-1) super-diagonal elements of T.XโXis REAL array, dimension(LDX,NRHS). The N by NRHS matrix X.LDXโLDXis INTEGER. The leading dimension of the array X.LDX >= max(N,1).BETAโBETAis REAL. The scalar beta.BETAmust be 0., 1., or -1.; otherwise, it is assumed to be 1.BโBis REAL array, dimension(LDB,NRHS). On entry, the N by NRHS matrix B. On exit, B is overwritten by the matrix expressionB := alpha * A * X + beta * B.LDBโLDBis INTEGER. The leading dimension of the array B.LDB >= max(N,1).
Author ๐ฆ
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slagv2 (real, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, real, dimension( 2 ) ALPHAR, real, dimension( 2 ) ALPHAI, real, dimension( 2 ) BETA, real CSL, real SNL, real CSR, real SNR)
SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. ๐
Purpose ๐
SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. This routine computes orthogonal (rotation) matrices given by CSL, SNL and CSR, SNR such that
- if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0 types), then
[ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ]
[ 0 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ]
[ b11 b12 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ]
[ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ],
- if the pencil (A,B) has a pair of complex conjugate eigenvalues, then
[ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ]
[ a21 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ]
[ b11 0 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ]
[ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ]
where b11 >= b22 > 0.
Parameters ๐ธ
AโAis REAL array, dimension(LDA, 2). On entry, the 2 x 2 matrix A. On exit, A is overwritten by the ``A-part'' of the generalized Schur form.LDAโLDAis INTEGER. The leading dimension of the array A.LDA >= 2.BโBis REAL array, dimension(LDB, 2). On entry, the upper triangular 2 x 2 matrix B. On exit, B is overwritten by the ``B-part'' of the generalized Schur form.LDBโLDBis INTEGER. The leading dimension of the array B.LDB >= 2.ALPHARโALPHARis REAL array, dimension(2).ALPHAIโALPHAIis REAL array, dimension(2).BETAโBETAis REAL array, dimension(2).(ALPHAR(k)+i*ALPHAI(k))/BETA(k)are the eigenvalues of the pencil (A,B),k=1,2,i = sqrt(-1). Note thatBETA(k)may be zero.CSLโCSLis REAL. The cosine of the left rotation matrix.SNLโSNLis REAL. The sine of the left rotation matrix.CSRโCSRis REAL. The cosine of the right rotation matrix.SNRโ
SNR is REAL โ literally the sine of the right rotation matrix ๐ that's the whole vibe, bestie.
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 slahqr โจ
subroutine slahqr (logical WANTT, logical WANTZ, integer N, integer ILO,
integer IHI, real, dimension( ldh, * ) H, integer LDH, real, dimension(
* ) WR, real, dimension( * ) WI, integer ILOZ, integer IHIZ, real,
dimension( ldz, * ) Z, integer LDZ, integer INFO)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. she's the workhorse fr fr ๐ฅ
Purpose: ๐ธ
SLAHQR is an auxiliary routine called by SHSEQR to update the eigenvalues and Schur decomposition already computed by SHSEQR, by dealing with the Hessenberg submatrix in rows and columns ILO to IHI. so she's not the main character, she's the main character's main character, iykyk ๐
Parameters ๐
WANTTโWANTTisLOGICAL= .TRUE.: the full Schur formTis required;= .FALSE.: only eigenvalues are required.
WANTZโWANTZisLOGICAL= .TRUE.: the matrix of Schur vectorsZis required;= .FALSE.: Schur vectors are not required.
NโNisINTEGER. The order of the matrixH.N >= 0.ILOโILOisINTEGERIHIโIHIisINTEGER. It is assumed thatHis already upper quasi-triangular in rows and columnsIHI+1:N, and thatH(ILO,ILO-1) = 0(unlessILO = 1). SLAHQR works primarily with the Hessenberg submatrix in rows and columnsILOtoIHI, but applies transformations to all ofHifWANTTis.TRUE..1 <= ILO <= max(1,IHI);IHI <= N.HโHisREALarray, dimension(LDH,N). On entry, the upper Hessenberg matrixH. On exit, ifINFOis zero and ifWANTTis.TRUE.,His upper quasi-triangular in rows and columnsILO:IHI, with any 2-by-2 diagonal blocks in standard form. IfINFOis zero andWANTTis.FALSE., the contents ofHare unspecified on exit. The output state ofHifINFOis nonzero is given below under the description ofINFO.LDHโLDHisINTEGER. The leading dimension of the arrayH.LDH >= max(1,N).WRโWRisREALarray, dimension(N)WIโWIisREALarray, dimension(N). The real and imaginary parts, respectively, of the computed eigenvaluesILOtoIHIare stored in the corresponding elements ofWRandWI. If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements ofWRandWI, say the i-th and (i+1)th, withWI(i) > 0andWI(i+1) < 0. IfWANTTis.TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned inH, withWR(i) = H(i,i), and, ifH(i:i+1,i:i+1)is a 2-by-2 diagonal block,WI(i) = sqrt(H(i+1,i)*H(i,i+1))andWI(i+1) = -WI(i). ๐ญ the complex conjugate pair thing is so satisfying when it lands rightILOZโILOZisINTEGERIHIZโIHIZisINTEGER. Specify the rows ofZto which transformations must be applied ifWANTZis.TRUE..1 <= ILOZ <= ILO;IHI <= IHIZ <= N.ZโZisREALarray, dimension(LDZ,N). IfWANTZis.TRUE., on entryZmust contain the current matrixZof transformations accumulated by SHSEQR, and on exitZhas been updated; transformations are applied only to the submatrixZ(ILOZ:IHIZ,ILO:IHI). IfWANTZis.FALSE.,Zis not referenced.LDZโLDZisINTEGER. The leading dimension of the arrayZ.LDZ >= max(1,N).INFOโINFOisINTEGER= 0: successful exit> 0: IfINFO = i, SLAHQR failed to compute all the eigenvaluesILOtoIHIin a total of 30 iterations per eigenvalue; elementsi+1:ihiofWRandWIcontain those eigenvalues which have been successfully computed.
If INFO > 0 and WANTT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H.
If INFO > 0 and WANTT is .TRUE., then on exit
(initial value of H)*U = U*(final value of H)
where U is an orthogonal matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI.
If INFO > 0 and WANTZ is .TRUE., then on exit
(final value of Z) = (initial value of Z)*U
where U is the orthogonal matrix in (*) (regardless of the value of WANTT.)
bestie listen ๐ a nonzeroINFOis not a glow-up moment, it means she hit 30 iterations per eigenvalue and tapped out. the eigenvalues she did get are safe inWR/WI, but the rest are out there unconverged. read theWANTT/WANTZcases above very carefully, they tell you exactly whatHandZmean on the way out. no cap.
Author ๐ซถ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details: ๐ฆ
02-96 Based on modifications by
David Day, Sandia National Laboratory, USA
12-04 Further modifications by
Ralph Byers, University of Kansas, USA
This is a modified version of SLAHQR from LAPACK version 3.0.
It is (1) more robust against overflow and underflow and
(2) adopts the more conservative Ahues & Tisseur stopping
criterion (LAWN 122, 1997).
โจ subroutine slahr2 โจ
subroutine slahr2 (integer N, integer K, integer NB, real, dimension( lda, *
) A, integer LDA, real, dimension( nb ) TAU, real, dimension( ldt, nb )
T, integer LDT, real, dimension( ldy, nb ) Y, integer LDY)
SLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. โจ
Purpose: ๐ธ
SLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1) matrix A so that elements below the k-th subdiagonal are zero. The reduction is performed by an orthogonal similarity transformation Q**T * A * Q. The routine returns the matrices V and T which determine Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.
This is an auxiliary routine called by SGEHRD.
Parameters ๐
NโNisINTEGER. The order of the matrixA.KโKisINTEGER. The offset for the reduction. Elements below the k-th subdiagonal in the firstNBcolumns are reduced to zero.K < N.NBโNBisINTEGER. The number of columns to be reduced.AโAisREALarray, dimension(LDA,N-K+1). On entry, the n-by-(n-k+1) general matrixA. On exit, the elements on and above the k-th subdiagonal in the firstNBcolumns are overwritten with the corresponding elements of the reduced matrix; the elements below the k-th subdiagonal, with the arrayTAU, represent the matrixQas a product of elementary reflectors. The other columns ofAare unchanged. See Further Details.LDAโLDAisINTEGER. The leading dimension of the arrayA.LDA >= max(1,N).TAUโTAUisREALarray, dimension(NB). The scalar factors of the elementary reflectors. See Further Details.TโTisREALarray, dimension(LDT,NB). The upper triangular matrixT.LDTโLDTisINTEGER. The leading dimension of the arrayT.LDT >= NB.YโYisREALarray, dimension(LDY,NB). The n-by-nb matrixY.LDYโLDYisINTEGER. The leading dimension of the arrayY.LDY >= N.
Author ๐ซถ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details: ๐ฆ
The matrix Q is represented as a product of nb elementary reflectors
Q = H(1) H(2) . . . H(nb).
Each H(i) has the form
H(i) = I - tau * v * v**T
where tau is a real scalar, and v is a real vector with
v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
A(i+k+1:n,i), and tau in TAU(i).
The elements of the vectors v together form the (n-k+1)-by-nb matrix
V which is needed, with T and Y, to apply the transformation to the
unreduced part of the matrix, using an update of the form:
A := (I - V*T*V**T) * (A - Y*V**T).
The contents of A on exit are illustrated by the following example
with n = 7, k = 3 and nb = 2:
( a a a a a )
( a a a a a )
( a a a a a )
( h h a a a )
( v1 h a a a )
( v1 v2 a a a )
( v1 v2 a a a )
where a denotes an element of the original matrix A, h denotes a
modified element of the upper Hessenberg matrix H, and vi denotes an
element of the vector defining H(i).
This subroutine is a slight modification of LAPACK-3.0's SLAHRD
incorporating improvements proposed by Quintana-Orti and Van de
Gejin. Note that the entries of A(1:K,2:NB) differ from those
returned by the original LAPACK-3.0's SLAHRD routine. (This
subroutine is not backward compatible with LAPACK-3.0's SLAHRD.)
heads up bestie ๐จ that last bit is a real compatibility warning โ the entries ofA(1:K,2:NB)are not what LAPACK-3.0'sSLAHRDhanded back. not backward compatible, periodt. don't assume the old output shape.
References: ๐
Gregorio Quintana-Orti and Robert van de Geijn, 'Improving the performance of reduction to Hessenberg form,' ACM Transactions on Mathematical Software, 32(2):180-194, June 2006.
โจ subroutine slaic1 โจ
subroutine slaic1 (integer JOB, integer J, real, dimension( j ) X, real
SEST, real, dimension( j ) W, real GAMMA, real SESTPR, real S, real C)
SLAIC1 applies one step of incremental condition estimation. ๐
Purpose: ๐ธ
SLAIC1 applies one step of incremental condition estimation in its simplest version:
So like, here's the whole vibe of what SLAIC1 is actually doing โ bestie, it's giving singular vector glow-up energy โจ
Let x, with twonorm(x) = 1, be an approximate singular vector of an j-by-j lower triangular matrix L, such that
twonorm(L*x) = sest
Then SLAIC1 computes sestpr, s, c such that the vector
[ s*x ]
xhat = [ c ]
is an approximate singular vector of
[ L 0 ]
Lhat = [ w**T gamma ]
in the sense that
twonorm(Lhat*xhat) = sestpr.
Depending on JOB, an estimate for the largest or smallest singular value is computed. ๐
Lowkey important: note that[s c]**Tandsestpr**2is an eigenpair of the system
diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
[ gamma ]
where alpha = x**T*w. No cap, that's the whole tea. ๐ต
Parameters ๐
JOBโJOB is INTEGER= 1: an estimate for the largest singular value is computed.= 2: an estimate for the smallest singular value is computed.
JโJ is INTEGER. Length ofXandW.XโX is REAL array, dimension (J). The j-vectorx.SESTโSEST is REAL. Estimated singular value of j by j matrixL.WโW is REAL array, dimension (J). The j-vectorw.GAMMAโGAMMA is REAL. The diagonal elementgamma.SESTPRโSESTPR is REAL. Estimated singular value of (j+1) by (j+1) matrixLhat.SโS is REAL. Sine needed in formingxhat.CโC is REAL. Cosine needed in formingxhat.
Author ๐
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
slaln2 โ the 1-by-1 / 2-by-2 linear system solver ๐ฆ
subroutine slaln2 (logical LTRANS, integer NA, integer NW, real SMIN, real
CA, real, dimension( lda, * ) A, integer LDA, real D1, real D2, real,
dimension( ldb, * ) B, integer LDB, real WR, real WI, real, dimension(
ldx, * ) X, integer LDX, real SCALE, real XNORM, integer INFO)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. Tiny but mighty, fr fr. ๐
Purpose โจ
SLALN2 solves a system of the form (ca A - w D ) X = s B or (ca A**T - w D) X = s B with possible scaling ('s') and perturbation of A. (A**T means A-transpose, bestie.)
A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA real diagonal matrix, w is a real or complex value, and X and B are NA x 1 matrices โ real if w is real, complex if w is complex. NA may be 1 or 2. ๐
If w is complex, X and B are represented as NA x 2 matrices, the first column of each being the real part and the second being the imaginary part. Classic.
's' is a scaling factor (<= 1), computed by SLALN2, which is so chosen that X can be computed without overflow. X is further scaled if necessary to assure that norm(ca A - w D)*norm(X) is less than overflow. ๐ซถ
If both singular values of (ca A - w D) are less than SMIN, SMIN*identity will be used instead of (ca A - w D). If only one singular value is less than SMIN, one element of (ca A - w D) will be perturbed enough to make the smallest singular value roughly SMIN. If both singular values are at least SMIN, (ca A - w D) will not be perturbed. In any case, the perturbation will be at most some small multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values are computed by infinity-norm approximations, and thus will only be correct to a factor of 2 or so. It's giving "approximately iconic." ๐ธ
Note: all input quantities are assumed to be smaller than overflow by a reasonable factor. (See BIGNUM.)
Parameters ๐
LTRANSโLTRANS is LOGICAL=.TRUE.: A-transpose will be used.=.FALSE.: A will be used (not transposed.)
NAโNA is INTEGER. The size of the matrixA. It may (only) be 1 or 2.NWโNW is INTEGER. 1 if'w'is real, 2 if'w'is complex. It may only be 1 or 2.SMINโSMIN is REAL. The desired lower bound on the singular values ofA. This should be a safe distance away from underflow or overflow, say, between(underflow/machine precision)and(machine precision * overflow ). (SeeBIGNUMandULP.)CAโCA is REAL. The coefficientc, whichAis multiplied by.AโA is REAL array, dimension (LDA,NA). The NA x NA matrixA.LDAโLDA is INTEGER. The leading dimension ofA. It must be at leastNA.D1โD1 is REAL. The 1,1 element in the diagonal matrixD.D2โD2 is REAL. The 2,2 element in the diagonal matrixD. Not used ifNA=1.BโB is REAL array, dimension (LDB,NW). The NA x NW matrixB(right-hand side). IfNW=2('w'is complex), column 1 contains the real part ofBand column 2 contains the imaginary part.LDBโLDB is INTEGER. The leading dimension ofB. It must be at leastNA.WRโWR is REAL. The real part of the scalar'w'.WIโWI is REAL. The imaginary part of the scalar'w'. Not used ifNW=1.XโX is REAL array, dimension (LDX,NW). The NA x NW matrixX(unknowns), as computed bySLALN2. IfNW=2('w'is complex), on exit, column 1 will contain the real part ofXand column 2 will contain the imaginary part.LDXโLDX is INTEGER. The leading dimension ofX. It must be at leastNA.SCALEโSCALE is REAL. The scale factor thatBmust be multiplied by to insure that overflow does not occur when computingX. Thus,(ca A - w D) Xwill beSCALE*B, notB(ignoring perturbations ofA.) It will be at most 1.XNORMโXNORM is REAL. The infinity-norm ofX, whenXis regarded as an NA x NW real matrix.INFOโINFO is INTEGER. An error flag. It will be set to zero if no error occurs, a negative number if an argument is in error, or a positive number ifca A - w Dhad to be perturbed. The possible values are:= 0: No error occurred, and(ca A - w D)did not have to be perturbed.= 1:(ca A - w D)had to be perturbed to make its smallest (or only) singular value greater thanSMIN.
NOTE: In the interests of speed, this routine does not check the inputs for errors. ๐ญ (she's speedy, not careful, we respect it)
Author ๐
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
slangt โ the tridiagonal norm calculator ๐
real function slangt (character NORM, integer N, real, dimension( * ) DL,
real, dimension( * ) D, real, dimension( * ) DU)
SLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. It's giving norm queen. ๐
Purpose โจ
SLANGT returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real tridiagonal matrix A.
Returns ๐
SLANGT = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. Periodt. ๐
Parameters ๐ฆ
NORMโNORM is CHARACTER*1. Specifies the value to be returned inSLANGTas described above.NโN is INTEGER. The order of the matrixA.N >= 0. WhenN = 0,SLANGTis set to zero.DLโDL is REAL array, dimension (N-1). The (n-1) sub-diagonal elements ofA.DโD is REAL array, dimension (N). The diagonal elements ofA.DUโDU is REAL array, dimension (N-1). The (n-1) super-diagonal elements ofA.
Author ๐
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
slanhs โ the upper Hessenberg norm girlie ๐ง
real function slanhs (character NORM, integer N, real, dimension( lda, * )
A, integer LDA, real, dimension( * ) WORK)
SLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix. Obsessed with her tbh. ๐
Purpose โจ
SLANHS returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hessenberg matrix A.
Returns ๐
SLANHS = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. Same energy, different matrix. โจ
Parameters ๐ฆ
NORMโNORM is CHARACTER*1. Specifies the value to be returned inSLANHSas described above.NโN is INTEGER. The order of the matrixA.N >= 0. WhenN = 0,SLANHSis set to zero.AAโAisREAL array, dimension (LDA,N). It's the n by n upper Hessenberg matrix A, and bestie listen โ the part ofAbelow the first sub-diagonal is not referenced, so we literally do not care about it. It's giving "left on read" energy. ๐LDAโLDAisINTEGER. The leading dimension of the arrayA. Constraint check:LDA >= max(N,1). That's the vibe check, no cap. โจWORKโWORKisREAL array, dimension (MAX(1,LWORK)), whereLWORK >= NwhenNORM = 'I'; otherwise,WORKis not referenced. So if you're not asking for the infinity norm,WORKis just... there. Existing. Not contributing. main character who got cut from the episode. ๐ญ
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ฆ slansb ๐ฆ
real function slansb (character NORM, character UPLO, integer N, integer K,
real, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) WORK)
SLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix. Okay so it's the norm-finder but for the banded symmetric girlies. ๐
Purpose:
SLANSB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n symmetric band matrix A, with k super-diagonals. ๐ซถ
Returns
SLANSB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And heads up bestie, max(abs(A(i,j))) is not a consistent matrix norm. It's the chaotic one. It said what it said. ๐
Parameters
NORMโNORMisCHARACTER*1. Specifies the value to be returned inSLANSBas described above. Pick your fighter. โจUPLOโUPLOisCHARACTER*1. Specifies whether the upper or lower triangular part of the band matrixAis supplied.= 'U':Upper triangular part is supplied= 'L':Lower triangular part is supplied
NโNisINTEGER. The order of the matrixA. Constraint:N >= 0. WhenN = 0,SLANSBis set to zero. Zero. Nothing. The void. ๐ธKโKisINTEGER. The number of super-diagonals or sub-diagonals of the band matrixA. Constraint:K >= 0.ABโABisREAL array, dimension (LDAB,N). The upper or lower triangle of the symmetric band matrixA, stored in the firstK+1rows ofAB. The j-th column ofAis stored in the j-th column of the arrayABas follows:- if
UPLO = 'U',AB(k+1+i-j,j) = A(i,j)formax(1,j-k)<=i<=j; - if
UPLO = 'L',AB(1+i-j,j) = A(i,j)forj<=i<=min(n,j+k).
- if
Memorize the index math, I'm so serious. This is the whole ballgame. ๐ฅ
LDABโLDABisINTEGER. The leading dimension of the arrayAB. Constraint:LDAB >= K+1.WORKโWORKisREAL array, dimension (MAX(1,LWORK)), whereLWORK >= NwhenNORM = 'I'or'1'or'O'; otherwise,WORKis not referenced. Note it's a different set of norms that needWORKhere than inslansy, so don't get cute and assume. ๐
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ slansp ๐
real function slansp (character NORM, character UPLO, integer N, real,
dimension( * ) AP, real, dimension( * ) WORK)
SLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form. Packed!! Like, no wasted space, she's efficient, she's thriving. ๐
Purpose:
SLANSP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix A, supplied in packed form. ๐ง
Returns
SLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And again: max(abs(A(i,j))) is not a consistent matrix norm. Do not build a theory of norms on it, I'm begging. ๐ญ
Parameters
NORMโNORMisCHARACTER*1. Specifies the value to be returned inSLANSPas described above.UPLOโUPLOisCHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrixAis supplied.= 'U':Upper triangular part ofAis supplied= 'L':Lower triangular part ofAis supplied
NโNisINTEGER. The order of the matrixA. Constraint:N >= 0. WhenN = 0,SLANSPis set to zero. Nothing to see here, bestie, literally. โจAPโAPisREAL array, dimension (N*(N+1)/2). The upper or lower triangle of the symmetric matrixA, packed columnwise in a linear array. The j-th column ofAis stored in the arrayAPas follows:- if
UPLO = 'U',AP(i + (j-1)*j/2) = A(i,j)for1<=i<=j; - if
UPLO = 'L',AP(i + (j-1)*(2n-j)/2) = A(i,j)forj<=i<=n.
- if
Two different formulas. Two. Depending on UPLO. Write them down. Tattoo them. Whatever works. ๐
WORKโWORKisREAL array, dimension (MAX(1,LWORK)), whereLWORK >= NwhenNORM = 'I'or'1'or'O'; otherwise,WORKis not referenced.
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ slantb ๐
real function slantb (character NORM, character UPLO, character DIAG,
integer N, integer K, real, dimension( ldab, * ) AB, integer LDAB, real,
dimension( * ) WORK)
SLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular band matrix. Triangular AND banded?? She's doing both. Iconic behavior. ๐
Purpose:
SLANTB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n triangular band matrix A, with ( k + 1 ) diagonals. ๐ฆ
Returns
SLANTB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And one more time for the people in the back: max(abs(A(i,j))) is not a consistent matrix norm. ๐ฅ
Parameters
NORMโNORMisCHARACTER*1. Specifies the value to be returned inSLANTBas described above.UPLOโUPLOisCHARACTER*1. Specifies whether the matrixAis upper or lower triangular.= 'U':Upper triangular= 'L':Lower triangular
DIAGโDIAGisCHARACTER*1. Specifies whether or not the matrixAis unit triangular.= 'N':Non-unit triangular= 'U':Unit triangular
NโNisINTEGER. The order of the matrixA. Constraint:N >= 0. WhenN = 0,SLANTBis set to zero. Zero city, population: this matrix. ๐ธKโKisINTEGER. The number of super-diagonals of the matrixAifUPLO = 'U', or the number of sub-diagonals of the matrixAifUPLO = 'L'. Constraint:K >= 0. Same letter, different direction depending onUPLOโ stay alert. ๐ABโABisREAL array, dimension (LDAB,N). The upper or lower triangular band matrixA, stored in the firstk+1rows ofAB. The j-th column ofAis stored in the j-th column of the arrayABas follows:- if
UPLO = 'U',AB(k+1+i-j,j) = A(i,j)formax(1,j-k)<=i<=j; - if
UPLO = 'L',AB(1+i-j,j) = A(i,j)forj<=i<=min(n,j+k).
- if
And here's the plot twist ๐ซข: note that when DIAG = 'U', the elements of the array AB corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one. So those slots? Ignored. Unbothered. Presumed to be 1. Do not read them, do not trust them, do not text them back. ๐
LDABโLDABisINTEGER. The leading dimension of the arrayAB. Constraint:LDAB >= K+1.WORKโWORKisREAL array, dimension (MAX(1,LWORK)), whereLWORK >= NwhenNORM = 'I'; otherwise,WORKis not referenced. Only'I'here. Not'1', not'O'. Just'I'. The way this one is different from the last two?? I need everyone to clock that. ๐
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ slantp ๐
real function slantp (character NORM, character UPLO, character DIAG,
integer N, real, dimension( * ) AP, real, dimension( * ) WORK)
SLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. Triangular. Packed. The crossover event nobody asked for but everybody needed. โญ
Purpose:
SLANTP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix A, supplied in packed form. ๐ง
Returns
SLANTP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). And say it with me one final time: max(abs(A(i,j))) is not a consistent matrix norm. Periodt. ๐
Parameters
NORMโNORMisCHARACTER*1. Specifies the value to be returned inSLANTPas described above.UPLOโUPLOisCHARACTER*1. Specifies whether the matrixAis upper or lower triangular.= 'U':Upper triangular= 'L':Lower triangular
DIAGโDIAGisCHARACTER*1. Specifies whether or not the matrixAis unit triangular.= 'N':Non-unit triangular= 'U':Unit triangular
NNโNisINTEGER. The order of the matrixA.N >= 0. WhenN = 0,SLANTPis set to zero. ๐APโAPisREALarray, dimension(N*(N+1)/2). The upper or lower triangular matrixA, packed columnwise in a linear array. The j-th column ofAis stored in the arrayAPas follows:- if
UPLO = 'U',AP(i + (j-1)*j/2) = A(i,j)for1<=i<=j; - if
UPLO = 'L',AP(i + (j-1)*(2n-j)/2) = A(i,j)forj<=i<=n. Note that whenDIAG = 'U', the elements of the arrayAPcorresponding to the diagonal elements of the matrixAare not referenced, but are assumed to be one. ๐ฆ
- if
WORKโWORKisREALarray, dimension(MAX(1,LWORK)), whereLWORK >= NwhenNORM = 'I'; otherwise,WORKis not referenced. โจ
Author ๐ธ
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
real function slantr ๐
real function slantr (character NORM, character UPLO, character DIAG,
integer M, integer N, real, dimension( lda, * ) A, integer LDA, real,
dimension( * ) WORK)
SLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix. ๐
Purpose โจ
SLANTR returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix A.
Returns ๐
SLANTR
SLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. ๐
Parameters ๐
NORMโNORMisCHARACTER*1. Specifies the value to be returned inSLANTRas described above.UPLOโUPLOisCHARACTER*1. Specifies whether the matrixAis upper or lower trapezoidal.= 'U': Upper trapezoidal= 'L': Lower trapezoidal Note thatAis triangular instead of trapezoidal ifM = N.
DIAGโDIAGisCHARACTER*1. Specifies whether or not the matrixAhas unit diagonal.= 'N': Non-unit diagonal= 'U': Unit diagonal
MโMisINTEGER. The number of rows of the matrixA.M >= 0, and ifUPLO = 'U',M <= N. WhenM = 0,SLANTRis set to zero.NโNisINTEGER. The number of columns of the matrixA.N >= 0, and ifUPLO = 'L',N <= M. WhenN = 0,SLANTRis set to zero.AโAisREALarray, dimension(LDA,N). The trapezoidal matrixA(Ais triangular ifM = N). IfUPLO = 'U', the leading m by n upper trapezoidal part of the arrayAcontains the upper trapezoidal matrix, and the strictly lower triangular part ofAis not referenced. IfUPLO = 'L', the leading m by n lower trapezoidal part of the arrayAcontains the lower trapezoidal matrix, and the strictly upper triangular part ofAis not referenced. Note that whenDIAG = 'U', the diagonal elements ofAare not referenced and are assumed to be one.LDAโLDAisINTEGER. The leading dimension of the arrayA.LDA >= max(M,1).WORKโWORKisREALarray, dimension(MAX(1,LWORK)), whereLWORK >= MwhenNORM = 'I'; otherwise,WORKis not referenced.
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slanv2 ๐ฆ
subroutine slanv2 (real A, real B, real C, real D, real RT1R, real RT1I,
real RT2R, real RT2I, real CS, real SN)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form. ๐
Purpose โจ
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form:
[ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
[ C D ] [ SN CS ] [ CC DD ] [-SN CS ]
where either
CC = 0so thatAAandDDare real eigenvalues of the matrix, orAA = DDandBB*CC < 0, so thatAA + or - sqrt(BB*CC)are complex conjugate eigenvalues.
Parameters ๐ธ
AโAisREAL.BโBisREAL.CโCisREAL.DโDisREAL. On entry, the elements of the input matrix. On exit, they are overwritten by the elements of the standardised Schur form.RT1RโRT1RisREAL.RT1IโRT1IisREAL.RT2RโRT2RisREAL.RT2IโRT2IisREAL. The real and imaginary parts of the eigenvalues. If the eigenvalues are a complex conjugate pair,RT1I > 0.CSโCSisREAL.SNโSNisREAL. Parameters of the rotation matrix.
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Further Details ๐ฅ
Modified by V. Sima, Research Institute for Informatics, Bucharest, Romania, to reduce the risk of cancellation errors, when computing real eigenvalues, and to ensure, if possible, that abs(RT1R) >= abs(RT2R).
subroutine slapll ๐
subroutine slapll (integer N, real, dimension( * ) X, integer INCX, real,
dimension( * ) Y, integer INCY, real SSMIN)
SLAPLL measures the linear dependence of two vectors. ๐
Purpose โจ
Given two column vectors X and Y, let
A = ( X Y ).
The subroutine first computes the QR factorization of A = Q*R, and then computes the SVD of the 2-by-2 upper triangular matrix R. The smaller singular value of R is returned in SSMIN, which is used as the measurement of the linear dependency of the vectors X and Y.
Parameters ๐
NโNisINTEGER. The length of the vectorsXandY.XโXisREALarray, dimension(1+(N-1)*INCX). On entry,Xcontains the N-vectorX. On exit,Xis overwritten.INCXโINCXisINTEGER. The increment between successive elements ofX.INCX > 0.YโYisREALarray, dimension(1+(N-1)*INCY). On entry,Ycontains the N-vectorY. On exit,Yis overwritten.INCYโINCYisINTEGER. The increment between successive elements ofY.INCY > 0.SSMINโSSMINisREAL. The smallest singular value of the N-by-2 matrixA = ( X Y ).
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slapmr ๐ฆ
subroutine slapmr (logical FORWRD, integer M, integer N, real, dimension(
ldx, * ) X, integer LDX, integer, dimension( * ) K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector. ๐
Purpose โจ
SLAPMR rearranges the rows of the M by N matrix X as specified by the permutation K(1),K(2),...,K(M) of the integers 1,...,M. If FORWRD = .TRUE., forward permutation:
X(K(I),*) is moved X(I,*) for I = 1,2,...,M.
If FORWRD = .FALSE., backward permutation:
X(I,*) is moved to X(K(I),*) for I = 1,2,...,M.
Parameters ๐ธ
FORWRDโFORWRDisLOGICAL.= .TRUE., forward permutation= .FALSE., backward permutation
MโMisINTEGER. The number of rows of the matrixX.M >= 0.NโNisINTEGER. The number of columns of the matrixX.N >= 0.XโXisREALarray, dimension(LDX,N). On entry, the M by N matrixX. On exit,Xcontains the permuted matrixX.LDXโLDXisINTEGER. The leading dimension of the arrayX,LDX >= MAX(1,M).KโKisINTEGERarray, dimension(M). On entry,Kcontains the permutation vector.Kis used as internal workspace, but reset to its original value on output.
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slapmt ๐
subroutine slapmt (logical FORWRD, integer M, integer N, real, dimension(
ldx, * ) X, integer LDX, integer, dimension( * ) K)
SLAPMT performs a forward or backward permutation of the columns of a matrix. ๐
Purpose โจ
SLAPMT rearranges the columns of the M by N matrix X as specified by the permutation K(1),K(2),...,K(N) of the integers 1,...,N. If FORWRD = .TRUE., forward permutation:
X(*,K(J)) is moved X(*,J) for J = 1,2,...,N.
If FORWRD = .FALSE., backward permutation:
X(*,J) is moved to X(*,K(J)) for J = 1,2,...,N.
Parameters ๐
FORWRDโFORWRDisLOGICAL.= .TRUE., forward permutation= .FALSE., backward permutation
MโMisINTEGER. The number of rows of the matrixX.M >= 0.NโNisINTEGER. The number of columns of the matrixX.N >= 0.XโXisREALarray, dimension(LDX,N). On entry, the M by N matrixX. On exit,Xcontains the permuted matrixX.LDXโLDXisINTEGER. The leading dimension of the arrayX,LDX >= MAX(1,M).KโKisINTEGERarray, dimension(N). On entry,Kcontains the permutation vector.Kis used as internal workspace, but reset to its original value on output.
Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ฏ Contributors
- G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
- X. Sun, Computer Science Dept., Duke University, USA
- Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.
๐ References
- LAPACK Working Note 176
๐ฆ slaqps โ QR steps but make them BLAS-3, bestie
subroutine slaqps (integer M, integer N, integer OFFSET, integer NB, integer
KB, real, dimension( lda, * ) A, integer LDA, integer, dimension( * )
JPVT, real, dimension( * ) TAU, real, dimension( * ) VN1, real,
dimension( * ) VN2, real, dimension( * ) AUXV, real, dimension( ldf, * )
F, integer LDF)
okay so SLAQPS does a step of QR factorization with column pivoting on a real m-by-n matrix A, and it flexes BLAS level 3 while doing it. it's giving "let the big BLAS routines carry" energy fr fr โจ
๐ฏ Purpose
SLAQPS computes a step of QR factorization with column pivoting of a real M-by-N matrix A by using Blas-3. It tries to factorize NB columns from A starting from the row OFFSET+1, and updates all of the matrix with Blas-3 xGEMM. ๐ค๐ช
ngl though โ in some cases, due to catastrophic cancellations ๐ญ, it cannot factorize NB columns. so the actual number of factorized columns gets handed back in KB. that's the plot twist, that's the drama. ๐
Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. periodt. ๐
๐ธ Parameters
MโMis INTEGER. the number of rows of the matrixA.M >= 0. slay.NโNis INTEGER. the number of columns of the matrixA.N >= 0OFFSETโOFFSETis INTEGER. the number of rows ofAthat have been factorized in previous steps.NBโNBis INTEGER. the number of columns to factorize. the dream, the goal ๐KBโKBis INTEGER. the number of columns actually factorized. reality check ๐ซAโAis REAL array, dimension(LDA,N). On entry, the M-by-N matrixA. On exit, blockA(OFFSET+1:M,1:KB)is the triangular factor obtained and blockA(1:OFFSET,1:N)has been accordingly pivoted, but no factorized. the rest of the matrix, blockA(OFFSET+1:M,KB+1:N), has been updated. โจLDAโLDAis INTEGER. the leading dimension of the arrayA.LDA >= max(1,M).JPVTโJPVTis INTEGER array, dimension(N).JPVT(I) = K <==> Column K of the full matrix A has been permuted into position I in AP.yes queen, tell them exactly where everything went ๐TAUโTAUis REAL array, dimension(KB). the scalar factors of the elementary reflectors.VN1โVN1is REAL array, dimension(N). the vector with the partial column norms.VN2โVN2is REAL array, dimension(N). the vector with the exact column norms.AUXVโAUXVis REAL array, dimension(NB). auxiliary vector. it's just there to help, we love a supportive icon ๐ซถFโFis REAL array, dimension(LDF,NB). MatrixF**T = L*Y**T*A.LDFโLDFis INTEGER. the leading dimension of the arrayF.LDF >= max(1,N).
โ๏ธ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ฏ Contributors
- G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
- X. Sun, Computer Science Dept., Duke University, USA
Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.
๐ References
- LAPACK Working Note 176
๐ฅ slaqr0 โ eigenvalues of a Hessenberg matrix + optional Schur tea
subroutine slaqr0 (logical WANTT, logical WANTZ, integer N, integer ILO,
integer IHI, real, dimension( ldh, * ) H, integer LDH, real, dimension(
* ) WR, real, dimension( * ) WI, integer ILOZ, integer IHIZ, real,
dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer
LWORK, integer INFO)
SLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. main character energy only ๐ฅ๐
๐ฏ Purpose
SLAQR0 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur form), and Z is the orthogonal matrix of Schur vectors. iconic, no notes ๐
optionally Z may be postmultiplied into an input orthogonal matrix Q so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. the way that algebra just works?? obsessed. โจ
๐ธ Parameters
WANTTโWANTTis LOGICAL.= .TRUE.: the full Schur formTis required;= .FALSE.: only eigenvalues are required.
WANTZโWANTZis LOGICAL.= .TRUE.: the matrix of Schur vectorsZis required;= .FALSE.: Schur vectors are not required.
NโNis INTEGER. the order of the matrixH.N >= 0.ILOโILOis INTEGERIHIโIHIis INTEGER. it is assumed thatHis already upper triangular in rows and columns1:ILO-1andIHI+1:Nand, ifILO > 1,H(ILO,ILO-1)is zero.ILOandIHIare normally set by a previous call toSGEBAL, and then passed toSGEHRDwhen the matrix output bySGEBALis reduced to Hessenberg form. otherwise,ILOandIHIshould be set to1andN, respectively. ifN > 0, then1 <= ILO <= IHI <= N. ifN = 0, thenILO = 1andIHI = 0.HโHis REAL array, dimension(LDH,N). on entry, the upper Hessenberg matrixH. on exit, ifINFO = 0andWANTTis.TRUE., thenHcontains the upper quasi-triangular matrixTfrom the Schur decomposition (the Schur form); 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, withH(i,i) = H(i+1,i+1)andH(i+1,i)*H(i,i+1) < 0. ifINFO = 0andWANTTis.FALSE., then the contents ofHare unspecified on exit. (the output value ofHwhenINFO > 0is given under the description ofINFObelow.) ๐
this subroutine may explicitly set H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. so it's zeroing stuff out, don't be shocked bestie ๐
LDHโLDHis INTEGER. the leading dimension of the arrayH.LDH >= max(1,N).WRโWRis REAL array, dimension(IHI)WIโWIis REAL array, dimension(IHI). the real and imaginary parts, respectively, of the computed eigenvalues ofH(ILO:IHI,ILO:IHI)are stored inWR(ILO:IHI)andWI(ILO:IHI). if two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements ofWRandWI, say the i-th and (i+1)th, withWI(i) > 0andWI(i+1) < 0. ifWANTTis.TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned inH, withWR(i) = H(i,i)and, ifH(i:i+1,i:i+1)is a 2-by-2 diagonal block,WI(i) = sqrt(-H(i+1,i)*H(i,i+1))andWI(i+1) = -WI(i). conjugate besties staying side by side, we love to see it ๐ซถILOZโILOZis INTEGERIHIZโIHIZis INTEGER. specify the rows ofZto which transformations must be applied ifWANTZis.TRUE..1 <= ILOZ <= ILO; IHI <= IHIZ <= N.ZZโ REAL array, dimension(LDZ,IHI). IfWANTZis.FALSE., bestie, Z is literally not referenced, it's just vibing. IfWANTZis.TRUE., thenZ(ILO:IHI,ILOZ:IHIZ)gets replaced byZ(ILO:IHI,ILOZ:IHIZ)*Uwhere U is the orthogonal Schur factor ofH(ILO:IHI,ILO:IHI). ๐ (And the output value of Z whenINFO > 0? That's covered down in the INFO description below, no skipping ahead.)LDZโ INTEGER. The leading dimension of the array Z. IfWANTZis.TRUE.thenLDZ >= MAX(1,IHIZ). Otherwise,LDZ >= 1. That's it, that's the whole rule. โจWORKโ REAL array, dimensionLWORK. On exit, ifLWORK = -1,WORK(1)returns an estimate of the optimal value forLWORK. Iconic behavior honestly.LWORKโ INTEGER. The dimension of the array WORK.LWORK >= max(1,N)is sufficient, but heads up: LWORK typically as large as6*Nmay be required for optimal performance. A workspace query to determine the optimal workspace size is recommended. ๐ซถ
If LWORK = -1, then SLAQR0 does a workspace query. In this case, SLAQR0 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. (A little peek-without-touching moment, we love a respectful algorithm. ๐)
INFOโ INTEGER. The receipts:= 0: successful exit ๐> 0: ifINFO = i,SLAQR0failed to compute all of the eigenvalues. Elements1:ilo-1andi+1:nof WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.) โ okay so it's not a disaster, it's a partial success situation, ngl.
If INFO > 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. ๐ฆ
If INFO > 0 and WANTT is .TRUE., then on exit
(*) (initial value of H)*U = U*(final value of H)
where U is an orthogonal matrix. The final value of H is upper Hessenberg and quasi-triangular in rows and columns INFO+1 through IHI. ๐ฅ
If INFO > 0 and WANTZ is .TRUE., then on exit
(final value of Z(ILO:IHI,ILOZ:IHIZ)
= (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
where U is the orthogonal matrix in (*) (regardless of the value of WANTT.) ๐
If INFO > 0 and WANTZ is .FALSE., then Z is not accessed. Respectfully left alone, as it should be.
Author ๐ธ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors: ๐
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
References: ๐
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
Performance, SIAM Journal of Matrix Analysis, volume 23, pages
929--947, 2002.
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm
Part II: Aggressive Early Deflation, SIAM Journal of Matrix
Analysis, volume 23, pages 948--973, 2002.
๐ slaqr1 ๐
subroutine slaqr1 (integer N, real, dimension( ldh, * ) H, integer LDH, real
SR1, real SI1, real SR2, real SI2, real, dimension( * ) V)
SLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. Slay.
Purpose: โจ
Given a 2-by-2 or 3-by-3 matrix H, SLAQR1 sets v to a scalar multiple of the first column of the product
(*) K = (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)
scaling to avoid overflows and most underflows. It is assumed that either
1) sr1 = sr2 and si1 = -si2
or
2) si1 = si2 = 0.
This is useful for starting double implicit shift bulges in the QR algorithm. Niche? Yes. Important? Also yes. ๐
Parameters
Nโ INTEGER. Order of the matrix H.Nmust be either 2 or 3. No in-between, bestie, it's giving very specific. ๐ โโ๏ธHโ REAL array, dimension(LDH,N). The 2-by-2 or 3-by-3 matrix H in (*).LDHโ INTEGER. The leading dimension of H as declared in the calling procedure.LDH >= N.SR1โ REAL.SI1โ REAL.SR2โ REAL.SI2โ REAL. The shifts in (*). โจVโ REAL array, dimension(N). A scalar multiple of the first column of the matrix K in (*). ๐
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors: ๐ธ
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
๐ slaqr2 ๐
subroutine slaqr2 (logical WANTT, logical WANTZ, integer N, integer KTOP,
integer KBOT, integer NW, real, dimension( ldh, * ) H, integer LDH,
integer ILOZ, integer IHIZ, real, dimension( ldz, * ) Z, integer LDZ,
integer NS, integer ND, real, dimension( * ) SR, real, dimension( * )
SI, real, dimension( ldv, * ) V, integer LDV, integer NH, real,
dimension( ldt, * ) T, integer LDT, integer NV, real, dimension( ldwv, *
) WV, integer LDWV, real, dimension( * ) WORK, integer LWORK)
SLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). The name is literally the whole plot. ๐ฅ
Purpose: โจ
SLAQR2 is identical to SLAQR3 except that it avoids recursion by calling SLAHQR instead of SLAQR4. (Same energy, different path, iykyk.)
Aggressive early deflation:
This subroutine accepts as input an upper Hessenberg matrix H and performs an orthogonal similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been overwritten by a new Hessenberg matrix that is a perturbation of an orthogonal similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. ๐ฆ (Manifesting those zeros fr fr.)
Parameters
WANTTโ LOGICAL. If.TRUE., then the Hessenberg matrix H is fully updated so that the quasi-triangular Schur factor may be computed (in cooperation with the calling subroutine). If.FALSE., then only enough of H is updated to preserve the eigenvalues. Minimal effort, maximum preservation. ๐WANTZโ LOGICAL. If.TRUE., then the orthogonal matrix Z is updated so so that the orthogonal Schur factor may be computed (in cooperation with the calling subroutine). If.FALSE., then Z is not referenced. (Yes, the "so so" is in the original and yes I'm keeping it, character for character, no cap.)Nโ INTEGER. The order of the matrix H and (ifWANTZis.TRUE.) the order of the orthogonal matrix Z.KTOPโ INTEGER. It is assumed that eitherKTOP = 1orH(KTOP,KTOP-1)=0.KBOTandKTOPtogether determine an isolated block along the diagonal of the Hessenberg matrix. ๐KBOTโ INTEGER. It is assumed without a check that eitherKBOT = NorH(KBOT+1,KBOT)=0.KBOTandKTOPtogether determine an isolated block along the diagonal of the Hessenberg matrix. (Note: without a check โ so don't come crying if you lie to it, bestie. ๐ญ)NWโ INTEGER. Deflation window size.1 <= NW <= (KBOT-KTOP+1).Hโ REAL array, dimension(LDH,N). On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by an orthogonal similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries. โญLDHโ INTEGER. Leading dimension of H just as declared in the calling subroutine.N <= LDH.ILOZโ INTEGER. ๐IHIZโ INTEGER. Specify the rows of Z to which transformations must be applied ifWANTZis.TRUE..1 <= ILOZ <= IHIZ <= N.Zโ REAL array, dimension(LDZ,N). IFWANTZis.TRUE., then on output, the orthogonal similarity transformation mentioned above has been accumulated intoZ(ILOZ:IHIZ,ILOZ:IHIZ)from the right. IfWANTZis.FALSE., then Z is unreferenced. ๐ซถLDZโ INTEGER. The leading dimension of Z just as declared in the calling subroutine.1 <= LDZ.NSโ INTEGER. The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine. โจNDโ INTEGER. The number of converged eigenvalues uncovered by this subroutine. The moment of discovery. ๐SRโ REAL array, dimension(KBOT).SIโ REAL array, dimension(KBOT). On output, the real and imaginary parts of approximate eigenvalues that may be used for shifts are stored inSR(KBOT-ND-NS+1)throughSR(KBOT-ND)andSI(KBOT-ND-NS+1)throughSI(KBOT-ND), respectively. The real and imaginary parts of converged eigenvalues are stored inSR(KBOT-ND+1)throughSR(KBOT)andSI(KBOT-ND+1)throughSI(KBOT), respectively. Get those indices right, bestie, they're doing a whole little choreography. ๐Vโ REAL array, dimension(LDV,NW). An NW-by-NW work array.LDVโ INTEGER. The leading dimension of V just as declared in the calling subroutine.NW <= LDV.NHNHโNHis INTEGER, bestie. It's the number of columns ofT. And the rule isNH >= NW, no exceptions. ๐TโTis REAL array, dimension(LDT,NW). That's it, that's the whole description, she's mysterious. ๐ฆLDTโLDTis INTEGER. It's the leading dimension ofT, just as declared in the calling subroutine. AndNW <= LDT. Periodt.NVโNVis INTEGER, the number of rows of the work arrayWVthat are available for workspace. Constraint check:NV >= NW. ๐ฅWVโWVis REAL array, dimension(LDWV,NW). Just vibing there as your workspace. โจLDWVโLDWVis INTEGER. The leading dimension ofWjust as declared in the calling subroutine. And here's the thing, the original saysNW <= LDVโ yes, it saysLDV, notLDWV. We are not fixing that, we are reporting it, iykyk. ๐WORKโWORKis REAL array, dimension(LWORK). On exit,WORK(1)is set to an estimate of the optimal value ofLWORKfor the given values ofN,NW,KTOPandKBOT. SoWORK(1)is literally the receipts. ๐งพLWORKโLWORKis INTEGER. The dimension of the work arrayWORK.LWORK = 2*NWsuffices, but greater efficiency may result from larger values ofLWORK. Bigger workspace, more slay, that's the tradeoff. ๐
And then the iconic bit: if LWORK = -1, then a workspace query is assumed. SLAQR2 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. So it's a lil peek behind the curtain, not an actual computation. ๐ธ
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Contributors:
- Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA ๐
subroutine slaqr3
subroutine slaqr3 (logical WANTT, logical WANTZ, integer N, integer KTOP,
integer KBOT, integer NW, real, dimension( ldh, * ) H, integer LDH,
integer ILOZ, integer IHIZ, real, dimension( ldz, * ) Z, integer LDZ,
integer NS, integer ND, real, dimension( * ) SR, real, dimension( * )
SI, real, dimension( ldv, * ) V, integer LDV, integer NH, real,
dimension( ldt, * ) T, integer LDT, integer NV, real, dimension( ldwv, *
) WV, integer LDWV, real, dimension( * ) WORK, integer LWORK)
SLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). Basically she's the bouncer at the club who checks which eigenvalues are actually done and escorts them out the door. ๐ชโจ
Purpose:
Aggressive early deflation:
SLAQR3 accepts as input an upper Hessenberg matrix H and performs an orthogonal similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been over-written by a new Hessenberg matrix that is a perturbation of an orthogonal similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries. ๐
"It is to be hoped" is doing so much heavy lifting in that sentence, and honestly? Same. ๐ญ
Parameters
WANTTโWANTTis LOGICAL. If.TRUE., then the Hessenberg matrixHis fully updated so that the quasi-triangular Schur factor may be computed (in cooperation with the calling subroutine). If.FALSE., then only enough ofHis updated to preserve the eigenvalues. So.TRUE.= full glow-up,.FALSE.= the bare minimum, which is real of her. ๐WANTZโWANTZis LOGICAL. If.TRUE., then the orthogonal matrixZis updated so so that the orthogonal Schur factor may be computed (in cooperation with the calling subroutine). If.FALSE., thenZis not referenced. (Yes, the original says "so so". She's excited. We're leaving it. ๐ซถ)NโNis INTEGER. The order of the matrixHand (ifWANTZis.TRUE.) the order of the orthogonal matrixZ.KTOPโKTOPis INTEGER. It is assumed that eitherKTOP = 1orH(KTOP,KTOP-1)=0.KBOTandKTOPtogether determine an isolated block along the diagonal of the Hessenberg matrix.KBOTโKBOTis INTEGER. It is assumed without a check that eitherKBOT = NorH(KBOT+1,KBOT)=0.KBOTandKTOPtogether determine an isolated block along the diagonal of the Hessenberg matrix.
Not me flagging that "without a check" โ girl she will simply trust you on this one. If you lie, that's between you and your segfault. ๐
NWโNWis INTEGER. Deflation window size.1 <= NW <= (KBOT-KTOP+1).HโHis REAL array, dimension(LDH,N). On input the initial N-by-N section ofHstores the Hessenberg matrix undergoing aggressive early deflation. On outputHhas been transformed by an orthogonal similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries.LDHโLDHis INTEGER. Leading dimension ofHjust as declared in the calling subroutine.N <= LDH.ILOZโILOZis INTEGER. (That's the whole description. Iconic minimalism. ๐)IHIZโIHIZis INTEGER. Specify the rows ofZto which transformations must be applied ifWANTZis.TRUE..1 <= ILOZ <= IHIZ <= N.ZโZis REAL array, dimension(LDZ,N). IFWANTZis.TRUE., then on output, the orthogonal similarity transformation mentioned above has been accumulated intoZ(ILOZ:IHIZ,ILOZ:IHIZ)from the right. IfWANTZis.FALSE., thenZis unreferenced.LDZโLDZis INTEGER. The leading dimension ofZjust as declared in the calling subroutine.1 <= LDZ.NSโNSis INTEGER. The number of unconverged (ie approximate) eigenvalues returned inSRandSIthat may be used as shifts by the calling subroutine. The "not quite there yet but we're working on it" girlies. ๐ฑNDโNDis INTEGER. The number of converged eigenvalues uncovered by this subroutine. These ones ate and left no crumbs. ๐ฝ๏ธSRโSRis REAL array, dimension(KBOT).SIโSIis REAL array, dimension(KBOT). On output, the real and imaginary parts of approximate eigenvalues that may be used for shifts are stored inSR(KBOT-ND-NS+1)throughSR(KBOT-ND)andSI(KBOT-ND-NS+1)throughSI(KBOT-ND), respectively. The real and imaginary parts of converged eigenvalues are stored inSR(KBOT-ND+1)throughSR(KBOT)andSI(KBOT-ND+1)throughSI(KBOT), respectively. Two zones, two vibes: shifts up front, converged in the back. ๐VโVis REAL array, dimension(LDV,NW). An NW-by-NW work array.LDVโLDVis INTEGER. The leading dimension ofVjust as declared in the calling subroutine.NW <= LDV.NHโNHis INTEGER. The number of columns ofT.NH >= NW.TโTis REAL array, dimension(LDT,NW).LDTโLDTis INTEGER. The leading dimension ofTjust as declared in the calling subroutine.NW <= LDT.NVโNVis INTEGER. The number of rows of work arrayWVavailable for workspace.NV >= NW.WVโWVis REAL array, dimension(LDWV,NW).LDWVโLDWVis INTEGER. The leading dimension ofWjust as declared in the calling subroutine.NW <= LDV.
Again: it saysLDV, notLDWV. Copying faithfully, no cap. ๐ซก
WORKโWORKis REAL array, dimension(LWORK). On exit,WORK(1)is set to an estimate of the optimal value ofLWORKfor the given values ofN,NW,KTOPandKBOT.LWORKโLWORKis INTEGER. The dimension of the work arrayWORK.LWORK = 2*NWsuffices, but greater efficiency may result from larger values ofLWORK.
If LWORK = -1, then a workspace query is assumed; SLAQR3 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Contributors:
- Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA ๐
subroutine slaqr4
subroutine slaqr4 (logical WANTT, logical WANTZ, integer N, integer ILO,
integer IHI, real, dimension( ldh, * ) H, integer LDH, real, dimension(
* ) WR, real, dimension( * ) WI, integer ILOZ, integer IHIZ, real,
dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer
LWORK, integer INFO)
SLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. The main character of the QR algorithm saga, fr fr. ๐
Purpose:
SLAQR4 implements one level of recursion for SLAQR0. It is a complete implementation of the small bulge multi-shift QR algorithm. It may be called by SLAQR0 and, for large enough deflation window size, it may be called by SLAQR3. This subroutine is identical to SLAQR0 except that it calls SLAQR2 instead of SLAQR3. Twin sisters, different friend groups. ๐ฏ
SLAQR4 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur form), and Z is the orthogonal matrix of Schur vectors.
Optionally Z may be postmultiplied into an input orthogonal matrix Q so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. Slay. โจ
Parameters
WANTTโWANTTis LOGICAL.= .TRUE.: the full Schur formTis required;= .FALSE.: only eigenvalues are required.
WANTZโWANTZis LOGICAL.= .TRUE.: the matrix of Schur vectorsZis required;= .FALSE.: Schur vectors are not required.
NโNis INTEGER. The order of the matrixH.N >= 0.ILOโILOis INTEGER. (Short, sweet, says nothing else. She's a woman of few words. ๐ )IHIโIHIis INTEGER.
so yeah โ IHI is INTEGER, and here's the rest of her whole deal:
It is assumed that H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N and, if ILO > 1, H(ILO,ILO-1) is zero. ILO and IHI are normally set by a previous call to SGEBAL, and then passed to SGEHRD when the matrix output by SGEBAL is reduced to Hessenberg form. Otherwise, ILO and IHI should be set to 1 and N, respectively. If N > 0, then 1 <= ILO <= IHI <= N. If N = 0, then ILO = 1 and IHI = 0.
HโHisREAL array, dimension(LDH,N)๐
On entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper quasi-triangular matrix T from the Schur decomposition (the Schur form); 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is .FALSE., then the contents of H are unspecified on exit. (The output value of H when INFO > 0 is given under the description of INFO below.)
This subroutine may explicitly set H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
LDHโLDHis INTEGER. The leading dimension of the arrayH.LDH >= max(1,N). ๐ฆWRโWRisREAL array, dimension(IHI)๐WIโWIisREAL array, dimension(IHI)
The real and imaginary parts, respectively, of the computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI) and WI(ILO:IHI). If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). โจ
ILOZโILOZis INTEGER ๐IHIZโIHIZis INTEGER
Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
ZโZisREAL array, dimension(LDZ,IHI)
If WANTZ is .FALSE., then Z is not referenced. If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the orthogonal Schur factor of H(ILO:IHI,ILO:IHI). (The output value of Z when INFO > 0 is given under the description of INFO below.)
LDZโLDZis INTEGER. The leading dimension of the arrayZ. ifWANTZis.TRUE.thenLDZ >= MAX(1,IHIZ). Otherwise,LDZ >= 1. ๐ธWORKโWORKisREAL array, dimensionLWORK
On exit, if LWORK = -1, WORK(1) returns an estimate of the optimal value for LWORK.
LWORKโLWORKis INTEGER
The dimension of the array WORK. LWORK >= max(1,N) is sufficient, but LWORK typically as large as 6*N may be required for optimal performance. A workspace query to determine the optimal workspace size is recommended.
If LWORK = -1, then SLAQR4 does a workspace query. In this case, SLAQR4 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed. ๐
INFOโINFOis INTEGER
INFO is INTEGER
= 0: successful exit
> 0: if INFO = i, SLAQR4 failed to compute all of
the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
and WI contain those eigenvalues which have been
successfully computed. (Failures are rare.)
If INFO > 0 and WANT is .FALSE., then on exit,
the remaining unconverged eigenvalues are the eigen-
values of the upper Hessenberg matrix rows and
columns ILO through INFO of the final, output
value of H.
If INFO > 0 and WANTT is .TRUE., then on exit
(*) (initial value of H)*U = U*(final value of H)
where U is a orthogonal matrix. The final
value of H is upper Hessenberg and triangular in
rows and columns INFO+1 through IHI.
If INFO > 0 and WANTZ is .TRUE., then on exit
(final value of Z(ILO:IHI,ILOZ:IHIZ)
= (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
where U is the orthogonal matrix in (*) (regard-
less of the value of WANTT.)
If INFO > 0 and WANTZ is .FALSE., then Z is not
accessed.
๐ Author ๐
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
๐ซถ Contributors: ๐ซถ
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
๐ References: ๐
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
Performance, SIAM Journal of Matrix Analysis, volume 23, pages
929--947, 2002.
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm
Part II: Aggressive Early Deflation, SIAM Journal of Matrix
Analysis, volume 23, pages 948--973, 2002.
๐ subroutine slaqr5 ๐
subroutine slaqr5 (logical WANTT, logical WANTZ, integer KACC22, integer N,
integer KTOP, integer KBOT, integer NSHFTS, real, dimension( * ) SR,
real, dimension( * ) SI, real, dimension( ldh, * ) H, integer LDH,
integer ILOZ, integer IHIZ, real, dimension( ldz, * ) Z, integer LDZ,
real, dimension( ldv, * ) V, integer LDV, real, dimension( ldu, * ) U,
integer LDU, integer NV, real, dimension( ldwv, * ) WV, integer LDWV,
integer NH, real, dimension( ldwh, * ) WH, integer LDWH)
SLAQR5 performs a single small-bulge multi-shift QR sweep. ๐
โจ Purpose: โจ
SLAQR5, called by SLAQR0, performs a
single small-bulge multi-shift QR sweep.
๐ Parameters ๐
WANTTโWANTTis LOGICAL.WANTT = .true.if the quasi-triangular Schur factor is being computed.WANTTis set to.false.otherwise. ๐WANTZโWANTZis LOGICAL.WANTZ = .true.if the orthogonal Schur factor is being computed.WANTZis set to.false.otherwise. ๐KACC22โKACC22is INTEGER with value 0, 1, or 2. Specifies the computation mode of far-from-diagonal orthogonal updates. ๐
= 0: SLAQR5 does not accumulate reflections and does not
use matrix-matrix multiply to update far-from-diagonal
matrix entries.
= 1: SLAQR5 accumulates reflections and uses matrix-matrix
multiply to update the far-from-diagonal matrix entries.
= 2: Same as KACC22 = 1. This option used to enable exploiting
the 2-by-2 structure during matrix multiplications, but
this is no longer supported.
NโNis INTEGER.Nis the order of the Hessenberg matrixHupon which this subroutine operates. ๐ธKTOPโKTOPis INTEGER ๐ฆKBOTโKBOTis INTEGER
These are the first and last rows and columns of an isolated diagonal block upon which the QR sweep is to be applied. It is assumed without a check that
either KTOP = 1 or H(KTOP,KTOP-1) = 0
and
either KBOT = N or H(KBOT+1,KBOT) = 0.
NSHFTSโNSHFTSis INTEGER.NSHFTSgives the number of simultaneous shifts.NSHFTSmust be positive and even. โจSRโSRisREAL array, dimension(NSHFTS)๐SIโSIisREAL array, dimension(NSHFTS)
SR contains the real parts and SI contains the imaginary parts of the NSHFTS shifts of origin that define the multi-shift QR sweep. On output SR and SI may be reordered. ๐ฅ
HโHisREAL array, dimension(LDH,N)
On input H contains a Hessenberg matrix. On output a multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied to the isolated diagonal block in rows and columns KTOP through KBOT.
LDHโLDHis INTEGER.LDHis the leading dimension ofHjust as declared in the calling procedure.LDH >= MAX(1,N). ๐ILOZโILOZis INTEGER ๐IHIZโIHIZis INTEGER. Specify the rows ofZto which transformations must be applied ifWANTZis.TRUE..1 <= ILOZ <= IHIZ <= NZโZisREAL array, dimension(LDZ,IHIZ)
If WANTZ = .TRUE., then the QR Sweep orthogonal similarity transformation is accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ = .FALSE., then Z is unreferenced. ๐
LDZโLDZis INTEGER.LDAis the leading dimension ofZjust as declared in the calling procedure.LDZ >= N.
ngl bestie it literally saysLDAhere in the original and we are not the ones who get to fix that ๐ no cap,LDAit is.
VโVisREAL array, dimension(LDV,NSHFTS/2)๐งLDVโLDVis INTEGER.LDVis the leading dimension ofVas declared in the calling procedure.LDV >= 3. โญUโUisREAL array, dimension(LDU,2*NSHFTS)๐LDULDUโ INTEGER, bestie. It's the leading dimension ofU, exactly like it was declared in the in the calling subroutine. Constraint:LDU >= 2*NSHFTS.NVโ INTEGER.NVis the number of rows inWVagailable for workspace. Constraint:NV >= 1. (yes it says "agailable", we don't fix other people's typos, we respect them ๐ )WVโ REAL array, dimension(LDWV,2*NSHFTS). It's the workspace, she's just vibing.LDWVโ INTEGER. The leading dimension ofWVas declared in the in the calling subroutine. Constraint:LDWV >= NV.NHโ INTEGER.NHis the number of columns in arrayWHavailable for workspace. Constraint:NH >= 1.WHโ REAL array, dimension(LDWH,NH). More workspace, more slay.LDWHโ INTEGER. Leading dimension ofWHjust as declared in the calling procedure. Constraint:LDWH >= 2*NSHFTS.
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ซถ Contributors
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
Lars Karlsson, Daniel Kressner, and Bruno Lang
Thijs Steel, Department of Computer science, KU Leuven, Belgium
๐ References
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume 23, pages 929--947, 2002.
Lars Karlsson, Daniel Kressner, and Bruno Lang, Optimally packed chains of bulges in multishift QR algorithms. ACM Trans. Math. Softw. 40, 2, Article 12 (February 2014).
๐ฆ subroutine slaqsb
subroutine slaqsb (character UPLO, integer N, integer KD, real, dimension(
ldab, * ) AB, integer LDAB, real, dimension( * ) S, real SCOND, real
AMAX, character EQUED)
Okay so SLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. It's giving "I inherited someone else's math homework and I'm making it look good" energy โจ
Purpose
SLAQSB equilibrates a symmetric band matrix A using the scaling factors in the vector S. That's it, that's the whole job. ๐ธ
Parameters
UPLOโUPLOisCHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrixAis stored.= 'U': Upper triangular= 'L': Lower triangular
NโNis INTEGER. The order of the matrixA. Constraint:N >= 0.KDโKDis INTEGER. The number of super-diagonals of the matrixAifUPLO = 'U', or the number of sub-diagonals ifUPLO = 'L'. Constraint:KD >= 0.ABโABis 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, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T of the band matrix A, in the same storage format as A. glow up complete ๐
LDABโLDABis INTEGER. The leading dimension of the arrayAB. Constraint:LDAB >= KD+1.SโSis REAL array, dimension(N). The scale factors forA.SCONDโSCONDis REAL. Ratio of the smallestS(i)to the largestS(i).AMAXโAMAXis REAL. Absolute value of largest matrix entry.EQUEDโEQUEDisCHARACTER*1. Specifies whether or not equilibration was done.= 'N': No equilibration.= 'Y': Equilibration was done, i.e.,Ahas been replaced bydiag(S) * A * diag(S).
๐ Internal Parameters
THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done.
LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ subroutine slaqsp
subroutine slaqsp (character UPLO, integer N, real, dimension( * ) AP, real,
dimension( * ) S, real SCOND, real AMAX, character EQUED)
SLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. Same energy as slaqsb but she's packed, not banded, and she knows it. ๐
Purpose
SLAQSP equilibrates a symmetric matrix A using the scaling factors in the vector S. Nothing more, nothing less, periodt.
Parameters
UPLOโUPLOisCHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrixAis stored.= 'U': Upper triangular= 'L': Lower triangular
NโNis INTEGER. The order of the matrixA. Constraint:N >= 0.APโAPis 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)*(2n-j)/2) = A(i,j) for j<=i<=n.
On exit, the equilibrated matrix: diag(S) * A * diag(S), in the same storage format as A. she got a whole new look and kept the same storage format, iconic ๐
SโSis REAL array, dimension(N). The scale factors forA.SCONDโSCONDis REAL. Ratio of the smallestS(i)to the largestS(i).AMAXโAMAXis REAL. Absolute value of largest matrix entry.EQUEDโEQUEDisCHARACTER*1. Specifies whether or not equilibration was done.= 'N': No equilibration.= 'Y': Equilibration was done, i.e.,Ahas been replaced bydiag(S) * A * diag(S).
๐ Internal Parameters
THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done.
LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ subroutine slaqtr
subroutine slaqtr (logical LTRAN, logical LREAL, integer N, real, dimension(
ldt, * ) T, integer LDT, real, dimension( * ) B, real W, real SCALE,
real, dimension( * ) X, real, dimension( * ) WORK, integer INFO)
SLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of special form, in real arithmetic. Yes bestie, complex math, but she's doing it in real arithmetic. Main character energy fr fr. ๐ญโจ
Purpose
SLAQTR solves the real quasi-triangular system
op(T)*p = scale*c, if LREAL = .TRUE.
or the complex quasi-triangular systems
op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE.
in real arithmetic, where T is upper quasi-triangular.
If LREAL = .FALSE., then the first diagonal block of T must be 1 by 1, B is the specially structured matrix
B = [ b(1) b(2) ... b(n) ]
[ w ]
[ w ]
[ . ]
[ w ]
op(A) = A or A**T, A**T denotes the transpose of matrix A.
On input, X = [ c ]. On output, X = [ p ].
[ d ] [ q ]
This subroutine is designed for the condition number estimation in routine STRSNA. iykyk ๐ง
Parameters
LTRANโLTRANis LOGICAL. On entry,LTRANspecifies the option of conjugate transpose:= .FALSE., op(T+i*B) = T+i*B,= .TRUE., op(T+i*B) = (T+i*B)**T.
LREALโLREALis LOGICAL. On entry,LREALspecifies the input matrix structure:= .FALSE., the input is complex= .TRUE., the input is real
NโNis INTEGER. On entry,Nspecifies the order ofT+i*B. Constraint:N >= 0.TโTis REAL array, dimension(LDT,N). On entry,Tcontains a matrix in Schur canonical form. IfLREAL = .FALSE., then the first diagonal block ofTmust be 1 by 1.LDTโLDTis INTEGER. The leading dimension of the matrixT. Constraint:LDT >= max(1,N).BโBis REAL array, dimension(N). On entry,Bcontains the elements to form the matrixBas described above. IfLREAL = .TRUE.,Bis not referenced. (she said "don't perceive me" ๐)WโWis REAL. On entry,Wis the diagonal element of the matrixB. IfLREAL = .TRUE.,Wis not referenced.SCALEโSCALEis REAL. On exit,SCALEis the scale factor.XโXis REAL array, dimension(2*N). On entry,Xcontains the right hand side of the system. On exit,Xis overwritten by the solution. she ate and left no crumbs ๐WORKโWORKis REAL array, dimension(N). Workspace, doing her little job, no notes.INFOโINFOis INTEGER. On exit,INFOis set to0: successful exit.1: the some diagonal 1 by 1 block has been perturbed by a small numberSMINto keep nonsingularity.2: the some diagonal 2 by 2 block has been perturbed by a small number inSLALN2to keep nonsingularity.
NOTE: In the interests of speed, this routine does not check the inputs for errors. Bestie. She is trusting you completely. Do not betray her. ๐
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ slar1v ๐
subroutine slar1v (integer N, integer B1, integer BN, real LAMBDA, real,
dimension( * ) D, real, dimension( * ) L, real, dimension( * ) LD, real,
dimension( * ) LLD, real PIVMIN, real GAPTOL, real, dimension( * ) Z,
logical WANTNC, integer NEGCNT, real ZTZ, real MINGMA, integer R,
integer, dimension( * ) ISUPPZ, real NRMINV, real RESID, real RQCORR,
real, dimension( * ) WORK)
SLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - ฮปI. ๐
๐ฏ Purpose
SLAR1V computes the (scaled) r-th column of the inverse of the sumbmatrix in rows B1 through BN of the tridiagonal matrix L D L**T - sigma I. When sigma is close to an eigenvalue, the computed vector is an accurate eigenvector. Usually, r corresponds to the index where the eigenvector is largest in magnitude. ๐ฆ
The following steps accomplish this computation:
- Stationary qd transform,
L D L**T - sigma I = L(+) D(+) L(+)**T, - Progressive qd transform,
L D L**T - sigma I = U(-) D(-) U(-)**T, - Computation of the diagonal elements of the inverse of
L D L**T - sigma Iby combining the above transforms, and choosingras the index where the diagonal of the inverse is (one of the) largest in magnitude. - Computation of the (scaled) r-th column of the inverse using the twisted factorization obtained by combining the top part of the the stationary and the bottom part of the progressive transform.
โจ Parameters
NโNis INTEGER. The order of the matrixL D L**T. ๐ธB1โB1is INTEGER. First index of the submatrix ofL D L**T.BNโBNis INTEGER. Last index of the submatrix ofL D L**T.LAMBDAโLAMBDAis REAL. The shift. In order to compute an accurate eigenvector,LAMBDAshould be a good approximation to an eigenvalue ofL D L**T. ๐LโLis REAL array, dimension(N-1). The(n-1)subdiagonal elements of the unit bidiagonal matrixL, in elements 1 toN-1.DโDis REAL array, dimension(N). Thendiagonal elements of the diagonal matrixD.LDโLDis REAL array, dimension(N-1). Then-1elementsL(i)*D(i).LLDโLLDis REAL array, dimension(N-1). Then-1elementsL(i)*L(i)*D(i).PIVMINโPIVMINis REAL. The minimum pivot in the Sturm sequence.GAPTOLโGAPTOLis REAL. Tolerance that indicates when eigenvector entries are negligible w.r.t. their contribution to the residual. ๐ZโZis REAL array, dimension(N). On input, all entries ofZmust be set to 0. On output,Zcontains the (scaled) r-th column of the inverse. The scaling is such thatZ(R)equals 1.WANTNCโWANTNCis LOGICAL. Specifies whetherNEGCNThas to be computed.NEGCNTโNEGCNTis INTEGER. IfWANTNCis.TRUE.thenNEGCNT= the number of pivots< pivminin the matrix factorizationL D L**T, andNEGCNT = -1otherwise.ZTZโZTZis REAL. The square of the 2-norm ofZ.MINGMAโMINGMAis REAL. The reciprocal of the largest (in magnitude) diagonal element of the inverse ofL D L**T - sigma I.RโRis INTEGER. The twist index for the twisted factorization used to computeZ. On input,0 <= R <= N. IfRis input as0,Ris set to the index where(L D L**T - sigma I)^{-1}is largest in magnitude. If1 <= R <= N,Ris unchanged. On output,Rcontains the twist index used to computeZ. Ideally,Rdesignates the position of the maximum entry in the eigenvector. ๐ISUPPZโISUPPZis INTEGER array, dimension(2). The support of the vector inZ, i.e., the vectorZis nonzero only in elementsISUPPZ(1)throughISUPPZ( 2 ).NRMINVโNRMINVis REAL.NRMINV = 1/SQRT( ZTZ )RESIDโRESIDis REAL. The residual of the FP vector.RESID = ABS( MINGMA )/SQRT( ZTZ )RQCORRโRQCORRis REAL. The Rayleigh Quotient correction toLAMBDA.RQCORR = MINGMA*TMPWORKโWORKis REAL array, dimension(4*N)
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ซถ Contributors
- Beresford Parlett, University of California, Berkeley, USA
- Jim Demmel, University of California, Berkeley, USA
- Inderjit Dhillon, University of Texas, Austin, USA
- Osni Marques, LBNL/NERSC, USA
- Christof Voemel, University of California, Berkeley, USA
๐ slar2v ๐
subroutine slar2v (integer N, real, dimension( * ) X, real, dimension( * )
Y, real, dimension( * ) Z, integer INCX, real, dimension( * ) C, real,
dimension( * ) S, integer INCC)
SLAR2V applies a vector of plane rotations with real cosines and real sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. โจ
๐ฏ Purpose
SLAR2V applies a vector of real plane rotations from both sides to a sequence of 2-by-2 real symmetric matrices, defined by the elements of the vectors x, y and z. For i = 1,2,...,n
( x(i) z(i) ) := ( c(i) s(i) ) ( x(i) z(i) ) ( c(i) -s(i) )
( z(i) y(i) ) ( -s(i) c(i) ) ( z(i) y(i) ) ( s(i) c(i) )
โจ Parameters
NโNis INTEGER. The number of plane rotations to be applied. ๐ธXโXis REAL array, dimension(1+(N-1)*INCX). The vectorx.YโYis REAL array, dimension(1+(N-1)*INCX). The vectory.ZโZis REAL array, dimension(1+(N-1)*INCX). The vectorz.INCXโINCXis INTEGER. The increment between elements ofX,YandZ.INCX > 0.CโCis REAL array, dimension(1+(N-1)*INCC). The cosines of the plane rotations.SโSis REAL array, dimension(1+(N-1)*INCC). The sines of the plane rotations.INCCโINCCis INTEGER. The increment between elements ofCandS.INCC > 0. ๐
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ฆ slarf ๐ฆ
subroutine slarf (character SIDE, integer M, integer N, real, dimension( * )
V, integer INCV, real TAU, real, dimension( ldc, * ) C, integer LDC,
real, dimension( * ) WORK)
SLARF applies an elementary reflector to a general rectangular matrix. ๐
๐ฏ Purpose
SLARF applies a real elementary reflector H to a real m by n matrix C, from either the left or the right. H is represented in the form
H = I - tau * v * v**T
where tau is a real scalar and v is a real vector.
If tau = 0, then H is taken to be the unit matrix. โญ
โจ Parameters
SIDEโSIDEis CHARACTER*1= 'L': formH * C= 'R': formC * H
MโMis INTEGER. The number of rows of the matrixC.NโNis INTEGER. The number of columns of the matrixC.VโVis REAL array, dimension(1 + (M-1)*abs(INCV))ifSIDE = 'L'- or
(1 + (N-1)*abs(INCV))ifSIDE = 'R'
The vector v in the representation of H. V is not used if TAU = 0. ๐
INCVโINCVis INTEGER. The increment between elements ofv.INCV <> 0.TAUโTAUis REAL. The valuetauin the representation ofH.CโCis REAL array, dimension(LDC,N). On entry, thembynmatrixC. On exit,Cis overwritten by the matrixH * CifSIDE = 'L', orC * HifSIDE = 'R'.LDCโLDCis INTEGER. The leading dimension of the arrayC.LDC >= max(1,M). ๐ฅWORKโWORKis REAL array, dimension(N)ifSIDE = 'L'- or
(M)ifSIDE = 'R'
๐ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ slarfb ๐
subroutine slarfb (character SIDE, character TRANS, character DIRECT,
character STOREV, integer M, integer N, integer K, real, dimension( ldv,
* ) V, integer LDV, real, dimension( ldt, * ) T, integer LDT, real,
dimension( ldc, * ) C, integer LDC, real, dimension( ldwork, * ) WORK,
integer LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix. ๐
๐ฏ Purpose
SLARFB applies a real block reflector H or its transpose H**T to a real m by n matrix C, from either the left or the right. โจ
โจ Parameters
SIDEโSIDEis CHARACTER*1= 'L': applyHorH**Tfrom the Left= 'R': applyHorH**Tfrom the Right`
TRANSโTRANSis CHARACTER*1= 'N': applyH(No transpose)= 'T': applyH**T(Transpose)
DIRECTโDIRECTis CHARACTER*1. Indicates howHis formed from a product of elementary reflectors= 'F':H = H(1) H(2) . . . H(k)(Forward)= 'B':H = H(k) . . . H(2) H(1)(Backward)
STOREVโSTOREVis CHARACTER*1. Indicates how the vectors which define the elementary reflectors are stored:= 'C': Columnwise= 'R': Rowwise
MโMis INTEGER. The number of rows of the matrixC.NโNis INTEGER. The number of columns of the matrixC.KโKโKis INTEGER. The order of the matrixT(which is literally just the number of elementary reflectors whose product defines the block reflector, no biggie ๐ ). IfSIDE = 'L', thenM >= K >= 0; ifSIDE = 'R', thenN >= K >= 0.VโVis REAL array, dimension
(LDV,K) if STOREV = 'C'
(LDV,M) if STOREV = 'R' and SIDE = 'L'
(LDV,N) if STOREV = 'R' and SIDE = 'R'
This is the matrix V, bestie. See Further Details for the whole vibe โจ
LDVโLDVis INTEGER, the leading dimension of the arrayV. IfSTOREV = 'C'andSIDE = 'L',LDV >= max(1,M); ifSTOREV = 'C'andSIDE = 'R',LDV >= max(1,N); ifSTOREV = 'R',LDV >= K.TโTis REAL array, dimension(LDT,K). The triangularkbykmatrixTin the representation of the block reflector. Iconic honestly ๐LDTโLDTis INTEGER, the leading dimension of the arrayT.LDT >= K.CโCis REAL array, dimension(LDC,N). On entry, thembynmatrixC. On exit,Cis overwritten byH*CorH**T*CorC*HorC*H**T. She gets rewritten, that's just her journey ๐ฆLDCโLDCis INTEGER, the leading dimension of the arrayC.LDC >= max(1,M).WORKโWORKis REAL array, dimension(LDWORK,K).LDWORKโLDWORKis INTEGER, the leading dimension of the arrayWORK. IfSIDE = 'L',LDWORK >= max(1,N); ifSIDE = 'R',LDWORK >= max(1,M).
โจ Author โจ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
๐ Further Details: ๐
The shape of the matrix V and the storage of the vectors which define the H(i) is best illustrated by the following example with n = 5 and k = 3. The elements equal to 1 are not stored; the corresponding array elements are modified but restored on exit. The rest of the array is not used. (Not stored, but still getting touched and then put back โ the audacity ๐)
DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
( v1 1 ) ( 1 v2 v2 v2 )
( v1 v2 1 ) ( 1 v3 v3 )
( v1 v2 v3 )
( v1 v2 v3 )
DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
V = ( v1 v2 v3 ) V = ( v1 v1 1 )
( v1 v2 v3 ) ( v2 v2 v2 1 )
( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
( 1 v3 )
( 1 )
๐ซถ subroutine slarfg ๐ซถ
subroutine slarfg (integer N, real ALPHA, real, dimension( * ) X, integer INCX, real TAU)
SLARFG generates an elementary reflector (Householder matrix). Okay so she's building you a Householder matrix, very demure, very mindful ๐
๐ Purpose: ๐
SLARFG generates a real elementary reflector H of order n, such that
H * ( alpha ) = ( beta ), H**T * H = I.
( x ) ( 0 )
where alpha and beta are scalars, and x is an (n-1)-element real vector. H is represented in the form
H = I - tau * ( 1 ) * ( 1 v**T ) ,
( v )
where tau is a real scalar and v is a real (n-1)-element vector.
If the elements of x are all zero, then tau = 0 and H is taken to be the unit matrix. (Nothing to see here, she's just the identity, respectfully ๐)
Otherwise 1 <= tau <= 2.
๐ธ Parameters ๐ธ
NโNis INTEGER. The order of the elementary reflector.ALPHAโALPHAis REAL. On entry, the valuealpha. On exit, it is overwritten with the valuebeta. Glow up, honestly โจXโXis REAL array, dimension
(1+(N-2)*abs(INCX))
On entry, the vector x. On exit, it is overwritten with the vector v.
INCXโINCXis INTEGER. The increment between elements ofX.INCX > 0. Periodt.TAUโTAUis REAL. The valuetau.
โจ Author โจ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
๐ฆ subroutine slarfgp ๐ฆ
subroutine slarfgp (integer N, real ALPHA, real, dimension( * ) X, integer INCX, real TAU)
SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. Same as before but the beta is non-negative โ she said no negativity in this household ๐
๐ Purpose: ๐
SLARFGP generates a real elementary reflector H of order n, such that
H * ( alpha ) = ( beta ), H**T * H = I.
( x ) ( 0 )
where alpha and beta are scalars, beta is non-negative, and x is an (n-1)-element real vector. H is represented in the form
H = I - tau * ( 1 ) * ( 1 v**T ) ,
( v )
where tau is a real scalar and v is a real (n-1)-element vector.
If the elements of x are all zero, then tau = 0 and H is taken to be the unit matrix.
๐ธ Parameters ๐ธ
NโNis INTEGER. The order of the elementary reflector.ALPHAโALPHAis REAL. On entry, the valuealpha. On exit, it is overwritten with the valuebeta.XโXis REAL array, dimension
(1+(N-2)*abs(INCX))
On entry, the vector x. On exit, it is overwritten with the vector v.
INCXโINCXis INTEGER. The increment between elements ofX.INCX > 0.TAUโTAUis REAL. The valuetau.
โจ Author โจ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
๐ subroutine slarft ๐
subroutine slarft (character DIRECT, character STOREV, integer N, integer K, real, dimension( ldv, * ) V, integer LDV, real, dimension( * ) TAU, real, dimension( ldt, * ) T, integer LDT)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH โ she's assembling the T, main character energy fr fr ๐
๐ Purpose: ๐
SLARFT forms the triangular factor T of a real block reflector H of order n, which is defined as a product of k elementary reflectors.
If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
If STOREV = 'C', the vector which defines the elementary reflector H(i) is stored in the i-th column of the array V, and
H = I - V * T * V**T
If STOREV = 'R', the vector which defines the elementary reflector H(i) is stored in the i-th row of the array V, and
H = I - V**T * T * V
๐ธ Parameters ๐ธ
DIRECTโDIRECTis CHARACTER*1. Specifies the order in which the elementary reflectors are multiplied to form the block reflector:= 'F':H = H(1) H(2) . . . H(k)(Forward)= 'B':H = H(k) . . . H(2) H(1)(Backward)
STOREVโSTOREVis CHARACTER*1. Specifies how the vectors which define the elementary reflectors are stored (see also Further Details):= 'C': columnwise= 'R': rowwise
NโNis INTEGER. The order of the block reflectorH.N >= 0.KโKis INTEGER. The order of the triangular factorT(= the number of elementary reflectors).K >= 1.VโVis REAL array, dimension
(LDV,K) if STOREV = 'C'
(LDV,N) if STOREV = 'R'
The matrix V. See further details.
LDVโLDVis INTEGER. The leading dimension of the arrayV. IfSTOREV = 'C',LDV >= max(1,N); ifSTOREV = 'R',LDV >= K.TAUโTAUis REAL array, dimension(K).TAU(i)must contain the scalar factor of the elementary reflectorH(i).TโTis REAL array, dimension(LDT,K). Thekbyktriangular factorTof the block reflector. IfDIRECT = 'F',Tis upper triangular; ifDIRECT = 'B',Tis lower triangular. The rest of the array is not used.LDTโLDTis INTEGER. The leading dimension of the arrayT.LDT >= K.
โจ Author โจ
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
๐ Further Details: ๐
The shape of the matrix V and the storage of the vectors which define the H(i) is best illustrated by the following example with n = 5 and k = 3. The elements equal to 1 are not stored.
DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
( v1 1 ) ( 1 v2 v2 v2 )
( v1 v2 1 ) ( 1 v3 v3 )
( v1 v2 v3 )
( v1 v2 v3 )
DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
V = ( v1 v2 v3 ) V = ( v1 v1 1 )
( v1 v2 v3 ) ( v2 v2 v2 1 )
( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
( 1 v3 )
( 1 )
โญ subroutine slarfx โญ
subroutine slarfx (character SIDE, integer M, integer N, real, dimension( * ) V, real TAU, real, dimension( ldc, * ) C, integer LDC, real, dimension( * ) WORK)
SLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order โค 10. She's got moves โ unrolls that loop when the reflector's order is โค 10 ๐ง
๐ Purpose: ๐
SLARFX applies a real elementary reflector H to a real m by n matrix C, from either the left or the right. H is represented in the form
H = I - tau * v * v**T
where tau is a real scalar and v is a real vector.
If tau = 0, then H is taken to be the unit matrix
This version uses inline code if H has order < 11.
๐ธ Parameters ๐ธ
SIDEโSIDEis CHARACTER*1= 'L': formH * C= 'R': formC * H
MโMis INTEGER. The number of rows of the matrixC.NNโN is INTEGER. It's the number of columns of the matrixC. ๐VโV is REAL array, dimension(M)ifSIDE = 'L', or(N)ifSIDE = 'R'. It's the vectorvin the representation ofH. โจTAUโTAU is REAL. It's the valuetauin the representation ofH.CโC is REAL array, dimension(LDC,N). On entry, thembynmatrixC. On exit,Cis overwritten by the matrixH * CifSIDE = 'L', orC * HifSIDE = 'R'. ๐LDCโLDC is INTEGER. The leading dimension of the arrayC.LDC >= (1,M).WORKโWORK is REAL array, dimension(N)ifSIDE = 'L', or(M)ifSIDE = 'R'. And heads up bestie โWORKis not referenced ifHhas order< 11. ๐
Author ๐ซถ
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
โจ SLARFY โจ
subroutine slarfy (character UPLO, integer N, real, dimension( * ) V,
integer INCV, real TAU, real, dimension( ldc, * ) C, integer LDC, real,
dimension( * ) WORK)
Purpose ๐ฆ
SLARFY applies an elementary reflector, or Householder matrix, H, to an n x n symmetric matrix C, from both the left and the right. It's giving both sides, periodt.
H is represented in the form
H = I - tau * v * v'
where tau is a scalar and v is a vector.
If tau is zero, then H is taken to be the unit matrix. No notes. ๐
Parameters ๐ธ
UPLOโUPLO is CHARACTER*1. Specifies whether the upper or lower triangular part of the symmetric matrixCis stored.= 'U': Upper triangle= 'L': Lower triangle
NโN is INTEGER. The number of rows and columns of the matrixC.N >= 0.VโV is REAL array, dimension(1 + (N-1)*abs(INCV)). The vectorvas described above. ๐INCVโINCV is INTEGER. The increment between successive elements ofv.INCVmust not be zero. Like, fr fr, don't.TAUโTAU is REAL. The valuetauas described above.CโC is REAL array, dimension(LDC, N). On entry, the matrixC. On exit,Cis overwritten byH * C * H'. ๐ it's gone, that's the whole point.LDCโLDC is INTEGER. The leading dimension of the arrayC.LDC >= max( 1, N ).WORKโWORK is REAL array, dimension(N).
Author ๐ซถ
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
โจ SLARGV โจ
subroutine slargv (integer N, real, dimension( * ) X, integer INCX, real,
dimension( * ) Y, integer INCY, real, dimension( * ) C, integer INCC)
SLARGV generates a vector of plane rotations with real cosines and real sines. Very demure, very rotational. ๐
Purpose ๐ซ
SLARGV generates a vector of real plane rotations, determined by elements of the real vectors x and y. For i = 1,2,...,n
( c(i) s(i) ) ( x(i) ) = ( a(i) )
( -s(i) c(i) ) ( y(i) ) = ( 0 )
Parameters ๐
NโN is INTEGER. The number of plane rotations to be generated.XโX is REAL array, dimension(1+(N-1)*INCX). On entry, the vectorx. On exit,x(i)is overwritten bya(i), fori = 1,...,n. ๐INCXโINCX is INTEGER. The increment between elements ofX.INCX > 0.YโY is REAL array, dimension(1+(N-1)*INCY). On entry, the vectory. On exit, the sines of the plane rotations.INCYโINCY is INTEGER. The increment between elements ofY.INCY > 0.CโC is REAL array, dimension(1+(N-1)*INCC). The cosines of the plane rotations. ๐INCCโINCC is INTEGER. The increment between elements ofC.INCC > 0.
Author ๐ซถ
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
โจ SLARRV โจ
subroutine slarrv (integer N, real VL, real VU, real, dimension( * ) D,
real, dimension( * ) L, real PIVMIN, integer, dimension( * ) ISPLIT,
integer M, integer DOL, integer DOU, real MINRGP, real RTOL1, real
RTOL2, real, dimension( * ) W, real, dimension( * ) WERR, real,
dimension( * ) WGAP, integer, dimension( * ) IBLOCK, integer, dimension(
* ) INDEXW, real, dimension( * ) GERS, real, dimension( ldz, * ) Z,
integer LDZ, integer, dimension( * ) ISUPPZ, real, dimension( * ) WORK,
integer, dimension( * ) IWORK, integer INFO)
SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. Bestie this one is the main character. ๐
Purpose ๐
SLARRV computes the eigenvectors of the tridiagonal matrix T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. The input eigenvalues should have been computed by SLARRE.
Parameters ๐
NโN is INTEGER. The order of the matrix.N >= 0.VLโVL is REAL. Lower bound of the interval that contains the desired eigenvalues.VL < VU. Needed to compute gaps on the left or right end of the extremal eigenvalues in the desired RANGE.VUโVU is REAL. Upper bound of the interval that contains the desired eigenvalues.VL < VU.
Note:VUis currently not used by this implementation ofSLARRV.VUis passed toSLARRVbecause it could be used to compute gaps on the right end of the extremal eigenvalues. However, with not much initial accuracy inLAMBDAandVU, the formula can lead to an overestimation of the right gap and thus to inadequately early RQI 'convergence'. This is currently prevented by forcing a small right gap. And so it turns out thatVUis currently not used by this implementation ofSLARRV. ๐ (iconic honestly, it's there for the plot)
DโD is REAL array, dimension(N). On entry, theNdiagonal elements of the diagonal matrixD. On exit,Dmay be overwritten.LโL is REAL array, dimension(N). On entry, the(N-1)subdiagonal elements of the unit bidiagonal matrixLare in elements1toN-1ofL(if the matrix is not split). At the end of each block is stored the corresponding shift as given bySLARRE. On exit,Lis overwritten. ๐PIVMINโPIVMIN is REAL. The minimum pivot allowed in the Sturm sequence.ISPLITโISPLIT is INTEGER array, dimension(N). The splitting points, at whichTbreaks up into blocks. The first block consists of rows/columns1toISPLIT( 1 ), the second of rows/columnsISPLIT( 1 )+1throughISPLIT( 2 ), etc. ๐ธMโM is INTEGER. The total number of input eigenvalues.0 <= M <= N.DOLโDOL is INTEGER.DOUโDOU is INTEGER. If the user wants to compute only selected eigenvectors from all the eigenvalues supplied, he can specify an index rangeDOL:DOU. Or else the settingDOL=1,DOU=Mshould be applied. Note thatDOLandDOUrefer to the order in which the eigenvalues are stored inW. If the user wants to compute only selected eigenpairs, then the columnsDOL-1toDOU+1of the eigenvector spaceZcontain the computed eigenvectors. All other columns ofZare set to zero. โจMINRGPโMINRGP is REAL.RTOL1โRTOL1 is REAL.RTOL2โRTOL2 is REAL. Parameters for bisection. An interval[LEFT,RIGHT]has converged if
RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
WโW is REAL array, dimension(N). The firstMelements ofWcontain the APPROXIMATE eigenvalues for which eigenvectors are to be computed. The eigenvalues should be grouped by split-off block and ordered from smallest to largest within the block (the output arrayWfromSLARREis expected here). Furthermore, they are with respect to the shift of the corresponding root representation for their block. On exit,Wholds the eigenvalues of the UNshifted matrix. ๐ฆWERRโWERR is REAL array, dimension(N). The firstMelements contain the semiwidth of the uncertainty interval of the corresponding eigenvalue inW.WGAPโWGAP is REAL array, dimension(N). The separation from the right neighbor eigenvalue inW.IBLOCKโIBLOCK is INTEGER array, dimension(N). The indices of the blocks (submatrices) associated with the corresponding eigenvalues inW;IBLOCK(i)=1if eigenvalueW(i)belongs to the first block from the top,=2ifW(i)belongs to the second block, etc.INDEXWโINDEXW is INTEGER array, dimension(N). The indices of the eigenvalues within each block (submatrix); for example,INDEXW(i)= 10andIBLOCK(i)=2imply that thei-th eigenvalueW(i)is the10-th eigenvalue in the second block. ๐GERSโGERS is REAL array, dimension(2*N). TheNGerschgorin intervals (thei-th Gerschgorin interval is(GERS(2*i-1), GERS(2*i))). The Gerschgorin intervals should be computed from the original UNshifted matrix.ZโZ is REAL array, dimension(LDZ, max(1,M) ). IfINFO = 0, the firstMcolumns ofZcontain the orthonormal eigenvectors of the matrixTcorresponding to the input eigenvalues, with thei-th column ofZholding the eigenvector associated withW(i).
Note: the user must ensure that at leastmax(1,M)columns are supplied in the arrayZ. ๐ซถ
LDZโ and the source text justโฆ stops right here mid-entry, no description, nothing. Not me inventing one, bestie. ๐ Onward!LDZโ okay bestie,LDZis anINTEGER. It's the leading dimension of the arrayZ. The rules?LDZ >= 1, and ifJOBZ = 'V', thenLDZ >= max(1,N). Slay. ๐ISUPPZโISUPPZis anINTEGERarray, dimension( 2*max(1,M) ). It's the support of the eigenvectors inZ, aka the indices that tell you which elements inZare nonzero. The I-th eigenvector is nonzero only in elementsISUPPZ( 2*I-1 )throughISUPPZ( 2*I ). Iconic. โจWORKโWORKis aREALarray, dimension(12*N). Get that workspace, bestie. ๐ ๏ธIWORKโIWORKis anINTEGERarray, dimension(7*N). Integer workspace, no cap. ๐งฎINFOโINFOis anINTEGER. Here's the tea, fr fr:= 0: successful exit. We love to see it. ๐> 0: A problem occurred inSLARRV. ๐ญ< 0: One of the called subroutines signaled an internal problem. Needs inspection of the corresponding parameterIINFOfor further information. ๐=-1: Problem inSLARRBwhen refining a child's eigenvalues. ๐=-2: Problem inSLARRFwhen computing the RRR of a child. When a child is inside a tight cluster, it can be difficult to find an RRR. A partial remedy from the user's point of view is to make the parameterMINRGPsmaller and recompile. However, as the orthogonality of the computed vectors is proportional to1/MINRGP, the user should be aware that he might be trading in precision when he decreasesMINRGP. > โ ๏ธ Bestie, this is a precision trade-off, so be careful!=-3: Problem inSLARRBwhen refining a single eigenvalue after the Rayleigh correction was rejected. ๐ต= 5: The Rayleigh Quotient Iteration failed to converge to full accuracy inMAXITRsteps. ๐ซ
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Contributors:
- Beresford Parlett, University of California, Berkeley, USA
- Jim Demmel, University of California, Berkeley, USA
- Inderjit Dhillon, University of Texas, Austin, USA
- Osni Marques, LBNL/NERSC, USA
- Christof Voemel, University of California, Berkeley, USA
subroutine slartv (integer N, real, dimension( * ) X, integer INCX, real, dimension( * ) Y, integer INCY, real, dimension( * ) C, real, dimension( * ) S, integer INCC)
SLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair of vectors. It's giving rotation energy. ๐
Purpose
SLARTV applies a vector of real plane rotations to elements of the real vectors x and y. For i = 1,2,...,n
( x(i) ) := ( c(i) s(i) ) ( x(i) )
( y(i) ) ( -s(i) c(i) ) ( y(i) )
Parameters
NโNisINTEGER. The number of plane rotations to be applied. ๐ขXโXisREALarray, dimension(1+(N-1)*INCX). The vectorx. ๐INCXโINCXisINTEGER. The increment between elements ofX.INCX > 0. ๐ถYโYisREALarray, dimension(1+(N-1)*INCY). The vectory. ๐INCYโINCYisINTEGER. The increment between elements ofY.INCY > 0. ๐CโCisREALarray, dimension(1+(N-1)*INCC). The cosines of the plane rotations. ๐SโSisREALarray, dimension(1+(N-1)*INCC). The sines of the plane rotations. ๐ฆINCCโINCCisINTEGER. The increment between elements ofCandS.INCC > 0. ๐ซ
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
subroutine slaswp (integer N, real, dimension( lda, * ) A, integer LDA, integer K1, integer K2, integer, dimension( * ) IPIV, integer INCX)
SLASWP performs a series of row interchanges on a general rectangular matrix. Row swap slay. ๐
Purpose
SLASWP performs a series of row interchanges on the matrix A. One row interchange is initiated for each of rows K1 through K2 of A. ๐
Parameters
NโNisINTEGER. The number of columns of the matrixA. ๐AโAisREALarray, dimension(LDA,N). On entry, the matrix of column dimensionNto which the row interchanges will be applied. On exit, the permuted matrix. ๐งฉLDAโLDAisINTEGER. The leading dimension of the arrayA. ๐K1โK1isINTEGER. The first element ofIPIVfor which a row interchange will be done. ๐ฅK2โK2isINTEGER.(K2-K1+1)is the number of elements ofIPIVfor which a row interchange will be done. ๐ฅIPIVโIPIVisINTEGERarray, dimension(K1+(K2-K1)*abs(INCX)). The vector of pivot indices. Only the elements in positionsK1throughK1+(K2-K1)*abs(INCX)ofIPIVare accessed.IPIV(K1+(K-K1)*abs(INCX)) = Limplies rowsKandLare to be interchanged. ๐INCXโINCXisINTEGER. The increment between successive values ofIPIV. IfINCXis negative, the pivots are applied in reverse order. โฉ๏ธ
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Further Details
Modified by R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA
subroutine slatbs (character UPLO, character TRANS, character DIAG, character NORMIN, integer N, integer KD, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)
SLATBS solves a triangular banded system of equations. It's giving linear algebra realness. ๐ฏ
Purpose
SLATBS solves one of the triangular systems
A *x = s*b or A**T*x = s*b
with scaling to prevent overflow, where A is an upper or lower triangular band matrix. Here A**T denotes the transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine STBSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned. ๐ฅ
Parameters
UPLOโUPLOisCHARACTER*1. Specifies whether the matrixAis upper or lower triangular.= 'U': Upper triangular= 'L': Lower triangular
TRANSโTRANSisCHARACTER*1. Specifies the operation applied toA.= 'N': SolveA * x = s*b(No transpose)= 'T': SolveA**T* x = s*b(Transpose)= 'C': SolveA**T* x = s*b(Conjugate transpose = Transpose)
DIAGโDIAGisCHARACTER*1. Specifies whether or not the matrixAis unit triangular.= 'N': Non-unit triangular= 'U': Unit triangular
NORMINโNORMINisCHARACTER*1. Specifies whetherCNORMhas been set or not.= 'Y':CNORMcontains the column norms on entry= 'N':CNORMis not set on entry. On exit, the norms will be computed and stored inCNORM.
NโNisINTEGER. The order of the matrixA.N >= 0.KDโKDisINTEGER. The number of subdiagonals or superdiagonals in the triangular matrixA.KD >= 0.ABโABisREALarray, dimension(LDAB,N). The upper or lower triangular band matrixA, stored in the firstKD+1rows of the array. The j-th column ofAis stored in the j-th column of the arrayABas follows: ifUPLO = 'U',AB(kd+1+i-j,j) = A(i,j)formax(1,j-kd)<=i<=j; ifUPLO = 'L',AB(1+i-j,j) = A(i,j)forj<=i<=min(n,j+kd).LDABโLDABisINTEGER. The leading dimension of the arrayAB.LDAB >= KD+1.XโXisREALarray, dimension(N). On entry, the right hand sidebof the triangular system. On exit,Xis overwritten by the solution vectorx.SCALEโSCALEisREAL. The scaling factorsfor the triangular system
A * x = s*b or A**T* x = s*b.
If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.
CNORMโCNORMisREALarray, dimension(N). IfNORMIN = 'Y',CNORMis an input argument andCNORM(j)contains the norm of the off-diagonal part of the j-th column ofA. IfTRANS = 'N',CNORM(j)must be greater than or equal to the infinity-norm, and ifTRANS = 'T'or'C',CNORM(j)must be greater than or equal to the 1-norm. IfNORMIN = 'N',CNORMis an output argument andCNORM(j)returns the 1-norm of the offdiagonal part of the j-th column ofA.INFOโINFOisINTEGER.= 0: successful exit< 0: ifINFO = -k, the k-th argument had an illegal value
Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
Further Details
A rough bound on x is computed; if that is less than overflow, STBSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation.
A columnwise scheme is used for solving A*x = b. The basic algorithm if A is lower triangular is
x[1:n] := b[1:n]
for j = 1, ..., n
x(j) := x(j) / A(j,j)
x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
end
Define bounds on the components of x after j iterations of the loop: M(j) = bound on x[1:j] G(j) = bound on x[j+1:n] Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
okay so picking up right where we were, bestie โ the bound for iteration j+1 is literally giving:
Then for iteration j+1 we have
M(j+1) <= G(j) / | A(j+1,j+1) |
G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
<= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
where CNORM(j+1) is greater than or equal to the infinity-norm of
column j+1 of A, not counting the diagonal. Hence
G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
1<=i<=j
and
|x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
1<=i< j
Since |x(j)| <= M(j), we use the Level 2 BLAS routine STBSV if the
reciprocal of the largest M(j), j=1,..,n, is larger than
max(underflow, 1/overflow).
The bound on x(j) is also used to determine when a step in the
columnwise method can be performed without fear of overflow. If
the computed bound is greater than a large constant, x is scaled to
prevent overflow, but if the bound overflows, x is set to 0, x(j) to
1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
Similarly, a row-wise scheme is used to solve A**T*x = b. The basic
algorithm for A upper triangular is
for j = 1, ..., n
x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
end
We simultaneously compute two bounds
G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
M(j) = bound on x(i), 1<=i<=j
The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
Then the bound on x(j) is
M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
<= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
1<=i<=j
and we can safely call STBSV if 1/M(n) and 1/G(n) are both greater
than max(underflow, 1/overflow).
and that's the whole vibe of the bound. โจ
๐ subroutine slatdf (integer IJOB, integer N, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) RHS, real RDSUM, real RDSCAL, integer, dimension( * ) IPIV, integer, dimension( * ) JPIV)
SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. ๐งฎ
Purpose:
SLATDF takes the LU factorization of the n-by-n matrix Z that SGETC2 computed, and it computes a contribution to the reciprocal Dif-estimate. How? By solving Z * x = b for x, and choosing the r.h.s. b so that the norm of x is as large as possible. On entry RHS = b holds the contribution from earlier solved sub-systems, and on return RHS = x. ๐
The factorization of Z returned by SGETC2 has the form Z = P*L*U*Q, where P and Q are permutation matrices. L is lower triangular with unit diagonal elements and U is upper triangular.
Parameters
IJOBโIJOBis INTEGERIJOB = 2: First compute an approximative null-vector e of Z using SGECON, e is normalized and solve forZx = +-e - fwith the sign giving the greater value of 2-norm(x). About 5 times as expensive as Default.IJOB .ne. 2: Local look ahead strategy where all entries of the r.h.s. b is chosen as either +1 or -1 (Default).NโNis INTEGER The number of columns of the matrix Z.ZโZis REAL array, dimension (LDZ,N) On entry, the LU part of the factorization of the n-by-n matrix Z computed by SGETC2:Z = P * L * U * QLDZโLDZis INTEGER The leading dimension of the array Z.LDA >= max(1, N).RHSโRHSis REAL array, dimension N. On entry, RHS contains contributions from other subsystems. On exit, RHS contains the solution of the subsystem with entries according to the value of IJOB (see above).RDSUMโRDSUMis REAL On entry, the sum of squares of computed contributions to the Dif-estimate under computation by STGSYL, where the scaling factor RDSCAL (see below) has been factored out. On exit, the corresponding sum of squares updated with the contributions from the current sub-system. IfTRANS = 'T'RDSUM is not touched. NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL.RDSCALโRDSCALis REAL On entry, scaling factor used to prevent overflow in RDSUM. On exit, RDSCAL is updated w.r.t. the current contributions in RDSUM. IfTRANS = 'T', RDSCAL is not touched. NOTE: RDSCAL only makes sense when STGSY2 is called by STGSYL.IPIVโIPIVis INTEGER array, dimension (N). The pivot indices; for1 <= i <= N, row i of the matrix has been interchanged with rowIPIV(i).JPIVโJPIVis INTEGER array, dimension (N). The pivot indices; for1 <= j <= N, column j of the matrix has been interchanged with columnJPIV(j).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
This routine is a further developed implementation of algorithm BSOLVE in [1] using complete pivoting in the LU factorization.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
[1] Bo Kagstrom and Lars Westin,
Generalized Schur Methods with Condition Estimators for
Solving the Generalized Sylvester Equation, IEEE Transactions
on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
[2] Peter Poromaa,
On Efficient and Robust Estimators for the Separation
between two Regular Matrix Pairs with Applications in
Condition Estimation. Report IMINF-95.05, Departement of
Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.
๐ธ subroutine slatps (character UPLO, character TRANS, character DIAG, character NORMIN, integer N, real, dimension( * ) AP, real, dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)
SLATPS solves a triangular system of equations with the matrix held in packed storage. ๐ฆ
Purpose:
SLATPS solves one of the triangular systems A *x = s*b or A**T*x = s*b with scaling to prevent overflow, where A is an upper or lower triangular matrix stored in packed form. Here A**T denotes the transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine STPSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
Parameters
UPLOโUPLOis CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangularTRANSโTRANSis CHARACTER*1 Specifies the operation applied to A. = 'N': SolveA * x = s*b(No transpose) = 'T': SolveA**T* x = s*b(Transpose) = 'C': SolveA**T* x = s*b(Conjugate transpose = Transpose)DIAGโDIAGis CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangularNORMINโNORMINis CHARACTER*1 Specifies whether CNORM has been set or not. = 'Y': CNORM contains the column norms on entry = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.NโNis INTEGER The order of the matrix A.N >= 0.APโAPis REAL array, dimension (N*(N+1)/2) The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: ifUPLO = 'U',AP(i + (j-1)*j/2) = A(i,j)for1<=i<=j; ifUPLO = 'L',AP(i + (j-1)*(2n-j)/2) = A(i,j)forj<=i<=n.XโXis REAL array, dimension (N) On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.SCALEโSCALEis REAL The scaling factor s for the triangular systemA * x = s*borA**T* x = s*b. IfSCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution toA*x = 0.CNORMโCNORMis REAL array, dimension (N) IfNORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. IfTRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and ifTRANS = 'T'or'C', CNORM(j) must be greater than or equal to the 1-norm. IfNORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.INFOโINFOis INTEGER = 0: successful exit < 0: ifINFO = -k, the k-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
A rough bound on x is computed; if that is less than overflow, STPSV
is called, otherwise, specific code is used which checks for possible
overflow or divide-by-zero at every operation.
A columnwise scheme is used for solving A*x = b. The basic algorithm
if A is lower triangular is
x[1:n] := b[1:n]
for j = 1, ..., n
x(j) := x(j) / A(j,j)
x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
end
Define bounds on the components of x after j iterations of the loop:
M(j) = bound on x[1:j]
G(j) = bound on x[j+1:n]
Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
Then for iteration j+1 we have
M(j+1) <= G(j) / | A(j+1,j+1) |
G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
<= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
where CNORM(j+1) is greater than or equal to the infinity-norm of column j+1 of A, not counting the diagonal. Hence
G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
1<=i<=j
and
|x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
1<=i< j
Since |x(j)| <= M(j), we use the Level 2 BLAS routine STPSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).
The bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. If the computed bound is greater than a large constant, x is scaled to prevent overflow, but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
Similarly, a row-wise scheme is used to solve A**T*x = b. The basic algorithm for A upper triangular is
for j = 1, ..., n
x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
end
We simultaneously compute two bounds
G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
M(j) = bound on x(i), 1<=i<=j
The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. Then the bound on x(j) is
M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
<= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
1<=i<=j
and we can safely call STPSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).
subroutine slatrs (character UPLO, character TRANS, character DIAG, character NORMIN, integer N, real, dimension( lda, * ) A, integer LDA, real, dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow. ๐
Purpose:
SLATRS solves one of the triangular systems
A *x = s*b or A**T*x = s*b
with scaling to prevent overflow. Here A is an upper or lower triangular matrix, A**T denotes the transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine STRSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned. ๐ซถ
Parameters
UPLOโUPLOisCHARACTER*1. Specifies whether the matrixAis upper or lower triangular.= 'U': Upper triangular= 'L': Lower triangular
TRANSโTRANSisCHARACTER*1. Specifies the operation applied toA.= 'N': SolveA * x = s*b(No transpose)= 'T': SolveA**T* x = s*b(Transpose)= 'C': SolveA**T* x = s*b(Conjugate transpose = Transpose)
DIAGโDIAGisCHARACTER*1. Specifies whether or not the matrixAis unit triangular.= 'N': Non-unit triangular= 'U': Unit triangular
NORMINโNORMINisCHARACTER*1. Specifies whetherCNORMhas been set or not.= 'Y':CNORMcontains the column norms on entry= 'N':CNORMis not set on entry. On exit, the norms will be computed and stored inCNORM.
NโNisINTEGER. The order of the matrixA.N >= 0.AโAisREALarray, dimension(LDA,N). The triangular matrixA. IfUPLO = 'U', the leading n by n upper triangular part of the arrayAcontains the upper triangular matrix, and the strictly lower triangular part ofAis not referenced. IfUPLO = 'L', the leading n by n lower triangular part of the arrayAcontains the lower triangular matrix, and the strictly upper triangular part ofAis not referenced. IfDIAG = 'U', the diagonal elements ofAare also not referenced and are assumed to be1.LDAโLDAisINTEGER. The leading dimension of the arrayA.LDA >= max (1,N).XโXisREALarray, dimension(N). On entry, the right hand sidebof the triangular system. On exit,Xis overwritten by the solution vectorx.SCALEโSCALEisREAL. The scaling factorsfor the triangular system
A * x = s*b or A**T* x = s*b.
If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.
CNORMโCNORMisREALarray, dimension(N).
If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm.
If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.
INFOโINFOisINTEGER.= 0: successful exit< 0: ifINFO = -k, the k-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
A rough bound on x is computed; if that is less than overflow, STRSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation. ๐
A columnwise scheme is used for solving A*x = b. The basic algorithm if A is lower triangular is
x[1:n] := b[1:n]
for j = 1, ..., n
x(j) := x(j) / A(j,j)
x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
end
Define bounds on the components of x after j iterations of the loop:
M(j) = bound on x[1:j]
G(j) = bound on x[j+1:n]
Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
Then for iteration j+1 we have
M(j+1) <= G(j) / | A(j+1,j+1) |
G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
<= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
where CNORM(j+1) is greater than or equal to the infinity-norm of column j+1 of A, not counting the diagonal. Hence
G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
1<=i<=j
and
|x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
1<=i< j
Since |x(j)| <= M(j), we use the Level 2 BLAS routine STRSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).
The bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. If the computed bound is greater than a large constant, x is scaled to prevent overflow, but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
Similarly, a row-wise scheme is used to solve A**T*x = b. The basic algorithm for A upper triangular is
for j = 1, ..., n
x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
end
We simultaneously compute two bounds
G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
M(j) = bound on x(i), 1<=i<=j
The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. Then the bound on x(j) is
M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
<= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
1<=i<=j
and we can safely call STRSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).
subroutine slauu2 (character UPLO, integer N, real, dimension( lda, * ) A, integer LDA, integer INFO)
SLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm). โจ
Purpose:
SLAUU2 computes the product U * U**T or L**T * L, where the triangular factor U or L is stored in the upper or lower triangular part of the array A.
If UPLO = 'U' or 'u' then the upper triangle of the result is stored, overwriting the factor U in A.
If UPLO = 'L' or 'l' then the lower triangle of the result is stored, overwriting the factor L in A.
This is the unblocked form of the algorithm, calling Level 2 BLAS. ๐
Parameters
UPLOโUPLOisCHARACTER*1. Specifies whether the triangular factor stored in the arrayAis upper or lower triangular:= 'U': Upper triangular= 'L': Lower triangular
NโNisINTEGER. The order of the triangular factorUorL.N >= 0.AโAisREALarray, dimension(LDA,N). On entry, the triangular factorUorL. On exit, ifUPLO = 'U', the upper triangle ofAis overwritten with the upper triangle of the productU * U**T; ifUPLO = 'L', the lower triangle ofAis overwritten with the lower triangle of the productL**T * L.LDAโLDAisINTEGER. The leading dimension of the arrayA.LDA >= max(1,N).INFOINFOisINTEGERโจ= 0: successful exit, we love that for you ๐< 0: ifINFO = -k, the k-th argument had an illegal value ๐ญ (it's snitching on exactly which one too, iconic)
๐ซถ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐
slauum
subroutine slauum (character UPLO, integer N, real, dimension( lda, * ) A, integer LDA, integer INFO)
SLAUUM computes the product U*U**T or L**T*L, where U and L are upper or lower triangular matrices (blocked algorithm). bestie this is the blocked glow-up, not the baby version.
๐ฏ Purpose
SLAUUM computes the product U * U**T or L**T * L, where the triangular factor U or L is stored in the upper or lower triangular part of the array A. โจ
If UPLO = 'U' or 'u' then the upper triangle of the result is stored, overwriting the factor U in A. If UPLO = 'L' or 'l' then the lower triangle of the result is stored, overwriting the factor L in A.
This is the blocked form of the algorithm, calling Level 3 BLAS. that's the fast lane, no cap ๐ฅ
๐ Parameters
UPLOโUPLOisCHARACTER*1. It specifies whether the triangular factor stored in the arrayAis upper or lower triangular:= 'U': Upper triangular= 'L': Lower triangular
NโNisINTEGER. The order of the triangular factorUorL.N >= 0.AโAisREAL array, dimension (LDA,N). On entry, the triangular factorUorL. On exit, ifUPLO = 'U', the upper triangle ofAis overwritten with the upper triangle of the productU * U**T; ifUPLO = 'L', the lower triangle ofAis overwritten with the lower triangle of the productL**T * L. ๐ฆLDAโLDAisINTEGER. The leading dimension of the arrayA.LDA >= max(1,N).INFOโINFOisINTEGER= 0: successful exit< 0: ifINFO = -k, the k-th argument had an illegal value
๐ซถ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ธ srscl
subroutine srscl (integer N, real SA, real, dimension( * ) SX, integer INCX)
SRSCL multiplies a vector by the reciprocal of a real scalar. dividing by a without ever dividing by a, she's clever like that ๐
๐ฏ Purpose
SRSCL multiplies an n-element real vector x by the real scalar 1/a. This is done without overflow or underflow as long as the final result x/a does not overflow or underflow. ๐ซถ
๐ Parameters
NโNisINTEGER. The number of components of the vectorx.SAโSAisREAL. The scalarawhich is used to divide each component ofx.SAmust be>= 0, or the subroutine will divide by zero. ๐ญ genuinely, do not do thatSXโSXisREAL array, dimension (1+(N-1)*abs(INCX)). The n-element vectorx.INCXโINCXisINTEGER. The increment between successive values of the vectorSX.> 0:SX(1) = X(1)andSX(1+(i-1)*INCX) = x(i),1< i<= n
๐ซถ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ stprfb
subroutine stprfb (character SIDE, character TRANS, character DIRECT,
character STOREV, integer M, integer N, integer K, integer L, real,
dimension( ldv, * ) V, integer LDV, real, dimension( ldt, * ) T, integer
LDT, real, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * )
B, integer LDB, real, dimension( ldwork, * ) WORK, integer LDWORK)
STPRFB applies a real 'triangular-pentagonal' block reflector to a real matrix, which is composed of two blocks. two blocks, one reflector, main character energy only ๐
๐ฏ Purpose
STPRFB applies a real 'triangular-pentagonal' block reflector H or its transpose H**T to a real matrix C, which is composed of two blocks A and B, either from the left or right. โจ
๐ Parameters
SIDEโSIDEisCHARACTER*1= 'L': applyHorH**Tfrom the Left= 'R': applyHorH**Tfrom the Right
TRANSโTRANSisCHARACTER*1= 'N': applyH(No transpose)= 'T': applyH**T(Transpose)
DIRECTโDIRECTisCHARACTER*1. Indicates howHis formed from a product of elementary reflectors= 'F':H = H(1) H(2) . . . H(k)(Forward)= 'B':H = H(k) . . . H(2) H(1)(Backward)
STOREVโSTOREVisCHARACTER*1. Indicates how the vectors which define the elementary reflectors are stored:= 'C': Columns= 'R': Rows
MโMisINTEGER. The number of rows of the matrixB.M >= 0.NโNisINTEGER. The number of columns of the matrixB.N >= 0.KโKisINTEGER. The order of the matrixT, i.e. the number of elementary reflectors whose product defines the block reflector.K >= 0.LโLisINTEGER. The order of the trapezoidal part ofV.K >= L >= 0. See Further Details.VโVisREAL array, dimension
(LDV,K) if STOREV = 'C'
(LDV,M) if STOREV = 'R' and SIDE = 'L'
(LDV,N) if STOREV = 'R' and SIDE = 'R'
The pentagonal matrix V, which contains the elementary reflectors H(1), H(2), ..., H(K). See Further Details.
LDVโLDVisINTEGER. The leading dimension of the arrayV. IfSTOREV = 'C'andSIDE = 'L',LDV >= max(1,M); ifSTOREV = 'C'andSIDE = 'R',LDV >= max(1,N); ifSTOREV = 'R',LDV >= K.TโTisREAL array, dimension (LDT,K). The triangular K-by-K matrixTin the representation of the block reflector.LDTโLDTisINTEGER. The leading dimension of the arrayT.LDT >= K.AโAisREAL array, dimension
(LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
On entry, the K-by-N or M-by-K matrix A. On exit, A is overwritten by the corresponding block of H*C or H**T*C or C*H or C*H**T. See Further Details.
LDAโLDAisINTEGER. The leading dimension of the arrayA. IfSIDE = 'L',LDA >= max(1,K); IfSIDE = 'R',LDA >= max(1,M).BโBisREAL array, dimension (LDB,N). On entry, the M-by-N matrixB. On exit,Bis overwritten by the corresponding block ofH*CorH**T*CorC*HorC*H**T. See Further Details.LDBโLDBisINTEGER. The leading dimension of the arrayB.LDB >= max(1,M).WORKโWORKisREAL array, dimension
(LDWORK,N) if SIDE = 'L',
(LDWORK,K) if SIDE = 'R'.
LDWORKโLDWORKisINTEGER. The leading dimension of the arrayWORK. IfSIDE = 'L',LDWORK >= K; ifSIDE = 'R',LDWORK >= M.
๐ซถ Author
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
๐ Further Details:
The matrix C is a composite matrix formed from blocks A and B. The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, and if SIDE = 'L', A is of size K-by-N.
If SIDE = 'R' and DIRECT = 'F', C = [A B].
If SIDE = 'L' and DIRECT = 'F', C = [A]
[B].
If SIDE = 'R' and DIRECT = 'B', C = [B A].
If SIDE = 'L' and DIRECT = 'B', C = [B]
[A].
The pentagonal matrix V is composed of a rectangular block V1 and a trapezoidal block V2. The size of the trapezoidal block is determined by the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; if L=0, there is no trapezoidal block, thus V = V1 is rectangular. ๐ฆ
If DIRECT = 'F' and STOREV = 'C': V = [V1]
[V2]
- V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
- V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
If DIRECT = 'B' and STOREV = 'C': V = [V2]
[V1]
- V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
- V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
- If
STOREV = 'C'andSIDE = 'L',Vis M-by-K withV2L-by-K. - If
STOREV = 'C'andSIDE = 'R',Vis N-by-K withV2L-by-K. - If
STOREV = 'R'andSIDE = 'L',Vis K-by-M withV2K-by-L. - If
STOREV = 'R'andSIDE = 'R',Vis K-by-N withV2K-by-L.
๐ซถ Author
Generated automatically by Doxygen for LAPACK from the source code.
and that's a wrap bestie ๐ this was LAPACK Version 3.11.0, dated Sun Nov 27 2022, from the page realOTHERauxiliary(3) โ you made it all the way to the end of part 21, which is lowkey iconic of you. go drink some water, hydrate that main character brain ๐โจ
๐ฆ 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