LAPACK 3.12.0
LAPACK: Linear Algebra PACKage

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
Download DORHR_COL + dependencies [TGZ] [ZIP] [TXT]
DORHR_COL takes an MbyN 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 MbyN matrix Q_out, with the property that Q_in = Q_out*S, where S is an NbyN 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).
[in]  M  M is INTEGER The number of rows of the matrix A. M >= 0. 
[in]  N  N is INTEGER The number of columns of the matrix A. M >= N >= 0. 
[in]  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.) 
[in,out]  A  A is DOUBLE PRECISION array, dimension (LDA,N) On entry: The array A contains an MbyN 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 lowertrapezoidal 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 "modified" LUdecomposition: Q_in  ( S ) = V * U ( 0 ) where 0 is a (MN)by(MN) zero matrix. 
[in]  LDA  LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M). 
[out]  T  T is 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 uppertriangular block reflectors used to define Q_out stored in compact form as a sequence of uppertriangular NBbyNB 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 uppertriangular blocks will be filled with zeros. See Further Details. 
[in]  LDT  LDT is INTEGER The leading dimension of the array T. LDT >= max(1,min(NB,N)). 
[out]  D  D is 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 i1 steps of “modified” Gaussian elimination. See Further Details. 
[out]  INFO  INFO is INTEGER = 0: successful exit < 0: if INFO = i, the ith argument had an illegal value 
The computed MbyM 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 WYrepresentation format in the corresponding blocks of matrices V (stored in A) and T. The MbyN unit lowertrapezoidal matrix V stored in the MbyN matrix A contains the column vectors V(i) in NBsize 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  (NOCB1)*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 uppertriangular block reflector TB(i) is computed. These blocks are stored as a sequence of uppertriangular column blocks in the NBbyN matrix T. The size of each TB(i) block is NBbyNB, except for the last block, whose size is LAST_NBbyLAST_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 MbyM factor Q_out is given as a product of NOCB orthogonal MbyM 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 WYrepresentation 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){MbyM} = I{MbyM}  VB(i){MbyINB} * TB(i){INBbyINB} * (VB(i))**T {INBbyM}, 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] "Reconstructing Householder vectors from tallskinny QR", G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen, E. Solomonik, J. Parallel Distrib. Comput., vol. 85, pp. 331, 2015.
November 2019, Igor Kozachenko, Computer Science Division, University of California, Berkeley
Definition at line 258 of file dorhr_col.f.