1*> \brief \b ZLARFB_GETT 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLARFB_GETT + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfb_gett.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfb_gett.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfb_gett.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLARFB_GETT( IDENT, M, N, K, T, LDT, A, LDA, B, LDB, 22* $ WORK, LDWORK ) 23* IMPLICIT NONE 24* 25* .. Scalar Arguments .. 26* CHARACTER IDENT 27* INTEGER K, LDA, LDB, LDT, LDWORK, M, N 28* .. 29* .. Array Arguments .. 30* COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ), 31* $ WORK( LDWORK, * ) 32* .. 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> ZLARFB_GETT applies a complex Householder block reflector H from the 40*> left to a complex (K+M)-by-N "triangular-pentagonal" matrix 41*> composed of two block matrices: an upper trapezoidal K-by-N matrix A 42*> stored in the array A, and a rectangular M-by-(N-K) matrix B, stored 43*> in the array B. The block reflector H is stored in a compact 44*> WY-representation, where the elementary reflectors are in the 45*> arrays A, B and T. See Further Details section. 46*> \endverbatim 47* 48* Arguments: 49* ========== 50* 51*> \param[in] IDENT 52*> \verbatim 53*> IDENT is CHARACTER*1 54*> If IDENT = not 'I', or not 'i', then V1 is unit 55*> lower-triangular and stored in the left K-by-K block of 56*> the input matrix A, 57*> If IDENT = 'I' or 'i', then V1 is an identity matrix and 58*> not stored. 59*> See Further Details section. 60*> \endverbatim 61*> 62*> \param[in] M 63*> \verbatim 64*> M is INTEGER 65*> The number of rows of the matrix B. 66*> M >= 0. 67*> \endverbatim 68*> 69*> \param[in] N 70*> \verbatim 71*> N is INTEGER 72*> The number of columns of the matrices A and B. 73*> N >= 0. 74*> \endverbatim 75*> 76*> \param[in] K 77*> \verbatim 78*> K is INTEGER 79*> The number or rows of the matrix A. 80*> K is also order of the matrix T, i.e. the number of 81*> elementary reflectors whose product defines the block 82*> reflector. 0 <= K <= N. 83*> \endverbatim 84*> 85*> \param[in] T 86*> \verbatim 87*> T is COMPLEX*16 array, dimension (LDT,K) 88*> The upper-triangular K-by-K matrix T in the representation 89*> of the block reflector. 90*> \endverbatim 91*> 92*> \param[in] LDT 93*> \verbatim 94*> LDT is INTEGER 95*> The leading dimension of the array T. LDT >= K. 96*> \endverbatim 97*> 98*> \param[in,out] A 99*> \verbatim 100*> A is COMPLEX*16 array, dimension (LDA,N) 101*> 102*> On entry: 103*> a) In the K-by-N upper-trapezoidal part A: input matrix A. 104*> b) In the columns below the diagonal: columns of V1 105*> (ones are not stored on the diagonal). 106*> 107*> On exit: 108*> A is overwritten by rectangular K-by-N product H*A. 109*> 110*> See Further Details section. 111*> \endverbatim 112*> 113*> \param[in] LDA 114*> \verbatim 115*> LDB is INTEGER 116*> The leading dimension of the array A. LDA >= max(1,K). 117*> \endverbatim 118*> 119*> \param[in,out] B 120*> \verbatim 121*> B is COMPLEX*16 array, dimension (LDB,N) 122*> 123*> On entry: 124*> a) In the M-by-(N-K) right block: input matrix B. 125*> b) In the M-by-N left block: columns of V2. 126*> 127*> On exit: 128*> B is overwritten by rectangular M-by-N product H*B. 129*> 130*> See Further Details section. 131*> \endverbatim 132*> 133*> \param[in] LDB 134*> \verbatim 135*> LDB is INTEGER 136*> The leading dimension of the array B. LDB >= max(1,M). 137*> \endverbatim 138*> 139*> \param[out] WORK 140*> \verbatim 141*> WORK is COMPLEX*16 array, 142*> dimension (LDWORK,max(K,N-K)) 143*> \endverbatim 144*> 145*> \param[in] LDWORK 146*> \verbatim 147*> LDWORK is INTEGER 148*> The leading dimension of the array WORK. LDWORK>=max(1,K). 149*> 150*> \endverbatim 151* 152* Authors: 153* ======== 154* 155*> \author Univ. of Tennessee 156*> \author Univ. of California Berkeley 157*> \author Univ. of Colorado Denver 158*> \author NAG Ltd. 159* 160*> \ingroup complex16OTHERauxiliary 161* 162*> \par Contributors: 163* ================== 164*> 165*> \verbatim 166*> 167*> November 2020, Igor Kozachenko, 168*> Computer Science Division, 169*> University of California, Berkeley 170*> 171*> \endverbatim 172* 173*> \par Further Details: 174* ===================== 175*> 176*> \verbatim 177*> 178*> (1) Description of the Algebraic Operation. 179*> 180*> The matrix A is a K-by-N matrix composed of two column block 181*> matrices, A1, which is K-by-K, and A2, which is K-by-(N-K): 182*> A = ( A1, A2 ). 183*> The matrix B is an M-by-N matrix composed of two column block 184*> matrices, B1, which is M-by-K, and B2, which is M-by-(N-K): 185*> B = ( B1, B2 ). 186*> 187*> Perform the operation: 188*> 189*> ( A_out ) := H * ( A_in ) = ( I - V * T * V**H ) * ( A_in ) = 190*> ( B_out ) ( B_in ) ( B_in ) 191*> = ( I - ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in ) 192*> ( V2 ) ( B_in ) 193*> On input: 194*> 195*> a) ( A_in ) consists of two block columns: 196*> ( B_in ) 197*> 198*> ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in )) 199*> ( B_in ) (( B1_in ) ( B2_in )) (( 0 ) ( B2_in )), 200*> 201*> where the column blocks are: 202*> 203*> ( A1_in ) is a K-by-K upper-triangular matrix stored in the 204*> upper triangular part of the array A(1:K,1:K). 205*> ( B1_in ) is an M-by-K rectangular ZERO matrix and not stored. 206*> 207*> ( A2_in ) is a K-by-(N-K) rectangular matrix stored 208*> in the array A(1:K,K+1:N). 209*> ( B2_in ) is an M-by-(N-K) rectangular matrix stored 210*> in the array B(1:M,K+1:N). 211*> 212*> b) V = ( V1 ) 213*> ( V2 ) 214*> 215*> where: 216*> 1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored; 217*> 2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix, 218*> stored in the lower-triangular part of the array 219*> A(1:K,1:K) (ones are not stored), 220*> and V2 is an M-by-K rectangular stored the array B(1:M,1:K), 221*> (because on input B1_in is a rectangular zero 222*> matrix that is not stored and the space is 223*> used to store V2). 224*> 225*> c) T is a K-by-K upper-triangular matrix stored 226*> in the array T(1:K,1:K). 227*> 228*> On output: 229*> 230*> a) ( A_out ) consists of two block columns: 231*> ( B_out ) 232*> 233*> ( A_out ) = (( A1_out ) ( A2_out )) 234*> ( B_out ) (( B1_out ) ( B2_out )), 235*> 236*> where the column blocks are: 237*> 238*> ( A1_out ) is a K-by-K square matrix, or a K-by-K 239*> upper-triangular matrix, if V1 is an 240*> identity matrix. AiOut is stored in 241*> the array A(1:K,1:K). 242*> ( B1_out ) is an M-by-K rectangular matrix stored 243*> in the array B(1:M,K:N). 244*> 245*> ( A2_out ) is a K-by-(N-K) rectangular matrix stored 246*> in the array A(1:K,K+1:N). 247*> ( B2_out ) is an M-by-(N-K) rectangular matrix stored 248*> in the array B(1:M,K+1:N). 249*> 250*> 251*> The operation above can be represented as the same operation 252*> on each block column: 253*> 254*> ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**H ) * ( A1_in ) 255*> ( B1_out ) ( 0 ) ( 0 ) 256*> 257*> ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**H ) * ( A2_in ) 258*> ( B2_out ) ( B2_in ) ( B2_in ) 259*> 260*> If IDENT != 'I': 261*> 262*> The computation for column block 1: 263*> 264*> A1_out: = A1_in - V1*T*(V1**H)*A1_in 265*> 266*> B1_out: = - V2*T*(V1**H)*A1_in 267*> 268*> The computation for column block 2, which exists if N > K: 269*> 270*> A2_out: = A2_in - V1*T*( (V1**H)*A2_in + (V2**H)*B2_in ) 271*> 272*> B2_out: = B2_in - V2*T*( (V1**H)*A2_in + (V2**H)*B2_in ) 273*> 274*> If IDENT == 'I': 275*> 276*> The operation for column block 1: 277*> 278*> A1_out: = A1_in - V1*T*A1_in 279*> 280*> B1_out: = - V2*T*A1_in 281*> 282*> The computation for column block 2, which exists if N > K: 283*> 284*> A2_out: = A2_in - T*( A2_in + (V2**H)*B2_in ) 285*> 286*> B2_out: = B2_in - V2*T*( A2_in + (V2**H)*B2_in ) 287*> 288*> (2) Description of the Algorithmic Computation. 289*> 290*> In the first step, we compute column block 2, i.e. A2 and B2. 291*> Here, we need to use the K-by-(N-K) rectangular workspace 292*> matrix W2 that is of the same size as the matrix A2. 293*> W2 is stored in the array WORK(1:K,1:(N-K)). 294*> 295*> In the second step, we compute column block 1, i.e. A1 and B1. 296*> Here, we need to use the K-by-K square workspace matrix W1 297*> that is of the same size as the as the matrix A1. 298*> W1 is stored in the array WORK(1:K,1:K). 299*> 300*> NOTE: Hence, in this routine, we need the workspace array WORK 301*> only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from 302*> the first step and W1 from the second step. 303*> 304*> Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I', 305*> more computations than in the Case (B). 306*> 307*> if( IDENT != 'I' ) then 308*> if ( N > K ) then 309*> (First Step - column block 2) 310*> col2_(1) W2: = A2 311*> col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2 312*> col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2 313*> col2_(4) W2: = T * W2 314*> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2 315*> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2 316*> col2_(7) A2: = A2 - W2 317*> else 318*> (Second Step - column block 1) 319*> col1_(1) W1: = A1 320*> col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1 321*> col1_(3) W1: = T * W1 322*> col1_(4) B1: = - V2 * W1 = - B1 * W1 323*> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1 324*> col1_(6) square A1: = A1 - W1 325*> end if 326*> end if 327*> 328*> Case (B), when V1 is an identity matrix, i.e. IDENT == 'I', 329*> less computations than in the Case (A) 330*> 331*> if( IDENT == 'I' ) then 332*> if ( N > K ) then 333*> (First Step - column block 2) 334*> col2_(1) W2: = A2 335*> col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2 336*> col2_(4) W2: = T * W2 337*> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2 338*> col2_(7) A2: = A2 - W2 339*> else 340*> (Second Step - column block 1) 341*> col1_(1) W1: = A1 342*> col1_(3) W1: = T * W1 343*> col1_(4) B1: = - V2 * W1 = - B1 * W1 344*> col1_(6) upper-triangular_of_(A1): = A1 - W1 345*> end if 346*> end if 347*> 348*> Combine these cases (A) and (B) together, this is the resulting 349*> algorithm: 350*> 351*> if ( N > K ) then 352*> 353*> (First Step - column block 2) 354*> 355*> col2_(1) W2: = A2 356*> if( IDENT != 'I' ) then 357*> col2_(2) W2: = (V1**H) * W2 358*> = (unit_lower_tr_of_(A1)**H) * W2 359*> end if 360*> col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2] 361*> col2_(4) W2: = T * W2 362*> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2 363*> if( IDENT != 'I' ) then 364*> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2 365*> end if 366*> col2_(7) A2: = A2 - W2 367*> 368*> else 369*> 370*> (Second Step - column block 1) 371*> 372*> col1_(1) W1: = A1 373*> if( IDENT != 'I' ) then 374*> col1_(2) W1: = (V1**H) * W1 375*> = (unit_lower_tr_of_(A1)**H) * W1 376*> end if 377*> col1_(3) W1: = T * W1 378*> col1_(4) B1: = - V2 * W1 = - B1 * W1 379*> if( IDENT != 'I' ) then 380*> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1 381*> col1_(6_a) below_diag_of_(A1): = - below_diag_of_(W1) 382*> end if 383*> col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1) 384*> 385*> end if 386*> 387*> \endverbatim 388*> 389* ===================================================================== 390 SUBROUTINE ZLARFB_GETT( IDENT, M, N, K, T, LDT, A, LDA, B, LDB, 391 $ WORK, LDWORK ) 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* 597 END 598