LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ zlarfb_gett()

subroutine zlarfb_gett ( character  IDENT,
integer  M,
integer  N,
integer  K,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldwork, * )  WORK,
integer  LDWORK 
)

ZLARFB_GETT

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

Purpose:
 ZLARFB_GETT applies a complex Householder block reflector H from the
 left to a complex (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 COMPLEX*16 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*16 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*16 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*16 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

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