unhr_col(3)

🌸 unhr_col(3) — Library Functions Manual 🌸

okay bestie, welcome to the unhr_col page — aka the {un,or}hr_col house of Householder reconstruction. This is the Library Functions Manual entry, and it is a lot, so buckle up. We're doing this in one long glorious scroll. 🎀

---

✨ NAME ✨

unhr_col — {un,or}hr_col: Householder reconstruction

---

✨ SYNOPSIS ✨

🦋 Functions

• subroutine cunhr_col (m, n, nb, a, lda, t, ldt, d, info) — CUNHR_COL
• subroutine dorhr_col (m, n, nb, a, lda, t, ldt, d, info) — DORHR_COL
• subroutine sorhr_col (m, n, nb, a, lda, t, ldt, d, info) — SORHR_COL
• subroutine zunhr_col (m, n, nb, a, lda, t, ldt, d, info) — ZUNHR_COL

---

✨ DETAILED DESCRIPTION ✨

(no extra prose here babes, the details live down in the function docs ⬇️)

---

✨ FUNCTION DOCUMENTATION ✨

💜 subroutine cunhr_col (integer m, integer n, integer nb, complex, dimension( lda, * ) a, integer lda, complex, dimension( ldt, * ) t, integer ldt, complex, dimension( * ) d, integer info)

🌟 CUNHR_COL 🌟

Purpose:

!>
!>  CUNHR_COL takes an M-by-N complex matrix Q_in with orthonormal columns
!>  as input, stored in A, and performs Householder Reconstruction (HR),
!>  i.e. reconstructs Householder vectors V(i) implicitly representing
!>  another M-by-N matrix Q_out, with the property that Q_in = Q_out*S,
!>  where S is an N-by-N diagonal matrix with diagonal entries
!>  equal to +1 or -1. The Householder vectors (columns V(i) of V) are
!>  stored in A on output, and the diagonal entries of S are stored in D.
!>  Block reflectors are also returned in T
!>  (same output format as CGEQRT).
!>

so like, imagine you hand it an M-by-N complex matrix Q_in with orthonormal columns, chillin' in A — she takes that and does Householder Reconstruction (HR, iconic acronym). Basically she rebuilds the Householder vectors V(i) that implicitly represent a whole other M-by-N matrix Q_out, and the whole point is:

Q_in = Q_out*S

where S is an N-by-N diagonal matrix whose diagonal entries are literally just +1 or -1. No in-between, no vibes. 💅 On output the Householder vectors (the columns V(i) of V) live in A, and the diagonal entries of S live in D. You also get block reflectors back in T — same output format as CGEQRT, no cap.

Parameters

• M

!>          M is INTEGER
!>          The number of rows of the matrix A. M >= 0.
!>

that's the row count bestie. M >= 0.

• N

!>          N is INTEGER
!>          The number of columns of the matrix A. M >= N >= 0.
!>

column count. note the constraint is M >= N >= 0 — N can't clown around and be bigger than M. 🙈

• NB

!>          NB is INTEGER
!>          The column block size to be used in the reconstruction
!>          of Householder column vector blocks in the array A and
!>          corresponding block reflectors in the array T. NB >= 1.
!>          (Note that if NB > N, then N is used instead of NB
!>          as the column block size.)
!>

the column block size used when reconstructing the Householder column vector blocks in A and the matching block reflectors in T. NB >= 1. ⚠️ and heads up: if NB > N, then N gets used instead of NB as the column block size — the code just quietly clamps, iykyk.

• A

!>          A is COMPLEX array, dimension (LDA,N)
!>
!>          On entry:
!>
!>             The array A contains an M-by-N orthonormal matrix Q_in,
!>             i.e the columns of A are orthogonal unit vectors.
!>
!>          On exit:
!>
!>             The elements below the diagonal of A represent the unit
!>             lower-trapezoidal matrix V of Householder column vectors
!>             V(i). The unit diagonal entries of V are not stored
!>             (same format as the output below the diagonal in A from
!>             CGEQRT). The matrix T and the matrix V stored on output
!>             in A implicitly define Q_out.
!>
!>             The elements above the diagonal contain the factor U
!>             of the  LU-decomposition:
!>                Q_in - ( S ) = V * U
!>                       ( 0 )
!>             where 0 is a (M-N)-by-(M-N) zero matrix.
!>

