Loading...
Searching...
No Matches

## clarfb_gett()

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]

Purpose:

      !>
      !> CLARFB_GETT applies a complex Householder block reflector H from the
      !> left to a complex (K+M)-by-N   matrix
      !> composed of two block matrices: an upper trapezoidal K-by-N matrix A
      !> stored in the array A, and a rectangular M-by-(N-K) matrix B, stored
      !> in the array B. The block reflector H is stored in a compact
      !> WY-representation, where the elementary reflectors are in the
      !> arrays A, B and T. See Further Details section.
      !>

Parameters

  [in]

  IDENT

      !>          IDENT is CHARACTER*1
      !>          If IDENT = not 'I', or not 'i', then V1 is unit
      !>             lower-triangular and stored in the left K-by-K 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 upper-triangular K-by-K 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 K-by-N upper-trapezoidal 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 K-by-N 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 M-by-(N-K) right block: input matrix B.
      !>            b) In the M-by-N left block: columns of V2.
      !>
      !>          On exit:
      !>            B is overwritten by rectangular M-by-N 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,N-K))
      !>

  [in]

  LDWORK

      !>          LDWORK is INTEGER
      !>          The leading dimension of the array WORK. LDWORK>=max(1,K).
      !>
      !>

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

Contributors:

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

Further Details:

      !>
      !>    (1) Description of the Algebraic Operation.
      !>
      !>    The matrix A is a K-by-N matrix composed of two column block
      !>    matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):
      !>    A = ( A1, A2 ).
      !>    The matrix B is an M-by-N matrix composed of two column block
      !>    matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):
      !>    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 K-by-K upper-triangular matrix stored in the
      !>                   upper triangular part of the array A(1:K,1:K).
      !>       (  B1_in )  is an M-by-K rectangular ZERO matrix and not stored.
      !>
      !>       ( A2_in )  is a K-by-(N-K) rectangular matrix stored
      !>                  in the array A(1:K,K+1:N).
      !>       ( B2_in )  is an M-by-(N-K) rectangular matrix stored
      !>                  in the array B(1:M,K+1:N).
      !>
      !>    b) V = ( V1 )
      !>           ( V2 )
      !>
      !>       where:
      !>       1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;
      !>       2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,
      !>          stored in the lower-triangular part of the array
      !>          A(1:K,1:K) (ones are not stored),
      !>       and V2 is an M-by-K 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 K-by-K upper-triangular 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 K-by-K square matrix, or a K-by-K
      !>                   upper-triangular matrix, if V1 is an
      !>                   identity matrix. AiOut is stored in
      !>                   the array A(1:K,1:K).
      !>       ( B1_out )  is an M-by-K rectangular matrix stored
      !>                   in the array B(1:M,K:N).
      !>
      !>       ( A2_out )  is a K-by-(N-K) rectangular matrix stored
      !>                   in the array A(1:K,K+1:N).
      !>       ( B2_out )  is an M-by-(N-K) 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 K-by-(N-K) rectangular workspace
      !>    matrix W2 that is of the same size as the matrix A2.
      !>    W2 is stored in the array WORK(1:K,1:(N-K)).
      !>
      !>    In the second step, we compute column block 1, i.e. A1 and B1.
      !>    Here, we need to use the K-by-K 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,N-K)) so it can hold both W2 from
      !>    the first step and W1 from the second step.
      !>
      !>    Case (A), when V1 is unit lower-triangular, 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) upper-triangular_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
      !>
      !>