PCLAHRD(l)

📦 version 1.5, in ubuntu jammy

💖 PCLAHRD(l) — the LAPACK auxiliary routine that eats subdiagonals for breakfast 💖

okay bestie, sit down, hydrate, because we're talking about PCLAHRD — a LAPACK auxiliary routine, version 1.5, and she is lowkey the backbone of your whole Hessenberg reduction era. no cap. 🦋

✨ NAME ✨

PCLAHRD — reduce the first NB columns of a complex 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. yes she really said all of that in one breath and she meant every word 💅

✨ SYNOPSIS ✨

SUBROUTINE PCLAHRD( 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( * )

    COMPLEX         A( * ), T( * ), TAU( * ), WORK( * ), Y( * )

like. look at that argument list. it's giving organized. it's giving parallel. 🎀

✨ PURPOSE ✨

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

the reduction gets done by a unitary similarity transformation Q' * A * Q. and then — because she's generous — the routine hands you back the matrices V and T that determine Q as a block reflector I - V*T*V', plus the matrix Y = A * V * T. three gifts. no receipt needed. 🫶

she's an auxiliary routine, meaning she's called by PCGEHRD — she's the supportive bestie in the group chat who does the real work behind the scenes and lets someone else take the credit. iconic behavior honestly.

and just so we're all on the same page: in the comments below, sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1). noted. memorized. moving on. 📝

✨ 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 ). and she has N >= 0 as a hard requirement, so don't come at her with negatives, she will simply not 💀
  • K — (global input) INTEGER. the offset for the reduction. the elements below the k-th subdiagonal in the first NB columns get reduced to zero. that's the whole assignment. that's the job. 🔥
  • NB — (global input) INTEGER. the number of columns to be reduced. how many columns are we cooking today, bestie? 🍓
  • A — (local input/local output) COMPLEX 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 — she does not touch what she was not asked to touch, we love a boundary-respecting queen 👑
  • IA — (global input) INTEGER. the row index in the global array A indicating the first row of sub( A ). (see Further Details, she has lore 📖)
  • 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 ID card. she's the vibe check. 🪪
  • TAU — (local output) COMPLEX array, dimension LOCc(JA+N-2). the scalar factors of the elementary reflectors (again: see Further Details). TAU is tied to the distributed matrix A — they're a package deal, inseparable, main character duo energy 💞
  • T — (local output) COMPLEX array, dimension (NB_A,NB_A). the upper triangular matrix T. triangle-shaped and thriving 🔺
  • Y — (local output) COMPLEX 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 heads up: LLD_Y >= LOCr(IA+N-1) — respect the leading dimension, she's not asking twice 🫠
  • 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. the other ID card 🪪
  • WORK — (local workspace) COMPLEX array, dimension (NB). just vibing. doing her little job. we stan a workspace 💜

✨ FURTHER DETAILS ✨

okay this is the part where it gets real pretty. 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 complex scalar, and v is a complex 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). 📌 memorize that indexing, bestie, it WILL be on the exam.

the elements of the vectors v together form the (n-k+1)-by-nb matrix V — which is needed, together 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').

like. that's a whole transformation. that's a whole glow-up. 💅✨

🌈 the illustration (yes, she drew you a picture) 🌈

the contents of A(ia:ia+n-1,ja:ja+n-k) on exit are illustrated by the following example 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 )

and now the legend, because we do not gatekeep knowledge in this house:

  • a — an element of the original matrix A(ia:ia+n-1,ja:ja+n-k) (untouched, unbothered, thriving 🌱)
  • h — a modified element of the upper Hessenberg matrix H (she got the makeover 💇‍♀️)
  • vi — an element of the vector defining H(i) (the reflection-vector girlies 🧚)
⚠️ real talk though: the reduction only overwrites what it says it overwrites, and the below-subdiagonal entries are where Q lives — so if you mess with them expecting them to still be your original data, that's on you, not on her. 🙈

that's the whole page, bestie. PCLAHRD, LAPACK auxiliary routine, version 1.5 — and yes, she's from 12 May 1997, a certified vintage icon still out here reducing subdiagonals with full main character energy 👑💖

📦 where this page comes from