On entry: A holds your M-by-N orthonormal matrix Q_in — meaning the columns of A are orthogonal unit vectors. 🦋

On exit: everything below the diagonal of A is the unit lower-trapezoidal matrix V made of Householder column vectors V(i). The unit diagonal entries of V are not stored (same format as the output below the diagonal in A from CGEQRT). Together, the matrix T and the matrix V sitting in A on output implicitly define Q_out. ✨

And the elements above the diagonal hold the factor U of the LU-decomposition:

   Q_in - ( S ) = V * U
          ( 0 )

where 0 is a (M-N)-by-(M-N) zero matrix.

• LDA

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

leading dimension of A. LDA >= max(1,M) — respect the bounds or it's a bad time. 💀

• T

!>          T is COMPLEX array,
!>          dimension (LDT, N)
!>
!>          Let NOCB = Number_of_output_col_blocks
!>                   = CEIL(N/NB)
!>
!>          On exit, T(1:NB, 1:N) contains NOCB upper-triangular
!>          block reflectors used to define Q_out stored in compact
!>          form as a sequence of upper-triangular NB-by-NB column
!>          blocks (same format as the output T in CGEQRT).
!>          The matrix T and the matrix V stored on output in A
!>          implicitly define Q_out. NOTE: The lower triangles
!>          below the upper-triangular blocks will be filled with
!>          zeros. See Further Details.
!>

T is a COMPLEX array of dimension (LDT, N). Here's the setup, bestie:

   Let NOCB = Number_of_output_col_blocks
            = CEIL(N/NB)

On exit, T(1:NB, 1:N) contains NOCB upper-triangular block reflectors that define Q_out, stored compactly as a sequence of upper-triangular NB-by-NB column blocks (same format as the output T in CGEQRT). The matrix T plus the matrix V stored in A on output implicitly define Q_out.

🚨 NOTE: the lower triangles below those upper-triangular blocks get filled with zeros. See Further Details, babe.

• LDT

!>          LDT is INTEGER
!>          The leading dimension of the array T.
!>          LDT >= max(1,min(NB,N)).
!>

leading dimension of T, and it's gotta be LDT >= max(1,min(NB,N)). 🫶

• D

!>          D is COMPLEX array, dimension min(M,N).
!>          The elements can be only plus or minus one.
!>
!>          D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where
!>          1 <= i <= min(M,N), and Q_in_i is Q_in after performing
!>          i-1 steps of “modified” Gaussian elimination.
!>          See Further Details.
!>

D is COMPLEX with dimension min(M,N). The elements can only be plus or minus one — that's the whole vibe check, no other values allowed. 💅

And it's built like this:

   D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where
   1 <= i <= min(M,N), and Q_in_i is Q_in after performing
   i-1 steps of “modified” Gaussian elimination.

(the "modified" has little curly quotes because she's fancy like that 🎀) — See Further Details.

• INFO

!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>

= 0 means successful exit, slay. < 0 means if INFO = -i, the i-th argument had an illegal value — go check your args. 😭

Further Details:

!>
!> The computed M-by-M unitary factor Q_out is defined implicitly as
!> a product of unitary matrices Q_out(i). Each Q_out(i) is stored in
!> the compact WY-representation format in the corresponding blocks of
!> matrices V (stored in A) and T.
!>
!> The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N
!> matrix A contains the column vectors V(i) in NB-size column
!> blocks VB(j). For example, VB(1) contains the columns
!> V(1), V(2), ... V(NB). NOTE: The unit entries on
!> the diagonal of Y are not stored in A.
!>
!> The number of column blocks is
!>
!>     NOCB = Number_of_output_col_blocks = CEIL(N/NB)
!>
!> where each block is of order NB except for the last block, which
!> is of order LAST_NB = N - (NOCB-1)*NB.
!>
!> For example, if M=6,  N=5 and NB=2, the matrix V is
!>
!>
!>     V = (    VB(1),   VB(2), VB(3) ) =
!>
!>       = (   1                      )
!>         ( v21    1                 )
!>         ( v31  v32    1            )
!>         ( v41  v42  v43   1        )
!>         ( v51  v52  v53  v54    1  )
!>         ( v61  v62  v63  v54   v65 )
!>
!>
!> For each of the column blocks VB(i), an upper-triangular block
!> reflector TB(i) is computed. These blocks are stored as
!> a sequence of upper-triangular column blocks in the NB-by-N
!> matrix T. The size of each TB(i) block is NB-by-NB, except
!> for the last block, whose size is LAST_NB-by-LAST_NB.
!>
!> For example, if M=6,  N=5 and NB=2, the matrix T is
!>
!>     T  = (    TB(1),    TB(2), TB(3) ) =
!>
!>        = ( t11  t12  t13  t14   t15  )
!>          (      t22       t24        )
!>
!>
!> The M-by-M factor Q_out is given as a product of NOCB
!> unitary M-by-M matrices Q_out(i).
!>
!>     Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),
!>
!> where each matrix Q_out(i) is given by the WY-representation
!> using corresponding blocks from the matrices V and T:
!>
!>     Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,
!>
!> where I is the identity matrix. Here is the formula with matrix
!> dimensions:
!>
!>  Q(i){M-by-M} = I{M-by-M} -
!>    VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},
!>
!> where INB = NB, except for the last block NOCB
!> for which INB=LAST_NB.
!>
!> =====
!> NOTE:
!> =====
!>
!> If Q_in is the result of doing a QR factorization
!> B = Q_in * R_in, then:
!>
!> B = (Q_out*S) * R_in = Q_out * (S * R_in) = Q_out * R_out.
!>
!> So if one wants to interpret Q_out as the result
!> of the QR factorization of B, then the corresponding R_out
!> should be equal to R_out = S * R_in, i.e. some rows of R_in
!> should be multiplied by -1.
!>
!> For the details of the algorithm, see [1].
!>
!> [1] ,
!>     G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
!>     E. Solomonik, J. Parallel Distrib. Comput.,
!>     vol. 85, pp. 3-31, 2015.
!>

okay so the tea ☕: the computed M-by-M unitary factor Q_out is defined implicitly as a product of unitary matrices Q_out(i). Each Q_out(i) is stored in compact WY-representation format in the corresponding blocks of the matrices V (stored in A) and T.

The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N matrix A contains the column vectors V(i) in NB-size column blocks VB(j). Like, VB(1) contains the columns V(1), V(2), ... V(NB). ⚠️ NOTE: the unit entries on the diagonal of Y are not stored in A.

Number of column blocks:

    NOCB = Number_of_output_col_blocks = CEIL(N/NB)

where every block is of order NB except the last one, which is of order LAST_NB = N - (NOCB-1)*NB.

Example time, bestie — if M=6, N=5 and NB=2, the matrix V is:

    V = (    VB(1),   VB(2), VB(3) ) =

      = (   1                      )
        ( v21    1                 )
        ( v31  v32    1            )
        ( v41  v42  v43   1        )
        ( v51  v52  v53  v54    1  )
        ( v61  v62  v63  v54   v65 )

For each of the column blocks VB(i), an upper-triangular block reflector TB(i) is computed. Those get stored as a sequence of upper-triangular column blocks in the NB-by-N matrix T. Each TB(i) block is NB-by-NB, except the last block, whose size is LAST_NB-by-LAST_NB.

Same example (M=6, N=5, NB=2), now for T:

    T  = (    TB(1),    TB(2), TB(3) ) =

       = ( t11  t12  t13  t14   t15  )
         (      t22       t24        )

