LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ 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 !> !>

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