/home/abui...rfb_gett.f(3)

πŸ“¦ version 3.12.0, in opensuse-leap 16.0

okay bestie, welcome to /home/abui...rfb_gett.f(3) β€” Library Functions Manual energy β€” and yes, we're keeping every detail intact. This is the DLARFB_GETT page, and it's giving very specific LAPACK realness πŸ’…βœ¨

πŸ’– NAME πŸ’–

/home/abuild/rpmbuild/BUILD/lapack-static-3.12.0-build/lapack-3.12.0/SRC/dlarfb_gett.f


✨ SYNOPSIS ✨

🧩 Functions/Subroutines 🧩

subroutine DLARFB_GETT (ident, m, n, k, t, ldt, a, lda, b, ldb, work,
    ldwork)
    DLARFB_GETT

πŸŽ€ Function/Subroutine Documentation πŸŽ€

subroutine DLARFB_GETT (character ident, integer m, integer n, integer k,
    double precision, dimension( ldt, * ) t, integer ldt, double precision,
    dimension( lda, * ) a, integer lda, double precision, dimension( ldb, *
    ) b, integer ldb, double precision, dimension( ldwork, * ) work, integer
    ldwork)

DLARFB_GETT

🌸 Purpose 🌸

okay so DLARFB_GETT applies a real Householder block reflector H from the left to a real (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 🌈

  • 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.
  • M β€” M is INTEGER.
    • The number of rows of the matrix B.
    • M >= 0.
  • N β€” N is INTEGER.
    • The number of columns of the matrices A and B.
    • N >= 0.
  • 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.
  • T β€” T is DOUBLE PRECISION array, dimension (LDT,K).
    • The upper-triangular K-by-K matrix T in the representation of the block reflector.
  • LDT β€” LDT is INTEGER.
    • The leading dimension of the array T.
    • LDT >= K.
  • A β€” A is DOUBLE PRECISION 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.
  • LDA β€” LDB is INTEGER. Yes bestie, this page really says LDB here β€” I'm keeping it exactly as written. The leading dimension of the array A. LDA >= max(1,K).
  • B β€” B is DOUBLE PRECISION 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.
  • LDB β€” LDB is INTEGER.
    • The leading dimension of the array B.
    • LDB >= max(1,M).
  • WORK β€” WORK is DOUBLE PRECISION array, dimension (LDWORK,max(K,N-K)).
  • 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 🧚

Okay bestie, now let's get into the juicy details of the subroutine documentation. It's giving very detailed, very iconic. πŸ’…

πŸ’– (1) Description of the Algebraic Operation πŸ’–

So, 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). So like, A = ( A1, A2 ). And 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). So B = ( B1, B2 ). Slay.

We perform the operation:

   ( A_out ) := H * ( A_in ) = ( I - V * T * V**T ) * ( A_in ) =
   ( B_out )        ( B_in )                          ( B_in )
              = ( I - ( V1 ) * T * ( V1**T, V2**T ) ) * ( 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 ) over ( V2 )

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 matrix stored in 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**T ) * ( A1_in )
   ( B1_out )        (     0 )                          (     0 )

   ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**T ) * ( A2_in )
   ( B2_out )        ( B2_in )                          ( B2_in )

If IDENT != 'I':

The computation for column block 1:

A1_out: = A1_in - V1*T*(V1**T)*A1_in

B1_out: = - V2*T*(V1**T)*A1_in

The computation for column block 2, which exists if N > K:

A2_out: = A2_in - V1*T*( (V1**T)*A2_in + (V2**T)*B2_in )

B2_out: = B2_in - V2*T*( (V1**T)*A2_in + (V2**T)*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**T)*B2_in )

B2_out: = B2_in - V2*T*( A2_in + (V2**T)*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 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).

Here's the exact algorithm for Case (A):

if( IDENT != 'I' ) then
 if ( N > K ) then
   (First Step - column block 2)
   col2_(1) W2: = A2
   col2_(2) W2: = (V1**T) * W2 = (unit_lower_tr_of_(A1)**T) * W2
   col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W1 = (unit_lower_tr_of_(A1)**T) * 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)

Here's the exact algorithm for Case (B):

if( IDENT == 'I' ) then
 if ( N > K ) then
   (First Step - column block 2)
   col2_(1) W2: = A2
   col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W2
                  = (unit_lower_tr_of_(A1)**T) * W2
  end if
  col2_(3)  W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W1
                = (unit_lower_tr_of_(A1)**T) * 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

And that's the whole algorithm, bestie! ✨ The definition is at line 390 of file dlarfb_gett.f.

Author: Generated automatically by Doxygen for LAPACK from the source code.


This iconic page was brought to you by LAPACK Version 3.12.0, from the file /home/abui...rfb_gett.f(3). Stay slay, bestie! πŸ’…πŸ’–

πŸ“¦ where this page comes from