The M-by-M factor Q_out is the product of NOCB unitary M-by-M matrices Q_out(i):

    Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),

where each matrix Q_out(i) comes from the WY-representation using the corresponding blocks from V and T:

    Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,

and I is the identity matrix, obviously. 💅 Here's the same formula with matrix dimensions spelled out:

 Q(i){M-by-M} = I{M-by-M} -
   VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},

where INB = NB, except for the last block NOCB where INB=LAST_NB.

⭐ THE NOTE SECTION — pay attention bestie ⭐ If Q_in is the result of doing a QR factorization B = Q_in * R_in, then: text B = (Q_out*S) * R_in = Q_out * (S * R_in) = Q_out * R_out. So if you wanna interpret Q_out as the result of the QR factorization of B, then the corresponding R_out should be equal to R_out = S * R_in — i.e. some rows of R_in should be multiplied by -1. 🙈

For the details of the algorithm, see [1]:

[1] ,
    G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
    E. Solomonik, J. Parallel Distrib. Comput.,
    vol. 85, pp. 3-31, 2015.

Author

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

Contributors:

!>
!> November   2019, Igor Kozachenko,
!>            Computer Science Division,
!>            University of California, Berkeley
!>
!>

Definition at line 258 of file cunhr_col.f.

---

🦋 subroutine dorhr_col (integer m, integer n, integer nb, double precision, dimension( lda, * ) a, integer lda, double precision, dimension( ldt, * ) t, integer ldt, double precision, dimension( * ) d, integer info)

🌟 DORHR_COL 🌟

Purpose:

!>
!>  DORHR_COL takes an M-by-N real matrix Q_in with orthonormal columns
!>  as input, stored in A, and performs Householder Reconstruction (HR),
!>  i.e. reconstructs Householder vectors V(i) implicitly representing
!>  another M-by-N matrix Q_out, with the property that Q_in = Q_out*S,
!>  where S is an N-by-N diagonal matrix with diagonal entries
!>  equal to +1 or -1. The Householder vectors (columns V(i) of V) are
!>  stored in A on output, and the diagonal entries of S are stored in D.
!>  Block reflectors are also returned in T
!>  (same output format as DGEQRT).
!>

same energy as her complex sister up there, except this one's for a real M-by-N matrix. 💚 DORHR_COL takes an M-by-N real matrix Q_in with orthonormal columns as input, stored in A, and performs Householder Reconstruction (HR) — reconstructing Householder vectors V(i) that implicitly represent another M-by-N matrix Q_out, with the property that

Q_in = Q_out*S

where S is an N-by-N diagonal matrix with diagonal entries equal to +1 or -1. On output the Householder vectors (columns V(i) of V) are stored in A, the diagonal entries of S are stored in D, and block reflectors are also returned in T — same output format as DGEQRT. 🫶

Parameters

• M

!>          M is INTEGER
!>          The number of rows of the matrix A. M >= 0.
!>

row count of A, M >= 0.

• N

!>          N is INTEGER
!>          The number of columns of the matrix A. M >= N >= 0.
!>

column count of A, and again M >= N >= 0. ✨

• NB

!>          NB is INTEGER
!>          The column block size to be used in the reconstruction
!>          of Householder column vector blocks in the array A and
!>          corresponding block reflectors in the array T. NB >= 1.
!>          (Note that if NB > N, then N is used instead of NB
!>          as the column block size.)
!>

column block size for reconstructing the Householder column vector blocks in A and the corresponding block reflectors in T. NB >= 1, and yes — if NB > N, then N is used instead of NB as the column block size. same clamp, different subroutine. 💅

• A —
• A — DOUBLE PRECISION array, dimension (LDA,N) 💖

On entry: the array A contains an M-by-N orthonormal matrix Q_in, i.e. the columns of A are orthogonal unit vectors. ✨

