 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dlarfb_gett()

 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:
``` DLARFB_GETT applies a real Householder block reflector H from the
left to a real (K+M)-by-N  "triangular-pentagonal" 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 DOUBLE PRECISION 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 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.``` [in] LDA ``` LDB is INTEGER The leading dimension of the array A. LDA >= max(1,K).``` [in,out] 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.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,M).``` [out] WORK ``` WORK is DOUBLE PRECISION array, dimension (LDWORK,max(K,N-K))``` [in] LDWORK ``` LDWORK is INTEGER The leading dimension of the array WORK. LDWORK>=max(1,K).```
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**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 )
( 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**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 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**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)

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```

Definition at line 390 of file dlarfb_gett.f.

392  IMPLICIT NONE
393 *
394 * -- LAPACK auxiliary routine --
395 * -- LAPACK is a software package provided by Univ. of Tennessee, --
396 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
397 *
398 * .. Scalar Arguments ..
399  CHARACTER IDENT
400  INTEGER K, LDA, LDB, LDT, LDWORK, M, N
401 * ..
402 * .. Array Arguments ..
403  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
404  \$ WORK( LDWORK, * )
405 * ..
406 *
407 * =====================================================================
408 *
409 * .. Parameters ..
410  DOUBLE PRECISION ONE, ZERO
411  parameter( one = 1.0d+0, zero = 0.0d+0 )
412 * ..
413 * .. Local Scalars ..
414  LOGICAL LNOTIDENT
415  INTEGER I, J
416 * ..
417 * .. EXTERNAL FUNCTIONS ..
418  LOGICAL LSAME
419  EXTERNAL lsame
420 * ..
421 * .. External Subroutines ..
422  EXTERNAL dcopy, dgemm, dtrmm
423 * ..
424 * .. Executable Statements ..
425 *
426 * Quick return if possible
427 *
428  IF( m.LT.0 .OR. n.LE.0 .OR. k.EQ.0 .OR. k.GT.n )
429  \$ RETURN
430 *
431  lnotident = .NOT.lsame( ident, 'I' )
432 *
433 * ------------------------------------------------------------------
434 *
435 * First Step. Computation of the Column Block 2:
436 *
437 * ( A2 ) := H * ( A2 )
438 * ( B2 ) ( B2 )
439 *
440 * ------------------------------------------------------------------
441 *
442  IF( n.GT.k ) THEN
443 *
444 * col2_(1) Compute W2: = A2. Therefore, copy A2 = A(1:K, K+1:N)
445 * into W2=WORK(1:K, 1:N-K) column-by-column.
446 *
447  DO j = 1, n-k
448  CALL dcopy( k, a( 1, k+j ), 1, work( 1, j ), 1 )
449  END DO
450
451  IF( lnotident ) THEN
452 *
453 * col2_(2) Compute W2: = (V1**T) * W2 = (A1**T) * W2,
454 * V1 is not an identy matrix, but unit lower-triangular
455 * V1 stored in A1 (diagonal ones are not stored).
456 *
457 *
458  CALL dtrmm( 'L', 'L', 'T', 'U', k, n-k, one, a, lda,
459  \$ work, ldwork )
460  END IF
461 *
462 * col2_(3) Compute W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
463 * V2 stored in B1.
464 *
465  IF( m.GT.0 ) THEN
466  CALL dgemm( 'T', 'N', k, n-k, m, one, b, ldb,
467  \$ b( 1, k+1 ), ldb, one, work, ldwork )
468  END IF
469 *
470 * col2_(4) Compute W2: = T * W2,
471 * T is upper-triangular.
472 *
473  CALL dtrmm( 'L', 'U', 'N', 'N', k, n-k, one, t, ldt,
474  \$ work, ldwork )
475 *
476 * col2_(5) Compute B2: = B2 - V2 * W2 = B2 - B1 * W2,
477 * V2 stored in B1.
478 *
479  IF( m.GT.0 ) THEN
480  CALL dgemm( 'N', 'N', m, n-k, k, -one, b, ldb,
481  \$ work, ldwork, one, b( 1, k+1 ), ldb )
482  END IF
483 *
484  IF( lnotident ) THEN
485 *
486 * col2_(6) Compute W2: = V1 * W2 = A1 * W2,
487 * V1 is not an identity matrix, but unit lower-triangular,
488 * V1 stored in A1 (diagonal ones are not stored).
489 *
490  CALL dtrmm( 'L', 'L', 'N', 'U', k, n-k, one, a, lda,
491  \$ work, ldwork )
492  END IF
493 *
494 * col2_(7) Compute A2: = A2 - W2 =
495 * = A(1:K, K+1:N-K) - WORK(1:K, 1:N-K),
496 * column-by-column.
497 *
498  DO j = 1, n-k
499  DO i = 1, k
500  a( i, k+j ) = a( i, k+j ) - work( i, j )
501  END DO
502  END DO
503 *
504  END IF
505 *
506 * ------------------------------------------------------------------
507 *
508 * Second Step. Computation of the Column Block 1:
509 *
510 * ( A1 ) := H * ( A1 )
511 * ( B1 ) ( 0 )
512 *
513 * ------------------------------------------------------------------
514 *
515 * col1_(1) Compute W1: = A1. Copy the upper-triangular
516 * A1 = A(1:K, 1:K) into the upper-triangular
517 * W1 = WORK(1:K, 1:K) column-by-column.
518 *
519  DO j = 1, k
520  CALL dcopy( j, a( 1, j ), 1, work( 1, j ), 1 )
521  END DO
522 *
523 * Set the subdiagonal elements of W1 to zero column-by-column.
524 *
525  DO j = 1, k - 1
526  DO i = j + 1, k
527  work( i, j ) = zero
528  END DO
529  END DO
530 *
531  IF( lnotident ) THEN
532 *
533 * col1_(2) Compute W1: = (V1**T) * W1 = (A1**T) * W1,
534 * V1 is not an identity matrix, but unit lower-triangular
535 * V1 stored in A1 (diagonal ones are not stored),
536 * W1 is upper-triangular with zeroes below the diagonal.
537 *
538  CALL dtrmm( 'L', 'L', 'T', 'U', k, k, one, a, lda,
539  \$ work, ldwork )
540  END IF
541 *
542 * col1_(3) Compute W1: = T * W1,
543 * T is upper-triangular,
544 * W1 is upper-triangular with zeroes below the diagonal.
545 *
546  CALL dtrmm( 'L', 'U', 'N', 'N', k, k, one, t, ldt,
547  \$ work, ldwork )
548 *
549 * col1_(4) Compute B1: = - V2 * W1 = - B1 * W1,
550 * V2 = B1, W1 is upper-triangular with zeroes below the diagonal.
551 *
552  IF( m.GT.0 ) THEN
553  CALL dtrmm( 'R', 'U', 'N', 'N', m, k, -one, work, ldwork,
554  \$ b, ldb )
555  END IF
556 *
557  IF( lnotident ) THEN
558 *
559 * col1_(5) Compute W1: = V1 * W1 = A1 * W1,
560 * V1 is not an identity matrix, but unit lower-triangular
561 * V1 stored in A1 (diagonal ones are not stored),
562 * W1 is upper-triangular on input with zeroes below the diagonal,
563 * and square on output.
564 *
565  CALL dtrmm( 'L', 'L', 'N', 'U', k, k, one, a, lda,
566  \$ work, ldwork )
567 *
568 * col1_(6) Compute A1: = A1 - W1 = A(1:K, 1:K) - WORK(1:K, 1:K)
569 * column-by-column. A1 is upper-triangular on input.
570 * If IDENT, A1 is square on output, and W1 is square,
571 * if NOT IDENT, A1 is upper-triangular on output,
572 * W1 is upper-triangular.
573 *
574 * col1_(6)_a Compute elements of A1 below the diagonal.
575 *
576  DO j = 1, k - 1
577  DO i = j + 1, k
578  a( i, j ) = - work( i, j )
579  END DO
580  END DO
581 *
582  END IF
583 *
584 * col1_(6)_b Compute elements of A1 on and above the diagonal.
585 *
586  DO j = 1, k
587  DO i = 1, j
588  a( i, j ) = a( i, j ) - work( i, j )
589  END DO
590  END DO
591 *
592  RETURN
593 *
594 * End of DLARFB_GETT
595 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: