PDLAHRD(l)

📦 version 1.5, in ubuntu jammy

💖 PDLAHRD — the LAPACK auxiliary routine that's lowkey carrying the whole reduction 💖

okay bestie, gather round, because today we're talking about a LAPACK auxiliary routine from way back (version 1.5, we love a vintage icon) and honestly? she's a supporting character with main character energy. no cap. 🎀

✨ NAME ✨

PDLAHRD reduces the first NB columns of a real general N-by-(N-K+1) distributed matrix A(IA:IA+N-1,JA:JA+N-K) so that elements below the k-th subdiagonal are zero. that's the whole gig. she zeroes the vibes below the subdiagonal and she does it flawlessly 💅

🦋 SYNOPSIS 🦋

SUBROUTINE PDLAHRD( N, K, NB, A, IA, JA, DESCA, TAU, T, Y, IY, JY,
                    DESCY, WORK )

    INTEGER         IA, IY, JA, JY, K, N, NB

    INTEGER         DESCA( * ), DESCY( * )

    DOUBLE          PRECISION A( * ), T( * ), TAU( * ), WORK( * ), Y( *
                    )

🌸 PURPOSE 🌸

so here's the tea ☕ — PDLAHRD reduces the first NB columns of a real general N-by-(N-K+1) distributed matrix A(IA:IA+N-1,JA:JA+N-K) so that elements below the k-th subdiagonal are zero. periodt.

the reduction is performed by an orthogonal similarity transformation Q' * A * Q. and then the routine hands you back the matrices V and T which determine Q as a block reflector I - V*T*V', plus the matrix Y = A * V * T. iconic behavior honestly 🫶

she's an auxiliary routine called by PDGEHRD. not the star of the show, but the show literally does not happen without her. in the following comments, sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1). iykyk ✨

👑 ARGUMENTS 👑

  • N (global input) INTEGER — the number of rows and columns to be operated on, i.e. the order of the distributed submatrix sub( A ). bestie, N >= 0. non-negotiable. 🌈
  • K (global input) INTEGER — the offset for the reduction. elements below the k-th subdiagonal in the first NB columns are reduced to zero. ✨
  • NB (global input) INTEGER — the number of columns to be reduced. that's how many we're doing. 💅
  • A (local input/local output) DOUBLE PRECISION pointer into the local memory to an array of dimension (LLD_A, LOCc(JA+N-K)). on entry, this array contains the local pieces of the N-by-(N-K+1) general distributed matrix A(IA:IA+N-1,JA:JA+N-K). on exit, the elements on and above the k-th subdiagonal in the first NB columns are overwritten with the corresponding elements of the reduced distributed matrix; the elements below the k-th subdiagonal, together with the array TAU, represent the matrix Q as a product of elementary reflectors. the other columns of A(IA:IA+N-1,JA:JA+N-K) are unchanged. see Further Details. 🦋
  • IA (global input) INTEGER — the row index in the global array A indicating the first row of sub( A ). 📍
  • JA (global input) INTEGER — the column index in the global array A indicating the first column of sub( A ). 📍
  • DESCA (global and local input) INTEGER array of dimension DLEN_. the array descriptor for the distributed matrix A. she's the map, the matrix is the territory 🗺️
  • TAU (local output) DOUBLE PRECISION array, dimension LOCc(JA+N-2). the scalar factors of the elementary reflectors (see Further Details). TAU is tied to the distributed matrix A. they're bonded, fr fr 🫶
  • T (local output) DOUBLE PRECISION array, dimension (NB_A,NB_A). the upper triangular matrix T. ⭐
  • Y (local output) DOUBLE PRECISION pointer into the local memory to an array of dimension (LLD_Y,NB_A). on exit, this array contains the local pieces of the N-by-NB distributed matrix Y. and listen: LLD_Y >= LOCr(IA+N-1). no cap. 💖
  • IY (global input) INTEGER — the row index in the global array Y indicating the first row of sub( Y ). 📍
  • JY (global input) INTEGER — the column index in the global array Y indicating the first column of sub( Y ). 📍
  • DESCY (global and local input) INTEGER array of dimension DLEN_. the array descriptor for the distributed matrix Y. 🗺️
  • WORK (local workspace) DOUBLE PRECISION array, dimension (NB). just the little scratchpad. do not perceive her 💀

🧚 FURTHER DETAILS 🧚

okay this is the part where it gets good. 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'

where tau is a real scalar, and v is a real vector with

v(1:i+k-1) = 0, v(i+k) = 1;

and v(i+k+1:n) is stored on exit in A(ia+i+k:ia+n-1,ja+i-1), and tau in TAU(ja+i-1). yes, all of that. every index. she's precise like that 💅

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(ia:ia+n-1,ja:ja+n-k) := (I-V*T*V')*(A(ia:ia+n-1,ja:ja+n-k)-Y*V').

and because visuals eat every time, here's what the contents of A(ia:ia+n-1,ja:ja+n-k) look like on exit, illustrated with n = 7, k = 3 and nb = 2:

( a   h   a   a   a )
( a   h   a   a   a )
( a   h   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(ia:ia+n-1,ja:ja+n-k), h denotes a modified element of the upper Hessenberg matrix H, and vi denotes an element of the vector defining H(i). 🎀

ok real talk bestie, this is a distributed-memory auxiliary routine — the descriptors DESCA and DESCY and the index arguments IA, JA, IY, JY are not decoration, they're the whole point. get them wrong and nothing means anything. 🙈

that's the page!! go forth and reduce some columns, you absolute icon 👑

this has been PDLAHRD from LAPACK version 1.5, dated 12 May 1997 — yes she's older than you and she's still slaying, we love a timeless queen 💖✨

📦 where this page comes from