On exit: the elements below the diagonal of A represent the unit lower-trapezoidal matrix V of Householder column vectors V(i). The unit diagonal entries of V are not stored (same format as the output below the diagonal in A from DGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out. 🫶

The elements above the diagonal contain the factor U of the LU-decomposition:

   Q_in - ( S ) = V * U
          ( 0 )

where 0 is a (M-N)-by-(M-N) zero matrix. no cap, that zero block is literally just zeros 💀

• LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). 📏
• T — DOUBLE PRECISION array, dimension (LDT, N) 🎀

Let NOCB = Number_of_output_col_blocks = CEIL(N/NB).

On exit, T(1:NB, 1:N) contains NOCB upper-triangular block reflectors used to define Q_out, stored in compact form as a sequence of upper-triangular NB-by-NB column blocks (same format as the output T in DGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out.

NOTE: the lower triangles below the upper-triangular blocks will be filled with zeros. See Further Details. 👀

• LDT — INTEGER. The leading dimension of the array T. LDT >= max(1,min(NB,N)). 💅
• D — DOUBLE PRECISION array, dimension min(M,N). The elements can be only plus or minus one. ⭐

D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where 1 <= i <= min(M,N), and Q_in_i is Q_in after performing i-1 steps of “modified” Gaussian elimination. See Further Details. 🧚

• INFO — INTEGER. Exit status, bestie:
  • = 0: successful exit 🌈
  • < 0: if INFO = -i, the i-th argument had an illegal value 😭

---

✨ Further Details ✨

The computed M-by-M orthogonal factor Q_out is defined implicitly as a product of orthogonal matrices Q_out(i). Each Q_out(i) is stored in the compact WY-representation format in the corresponding blocks of matrices V (stored in A) and T. 🔥

The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N matrix A contains the column vectors V(i) in NB-size column blocks VB(j). For example, VB(1) contains the columns V(1), V(2), ... V(NB).

NOTE: the unit entries on the diagonal of Y are not stored in A. (yes, Y, that's what it says — iykyk 🫶)

The number of column blocks is

    NOCB = Number_of_output_col_blocks = CEIL(N/NB)

where each block is of order NB except for the last block, which is of order LAST_NB = N - (NOCB-1)*NB.

For example, if M=6, N=5 and NB=2, the matrix V is

    V = (    VB(1),   VB(2), VB(3) ) =

      = (   1                      )
        ( v21    1                 )
        ( v31  v32    1            )
        ( v41  v42  v43   1        )
        ( v51  v52  v53  v54    1  )
        ( v61  v62  v63  v54   v65 )

For each of the column blocks VB(i), an upper-triangular block reflector TB(i) is computed. These blocks are stored as a sequence of upper-triangular column blocks in the NB-by-N matrix T. The size of each TB(i) block is NB-by-NB, except for the last block, whose size is LAST_NB-by-LAST_NB.

For example, if M=6, N=5 and NB=2, the matrix T is

    T  = (    TB(1),    TB(2), TB(3) ) =

       = ( t11  t12  t13  t14   t15  )
         (      t22       t24        )

The M-by-M factor Q_out is given as a product of NOCB orthogonal M-by-M matrices Q_out(i).

    Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),

where each matrix Q_out(i) is given by the WY-representation using corresponding blocks from the matrices V and T:

    Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,

where I is the identity matrix. Here is the formula with matrix dimensions:

 Q(i){M-by-M} = I{M-by-M} -
   VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},

where INB = NB, except for the last block NOCB for which INB=LAST_NB.

=====
NOTE:
=====

If Q_in is the result of doing a QR factorization B = Q_in * R_in, then:

B = (Q_out*S) * R_in = Q_out * (S * R_in) = Q_out * R_out.

So if one wants to interpret Q_out as the result of the QR factorization of B, then the corresponding R_out should be equal to R_out = S * R_in, i.e. some rows of R_in should be multiplied by -1. 😮‍💨

For the details of the algorithm, see [1].

[1] , G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen, E. Solomonik, J. Parallel Distrib. Comput., vol. 85, pp. 3-31, 2015. 📚

