ZLATRS(3)

💖 ZLATRS(3) — the triangular-solve-with-scaling glow-up 💖

okay so bestie, buckle in, because this one is a LAPACK routine from the zlatrs.f file and it's honestly kind of iconic. it lives in the Library Functions Manual corner of the world (the header up top even says /home/abui...RC/zlatrs.f(3)), and yes, the full path situation is giving very much /home/abuild/rpmbuild/BUILD/lapack-static-3.12.0-build/lapack-3.12.0/SRC/zlatrs.f. don't worry about it. 🫶

---

✨ SYNOPSIS ✨

Functions/Subroutines

subroutine ZLATRS (uplo, trans, diag, normin, n, a, lda, x, scale,
    cnorm, info)

ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow. no cap, that's literally the whole gig. 💅

---

🦋 Function/Subroutine Documentation 🦋

subroutine ZLATRS (character uplo, character trans, character diag,
    character normin, integer n, complex*16, dimension( lda, * ) a, integer
    lda, complex*16, dimension( * ) x, double precision scale, double
    precision, dimension( * ) cnorm, integer info)

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

---

🎀 Purpose 🎀

ok so here's the tea, straight from the source:

fortran ZLATRS solves one of the triangular systems A * x = s*b, A**T * x = s*b, or A**H * x = s*b, with scaling to prevent overflow. Here A is an upper or lower triangular matrix, A**T denotes the transpose of A, A**H denotes the conjugate transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine ZTRSV 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.

like, it's out here solving A * x = s*b variants and keeping everything from yeeting itself into overflow territory. main character energy fr fr. 👑

---

💜 Parameters 💜

UPLO

fortran UPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangular

• 'U' — upper triangular, the classic ✨
• 'L' — lower triangular, also valid, we don't judge 💅

TRANS

fortran TRANS is CHARACTER*1 Specifies the operation applied to A. = 'N': Solve A * x = s*b (No transpose) = 'T': Solve A**T * x = s*b (Transpose) = 'C': Solve A**H * x = s*b (Conjugate transpose)

• 'N' — no transpose, just vibes
• 'T' — transpose it bestie
• 'C' — conjugate transpose, the fancy one 🦋

DIAG

fortran DIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangular

• 'N' — non-unit triangular
• 'U' — unit triangular (diagonal is assumed to be 1s)

NORMIN

fortran NORMIN is CHARACTER*1 Specifies whether CNORM has been set or not. = 'Y': CNORM contains the column norms on entry = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.

so 'Y' means the norms are already fed in, and 'N' means this routine is gonna compute them for you and stash them in CNORM on the way out. teamwork!! 🫶

N

fortran N is INTEGER The order of the matrix A. N >= 0.

the size of the whole situation. and yes, N >= 0, we're being strict about it.

A

fortran A is COMPLEX*16 array, dimension (LDA,N) The triangular matrix A. If UPLO = 'U', the leading n by n upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1.

so depending on UPLO, one triangle of A is the moment and the other is literally not referenced. and if DIAG = 'U', the diagonal is skipped too and assumed to be 1. respectful of your time honestly. ✨

LDA

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

the leading dimension bestie. LDA >= max (1,N). don't fumble it.

X

fortran X is COMPLEX*16 array, dimension (N) On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.

comes in as the right hand side b, leaves as the solution x. glow-up of the century, honestly. 😭

SCALE

fortran SCALE is DOUBLE PRECISION The scaling factor s for the triangular system A * x = s*b, A**T * x = s*b, or A**H * x = s*b. If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.

the scaling factor s. and heads up: if SCALE = 0, that means A is singular or badly scaled and x is an exact or approximate solution to A*x = 0. not a drill. 💀

CNORM

fortran CNORM is DOUBLE PRECISION array, dimension (N) If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm. If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.

okay so it flips roles depending on NORMIN:

• NORMIN = 'Y' → it's an input, CNORM(j) holds the norm of the off-diagonal part of the j-th column of A. if TRANS = 'N' it must be >= the infinity-norm; if TRANS = 'T' or 'C' it must be >= the 1-norm.
• NORMIN = 'N' → it's an output, and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.

INFO

fortran INFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal value

• 0 — slayed, successful exit ✨
• < 0 — if INFO = -k, the k-th argument had an illegal value. so like, check your args bestie 🙈

---

🌸 Author 🌸

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

the whole squad contributed. iconic collaboration. 🫶

---

🌈 Further Details: 🌈

this is the part where it gets lowkey math-gorgeous, so stay with me:

fortran A rough bound on x is computed; if that is less than overflow, ZTRSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation.

so it peeks at a rough bound first — cheap check — and if that's safe, it just calls ZTRSV. if not, it goes into careful mode and checks for overflow or divide-by-zero at every single operation. thorough queen. 👑

fortran 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

columnwise scheme for A*x = b, and that's the lower-triangular basic algorithm right there. ⭐

fortran 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}.

the two bound besties: M(j) for x[1:j], G(j) for x[j+1:n], starting at M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

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

fortran 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

so that's the whole bound chain, and here's the payoff:

fortran Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTRSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow).

i.e. if the reciprocal of the largest M(j) across j=1,..,n beats max(underflow, 1/overflow), ZTRSV gets the call. 💅

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

so if the bound is bigger than a large constant → scale x to prevent overflow. but if the bound itself overflows → x becomes 0, x(j) becomes 1, scale becomes 0, and boom, you've got a non-trivial solution to A*x = 0. it's giving emergency rescue mode. 😭

fortran Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *x = b. The basic algorithm for A upper triangular is for j = 1, ..., n x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) end

row-wise scheme for the transpose/conjugate-transpose cases. same energy, different direction. 🦋

fortran We simultaneously compute two bounds G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j M(j) = bound on x(i), 1<=i<=j The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. Then the bound on x(j) is M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 1<=i<=j and we can safely call ZTRSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).

two bounds running in parallel with the constraints G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1, and the green light for ZTRSV is when both 1/M(n) and 1/G(n) are greater than max(underflow, 1/overflow). periodt. 💅

---

🧚 Definition location 🧚

Definition at line 237 of file zlatrs.f.

---

⭐ Author ⭐

Generated automatically by Doxygen for LAPACK from the source code.

---

and that's the whole page bestie! 🎀 made with love by the folks at LAPACK, sitting pretty at Version 3.12.0, and dated right alongside /home/abui...RC/zlatrs.f(3) — mwah, go solve some triangular systems responsibly 💖✨