LAPACK 3.12.0
LAPACK: Linear Algebra PACKage

subroutine clarfb_gett  (  character  ident, 
integer  m,  
integer  n,  
integer  k,  
complex, dimension( ldt, * )  t,  
integer  ldt,  
complex, dimension( lda, * )  a,  
integer  lda,  
complex, dimension( ldb, * )  b,  
integer  ldb,  
complex, dimension( ldwork, * )  work,  
integer  ldwork  
) 
CLARFB_GETT
Download CLARFB_GETT + dependencies [TGZ] [ZIP] [TXT]
CLARFB_GETT applies a complex Householder block reflector H from the left to a complex (K+M)byN "triangularpentagonal" matrix composed of two block matrices: an upper trapezoidal KbyN matrix A stored in the array A, and a rectangular Mby(NK) matrix B, stored in the array B. The block reflector H is stored in a compact WYrepresentation, where the elementary reflectors are in the arrays A, B and T. See Further Details section.
[in]  IDENT  IDENT is CHARACTER*1 If IDENT = not 'I', or not 'i', then V1 is unit lowertriangular and stored in the left KbyK block of the input matrix A, If IDENT = 'I' or 'i', then V1 is an identity matrix and not stored. See Further Details section. 
[in]  M  M is INTEGER The number of rows of the matrix B. M >= 0. 
[in]  N  N is INTEGER The number of columns of the matrices A and B. N >= 0. 
[in]  K  K is INTEGER The number or rows of the matrix A. K is also order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. 0 <= K <= N. 
[in]  T  T is COMPLEX array, dimension (LDT,K) The uppertriangular KbyK matrix T in the representation of the block reflector. 
[in]  LDT  LDT is INTEGER The leading dimension of the array T. LDT >= K. 
[in,out]  A  A is COMPLEX array, dimension (LDA,N) On entry: a) In the KbyN uppertrapezoidal part A: input matrix A. b) In the columns below the diagonal: columns of V1 (ones are not stored on the diagonal). On exit: A is overwritten by rectangular KbyN product H*A. See Further Details section. 
[in]  LDA  LDB is INTEGER The leading dimension of the array A. LDA >= max(1,K). 
[in,out]  B  B is COMPLEX array, dimension (LDB,N) On entry: a) In the Mby(NK) right block: input matrix B. b) In the MbyN left block: columns of V2. On exit: B is overwritten by rectangular MbyN product H*B. See Further Details section. 
[in]  LDB  LDB is INTEGER The leading dimension of the array B. LDB >= max(1,M). 
[out]  WORK  WORK is COMPLEX array, dimension (LDWORK,max(K,NK)) 
[in]  LDWORK  LDWORK is INTEGER The leading dimension of the array WORK. LDWORK>=max(1,K). 
November 2020, Igor Kozachenko, Computer Science Division, University of California, Berkeley
(1) Description of the Algebraic Operation. The matrix A is a KbyN matrix composed of two column block matrices, A1, which is KbyK, and A2, which is Kby(NK): A = ( A1, A2 ). The matrix B is an MbyN matrix composed of two column block matrices, B1, which is MbyK, and B2, which is Mby(NK): B = ( B1, B2 ). Perform the operation: ( A_out ) := H * ( A_in ) = ( I  V * T * V**H ) * ( A_in ) = ( B_out ) ( B_in ) ( B_in ) = ( I  ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in ) ( V2 ) ( B_in ) On input: a) ( A_in ) consists of two block columns: ( B_in ) ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in )) ( B_in ) (( B1_in ) ( B2_in )) (( 0 ) ( B2_in )), where the column blocks are: ( A1_in ) is a KbyK uppertriangular matrix stored in the upper triangular part of the array A(1:K,1:K). ( B1_in ) is an MbyK rectangular ZERO matrix and not stored. ( A2_in ) is a Kby(NK) rectangular matrix stored in the array A(1:K,K+1:N). ( B2_in ) is an Mby(NK) rectangular matrix stored in the array B(1:M,K+1:N). b) V = ( V1 ) ( V2 ) where: 1) if IDENT == 'I',V1 is a KbyK identity matrix, not stored; 2) if IDENT != 'I',V1 is a KbyK unit lowertriangular matrix, stored in the lowertriangular part of the array A(1:K,1:K) (ones are not stored), and V2 is an MbyK rectangular stored the array B(1:M,1:K), (because on input B1_in is a rectangular zero matrix that is not stored and the space is used to store V2). c) T is a KbyK uppertriangular matrix stored in the array T(1:K,1:K). On output: a) ( A_out ) consists of two block columns: ( B_out ) ( A_out ) = (( A1_out ) ( A2_out )) ( B_out ) (( B1_out ) ( B2_out )), where the column blocks are: ( A1_out ) is a KbyK square matrix, or a KbyK uppertriangular matrix, if V1 is an identity matrix. AiOut is stored in the array A(1:K,1:K). ( B1_out ) is an MbyK rectangular matrix stored in the array B(1:M,K:N). ( A2_out ) is a Kby(NK) rectangular matrix stored in the array A(1:K,K+1:N). ( B2_out ) is an Mby(NK) rectangular matrix stored in the array B(1:M,K+1:N). The operation above can be represented as the same operation on each block column: ( A1_out ) := H * ( A1_in ) = ( I  V * T * V**H ) * ( A1_in ) ( B1_out ) ( 0 ) ( 0 ) ( A2_out ) := H * ( A2_in ) = ( I  V * T * V**H ) * ( A2_in ) ( B2_out ) ( B2_in ) ( B2_in ) If IDENT != 'I': The computation for column block 1: A1_out: = A1_in  V1*T*(V1**H)*A1_in B1_out: =  V2*T*(V1**H)*A1_in The computation for column block 2, which exists if N > K: A2_out: = A2_in  V1*T*( (V1**H)*A2_in + (V2**H)*B2_in ) B2_out: = B2_in  V2*T*( (V1**H)*A2_in + (V2**H)*B2_in ) If IDENT == 'I': The operation for column block 1: A1_out: = A1_in  V1*T*A1_in B1_out: =  V2*T*A1_in The computation for column block 2, which exists if N > K: A2_out: = A2_in  T*( A2_in + (V2**H)*B2_in ) B2_out: = B2_in  V2*T*( A2_in + (V2**H)*B2_in ) (2) Description of the Algorithmic Computation. In the first step, we compute column block 2, i.e. A2 and B2. Here, we need to use the Kby(NK) rectangular workspace matrix W2 that is of the same size as the matrix A2. W2 is stored in the array WORK(1:K,1:(NK)). In the second step, we compute column block 1, i.e. A1 and B1. Here, we need to use the KbyK square workspace matrix W1 that is of the same size as the as the matrix A1. W1 is stored in the array WORK(1:K,1:K). NOTE: Hence, in this routine, we need the workspace array WORK only of size WORK(1:K,1:max(K,NK)) so it can hold both W2 from the first step and W1 from the second step. Case (A), when V1 is unit lowertriangular, i.e. IDENT != 'I', more computations than in the Case (B). if( IDENT != 'I' ) then if ( N > K ) then (First Step  column block 2) col2_(1) W2: = A2 col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2 col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2 col2_(4) W2: = T * W2 col2_(5) B2: = B2  V2 * W2 = B2  B1 * W2 col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2 col2_(7) A2: = A2  W2 else (Second Step  column block 1) col1_(1) W1: = A1 col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1 col1_(3) W1: = T * W1 col1_(4) B1: =  V2 * W1 =  B1 * W1 col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1 col1_(6) square A1: = A1  W1 end if end if Case (B), when V1 is an identity matrix, i.e. IDENT == 'I', less computations than in the Case (A) if( IDENT == 'I' ) then if ( N > K ) then (First Step  column block 2) col2_(1) W2: = A2 col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2 col2_(4) W2: = T * W2 col2_(5) B2: = B2  V2 * W2 = B2  B1 * W2 col2_(7) A2: = A2  W2 else (Second Step  column block 1) col1_(1) W1: = A1 col1_(3) W1: = T * W1 col1_(4) B1: =  V2 * W1 =  B1 * W1 col1_(6) uppertriangular_of_(A1): = A1  W1 end if end if Combine these cases (A) and (B) together, this is the resulting algorithm: if ( N > K ) then (First Step  column block 2) col2_(1) W2: = A2 if( IDENT != 'I' ) then col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2 end if col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2] col2_(4) W2: = T * W2 col2_(5) B2: = B2  V2 * W2 = B2  B1 * W2 if( IDENT != 'I' ) then col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2 end if col2_(7) A2: = A2  W2 else (Second Step  column block 1) col1_(1) W1: = A1 if( IDENT != 'I' ) then col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1 end if col1_(3) W1: = T * W1 col1_(4) B1: =  V2 * W1 =  B1 * W1 if( IDENT != 'I' ) then col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1 col1_(6_a) below_diag_of_(A1): =  below_diag_of_(W1) end if col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1)  up_tr_of_(W1) end if
Definition at line 390 of file clarfb_gett.f.