---

👑 Author 👑

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

---

🫶 Contributors: 🫶

November 2019, Igor Kozachenko, Computer Science Division, University of California, Berkeley

---

Definition at line 258 of file dorhr_col.f.

---

subroutine sorhr_col (integer m, integer n, integer nb, real, dimension(
    lda, * ) a, integer lda, real, dimension( ldt, * ) t, integer ldt, real,
    dimension( * ) d, integer info)

SORHR_COL 🩵

Purpose

SORHR_COL takes an M-by-N real matrix Q_in with orthonormal columns as input, stored in A, and performs Householder Reconstruction (HR), i.e. reconstructs Householder vectors V(i) implicitly representing another M-by-N matrix Q_out, with the property that Q_in = Q_out*S, where S is an N-by-N diagonal matrix with diagonal entries equal to +1 or -1. The Householder vectors (columns V(i) of V) are stored in A on output, and the diagonal entries of S are stored in D. Block reflectors are also returned in T (same output format as SGEQRT). ✨

Parameters

• M — INTEGER. The number of rows of the matrix A. M >= 0. 📐
• N — INTEGER. The number of columns of the matrix A. M >= N >= 0. 📏
• NB — INTEGER. The column block size to be used in the reconstruction of Householder column vector blocks in the array A and corresponding block reflectors in the array T. NB >= 1.

(Note that if NB > N, then N is used instead of NB as the column block size.) 🎀

• A — REAL array, dimension (LDA,N) 💖

On entry: the array A contains an M-by-N orthonormal matrix Q_in, i.e. the columns of A are orthogonal unit vectors. ✨

On exit: the elements below the diagonal of A represent the unit lower-trapezoidal matrix V of Householder column vectors V(i). The unit diagonal entries of V are not stored (same format as the output below the diagonal in A from SGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out. 🫶

The elements above the diagonal contain the factor U of the LU-decomposition:

   Q_in - ( S ) = V * U
          ( 0 )

where 0 is a (M-N)-by-(M-N) zero matrix. literally all zeros, bestie 💀

• LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). 📏
• T — REAL array, dimension (LDT, N) 🎀

Let NOCB = Number_of_output_col_blocks = CEIL(N/NB).

On exit, T(1:NB, 1:N) contains NOCB upper-triangular block reflectors used to define Q_out, stored in compact form as a sequence of upper-triangular NB-by-NB column blocks (same format as the output T in SGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out.

NOTE: the lower triangles below the upper-triangular blocks will be filled with zeros. See Further Details. 👀

• LDT — INTEGER. The leading dimension of the array T. LDT >= max(1,min(NB,N)). 💅
• D — REAL array, dimension min(M,N). The elements can be only plus or minus one. ⭐

D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where 1 <= i <= min(M,N), and Q_in_i is Q_in after performing i-1 steps of “modified” Gaussian elimination. See Further Details. 🧚

• INFO — INTEGER. Exit status, bestie:
  • = 0: successful exit 🌈
  • < 0: if INFO = -i, the i-th argument had an illegal value 😭

✨ Further Details ✨

okay so the computed M-by-M orthogonal factor Q_out is defined implicitly as a product of orthogonal matrices Q_out(i). Each Q_out(i) is stored in the compact WY-representation format in the corresponding blocks of matrices V (stored in A) and T. 💅

The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N matrix A contains the column vectors V(i) in NB-size column blocks VB(j). For example, VB(1) contains the columns V(1), V(2), ... V(NB). NOTE: the unit entries on the diagonal of Y are not stored in A (yes bestie, it says Y here — we keep it moving). 🫶

The number of column blocks is

NOCB = Number_of_output_col_blocks = CEIL(N/NB)

where each block is of order NB except for the last block, which is of order LAST_NB = N - (NOCB-1)*NB.

For example, if M=6, N=5 and NB=2, the matrix V is

    V = (    VB(1),   VB(2), VB(3) ) =

      = (   1                      )
        ( v21    1                 )
        ( v31  v32    1            )
        ( v41  v42  v43   1        )
        ( v51  v52  v53  v54    1  )
        ( v61  v62  v63  v54   v65 )

For each of the column blocks VB(i), an upper-triangular block reflector TB(i) is computed. These blocks are stored as a sequence of upper-triangular column blocks in the NB-by-N matrix T. The size of each TB(i) block is NB-by-NB, except for the last block, whose size is LAST_NB-by-LAST_NB.

For example, if M=6, N=5 and NB=2, the matrix T is

    T  = (    TB(1),    TB(2), TB(3) ) =

       = ( t11  t12  t13  t14   t15  )
         (      t22       t24        )

The M-by-M factor Q_out is given as a product of NOCB orthogonal M-by-M matrices Q_out(i).

    Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),

