LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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

Download DLARFB_GETT + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> 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
[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).
!>
!>
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**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 388 of file dlarfb_gett.f.

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