where each matrix Q_out(i) is given by the WY-representation using corresponding blocks from the matrices V and T:

    Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,

where I is the identity matrix. Here is the formula with matrix dimensions:

 Q(i){M-by-M} = I{M-by-M} -
   VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},

where INB = NB, except for the last block NOCB for which INB=LAST_NB.

💀 NOTE

If Q_in is the result of doing a QR factorization B = Q_in * R_in, then:

B = (Q_out*S) * R_in = Q_out * (S * R_in) = Q_out * R_out.

So if one wants to interpret Q_out as the result of the QR factorization of B, then the corresponding R_out should be equal to R_out = S * R_in, i.e. some rows of R_in should be multiplied by -1. periodt. ✨

For the details of the algorithm, see [1].

[1] ,
    G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
    E. Solomonik, J. Parallel Distrib. Comput.,
    vol. 85, pp. 3-31, 2015.

👑 Author

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

🫶 Contributors

• November 2019, Igor Kozachenko, Computer Science Division, University of California, Berkeley

Definition at line 258 of file sorhr_col.f.

🌸 subroutine zunhr_col

subroutine zunhr_col (integer m, integer n, integer nb, complex*16,
    dimension( lda, * ) a, integer lda, complex*16, dimension( ldt, * ) t,
    integer ldt, complex*16, dimension( * ) d, integer info)

ZUNHR_COL

✨ Purpose

ZUNHR_COL takes an M-by-N complex matrix Q_in with orthonormal columns as input, stored in A, and performs Householder Reconstruction (HR), i.e. reconstructs Householder vectors V(i) implicitly representing another M-by-N matrix Q_out, with the property that Q_in = Q_out*S, where S is an N-by-N diagonal matrix with diagonal entries equal to +1 or -1. The Householder vectors (columns V(i) of V) are stored in A on output, and the diagonal entries of S are stored in D. Block reflectors are also returned in T (same output format as ZGEQRT). 🎀

🎀 Parameters

• M — M is INTEGER. The number of rows of the matrix A. M >= 0.
• N — N is INTEGER. The number of columns of the matrix A. M >= N >= 0.
• NB — NB is INTEGER. The column block size to be used in the reconstruction of Householder column vector blocks in the array A and corresponding block reflectors in the array T. NB >= 1. (Note that if NB > N, then N is used instead of NB as the column block size.)
• A — A is COMPLEX*16 array, dimension (LDA,N).

On entry:

The array A contains an M-by-N orthonormal matrix Q_in, i.e. the columns of A are orthogonal unit vectors.

On exit:

The elements below the diagonal of A represent the unit lower-trapezoidal matrix V of Householder column vectors V(i). The unit diagonal entries of V are not stored (same format as the output below the diagonal in A from ZGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out.

The elements above the diagonal contain the factor U of the LU-decomposition:

   Q_in - ( S ) = V * U
          ( 0 )

where 0 is a (M-N)-by-(M-N) zero matrix.

• LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
• T — T is COMPLEX*16 array, dimension (LDT, N).

Let NOCB = Number_of_output_col_blocks = CEIL(N/NB).

On exit, T(1:NB, 1:N) contains NOCB upper-triangular block reflectors used to define Q_out, stored in compact form as a sequence of upper-triangular NB-by-NB column blocks (same format as the output T in ZGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out. NOTE: the lower triangles below the upper-triangular blocks will be filled with zeros. See Further Details.

• LDT — LDT is INTEGER. The leading dimension of the array T. LDT >= max(1,min(NB,N)).
• D — D is COMPLEX*16 array, dimension min(M,N). The elements can be only plus or minus one.

D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where 1 <= i <= min(M,N), and Q_in_i is Q_in after performing i-1 steps of “modified” Gaussian elimination. See Further Details.

• INFO — INFO is INTEGER.

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

🧚 Further Details

okay so back to the glow-up, bestie 💅 — the computed M-by-M unitary factor Q_out is defined implicitly as a product of unitary matrices Q_out(i). Each Q_out(i) is stored in the compact WY-representation format in the corresponding blocks of matrices V (stored in A) and T.

The M-by-N unit lower-trapezoidal matrix V, stored in the M-by-N matrix A, contains the column vectors V(i) in NB-size column blocks VB(j). Example: VB(1) contains the columns V(1), V(2), ... V(NB). NOTE: the unit entries on the diagonal of Y are not stored in A — and yes, the source says Y there, not V, so I’m not gaslighting you 🙈.

The number of column blocks is

NOCB = Number_of_output_col_blocks = CEIL(N/NB)

where each block is of order NB except for the last block, which is of order LAST_NB = N - (NOCB-1)*NB.

For example, if M=6, N=5 and NB=2, the matrix V is

     V = (    VB(1),   VB(2), VB(3) ) =

       = (   1                      )
         ( v21    1                 )
         ( v31  v32    1            )
         ( v41  v42  v43   1        )
         ( v51  v52  v53  v54    1  )
         ( v61  v62  v63  v54   v65 )

For each of the column blocks VB(i), an upper-triangular block reflector TB(i) is computed. These blocks are stored as a sequence of upper-triangular column blocks in the NB-by-N matrix T. The size of each TB(i) block is NB-by-NB, except for the last block, whose size is LAST_NB-by-LAST_NB.

For example, if M=6, N=5 and NB=2, the matrix T is

     T  = (    TB(1),    TB(2), TB(3) ) =

        = ( t11  t12  t13  t14   t15  )
          (      t22       t24        )

The M-by-M factor Q_out is given as a product of NOCB unitary M-by-M matrices Q_out(i).

     Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),

where each matrix Q_out(i) is given by the WY-representation using corresponding blocks from the matrices V and T:

     Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,

where I is the identity matrix. Here is the formula with matrix dimensions:

  Q(i){M-by-M} = I{M-by-M} -
    VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},

where INB = NB, except for the last block NOCB for which INB=LAST_NB.

⚠️ NOTE — the source literally does the ===== / NOTE: / ===== thing, so here’s that energy, bestie 😭

If Q_in is the result of doing a QR factorization B = Q_in * R_in, then:

 B = (Q_out*S) * R_in = Q_out * (S * R_in) = Q_out * R_out.

So if one wants to interpret Q_out as the result of the QR factorization of B, then the corresponding R_out should be equal to R_out = S * R_in, i.e. some rows of R_in should be multiplied by -1.

For the details of the algorithm, see [1].

[1] ,
    G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
    E. Solomonik, J. Parallel Distrib. Comput.,
    vol. 85, pp. 3-31, 2015.

---

👑 Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

---

🌸 Contributors:

November 2019, Igor Kozachenko, Computer Science Division, University of California, Berkeley

---

Definition at line 258 of file zunhr_col.f.

🪄 Author

Generated automatically by Doxygen for LAPACK from the source code.

---

That’s the whole page, bestie 🫶 — the footer wraps it with LAPACK Version 3.12.0 for unhr_col(3), plus the date line, and it was generated automatically by Doxygen for LAPACK from the source code. Mwah 💋✨