1 *> \brief \b CLARFB 2 * 3 * =========== DOCUMENTATION =========== 4 * 5 * Online html documentation available at 6 * http://www.netlib.org/lapack/explore-html/ 7 * 8 *> \htmlonly 9 *> Download CLARFB + dependencies 10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfb.f"> 11 *> [TGZ]</a> 12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfb.f"> 13 *> [ZIP]</a> 14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfb.f"> 15 *> [TXT]</a> 16 *> \endhtmlonly 17 * 18 * Definition: 19 * =========== 20 * 21 * SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, 22 * T, LDT, C, LDC, WORK, LDWORK ) 23 * 24 * .. Scalar Arguments .. 25 * CHARACTER DIRECT, SIDE, STOREV, TRANS 26 * INTEGER K, LDC, LDT, LDV, LDWORK, M, N 27 * .. 28 * .. Array Arguments .. 29 * COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), 30 * $ WORK( LDWORK, * ) 31 * .. 32 * 33 * 34 *> \par Purpose: 35 * ============= 36 *> 37 *> \verbatim 38 *> 39 *> CLARFB applies a complex block reflector H or its transpose H**H to a 40 *> complex M-by-N matrix C, from either the left or the right. 41 *> \endverbatim 42 * 43 * Arguments: 44 * ========== 45 * 46 *> \param[in] SIDE 47 *> \verbatim 48 *> SIDE is CHARACTER*1 49 *> = 'L': apply H or H**H from the Left 50 *> = 'R': apply H or H**H from the Right 51 *> \endverbatim 52 *> 53 *> \param[in] TRANS 54 *> \verbatim 55 *> TRANS is CHARACTER*1 56 *> = 'N': apply H (No transpose) 57 *> = 'C': apply H**H (Conjugate transpose) 58 *> \endverbatim 59 *> 60 *> \param[in] DIRECT 61 *> \verbatim 62 *> DIRECT is CHARACTER*1 63 *> Indicates how H is formed from a product of elementary 64 *> reflectors 65 *> = 'F': H = H(1) H(2) . . . H(k) (Forward) 66 *> = 'B': H = H(k) . . . H(2) H(1) (Backward) 67 *> \endverbatim 68 *> 69 *> \param[in] STOREV 70 *> \verbatim 71 *> STOREV is CHARACTER*1 72 *> Indicates how the vectors which define the elementary 73 *> reflectors are stored: 74 *> = 'C': Columnwise 75 *> = 'R': Rowwise 76 *> \endverbatim 77 *> 78 *> \param[in] M 79 *> \verbatim 80 *> M is INTEGER 81 *> The number of rows of the matrix C. 82 *> \endverbatim 83 *> 84 *> \param[in] N 85 *> \verbatim 86 *> N is INTEGER 87 *> The number of columns of the matrix C. 88 *> \endverbatim 89 *> 90 *> \param[in] K 91 *> \verbatim 92 *> K is INTEGER 93 *> The order of the matrix T (= the number of elementary 94 *> reflectors whose product defines the block reflector). 95 *> \endverbatim 96 *> 97 *> \param[in] V 98 *> \verbatim 99 *> V is COMPLEX array, dimension 100 *> (LDV,K) if STOREV = 'C' 101 *> (LDV,M) if STOREV = 'R' and SIDE = 'L' 102 *> (LDV,N) if STOREV = 'R' and SIDE = 'R' 103 *> The matrix V. See Further Details. 104 *> \endverbatim 105 *> 106 *> \param[in] LDV 107 *> \verbatim 108 *> LDV is INTEGER 109 *> The leading dimension of the array V. 110 *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); 111 *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); 112 *> if STOREV = 'R', LDV >= K. 113 *> \endverbatim 114 *> 115 *> \param[in] T 116 *> \verbatim 117 *> T is COMPLEX array, dimension (LDT,K) 118 *> The triangular K-by-K matrix T in the representation of the 119 *> block reflector. 120 *> \endverbatim 121 *> 122 *> \param[in] LDT 123 *> \verbatim 124 *> LDT is INTEGER 125 *> The leading dimension of the array T. LDT >= K. 126 *> \endverbatim 127 *> 128 *> \param[in,out] C 129 *> \verbatim 130 *> C is COMPLEX array, dimension (LDC,N) 131 *> On entry, the M-by-N matrix C. 132 *> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H. 133 *> \endverbatim 134 *> 135 *> \param[in] LDC 136 *> \verbatim 137 *> LDC is INTEGER 138 *> The leading dimension of the array C. LDC >= max(1,M). 139 *> \endverbatim 140 *> 141 *> \param[out] WORK 142 *> \verbatim 143 *> WORK is COMPLEX array, dimension (LDWORK,K) 144 *> \endverbatim 145 *> 146 *> \param[in] LDWORK 147 *> \verbatim 148 *> LDWORK is INTEGER 149 *> The leading dimension of the array WORK. 150 *> If SIDE = 'L', LDWORK >= max(1,N); 151 *> if SIDE = 'R', LDWORK >= max(1,M). 152 *> \endverbatim 153 * 154 * Authors: 155 * ======== 156 * 157 *> \author Univ. of Tennessee 158 *> \author Univ. of California Berkeley 159 *> \author Univ. of Colorado Denver 160 *> \author NAG Ltd. 161 * 162 *> \date November 2011 163 * 164 *> \ingroup complexOTHERauxiliary 165 * 166 *> \par Further Details: 167 * ===================== 168 *> 169 *> \verbatim 170 *> 171 *> The shape of the matrix V and the storage of the vectors which define 172 *> the H(i) is best illustrated by the following example with n = 5 and 173 *> k = 3. The elements equal to 1 are not stored; the corresponding 174 *> array elements are modified but restored on exit. The rest of the 175 *> array is not used. 176 *> 177 *> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 178 *> 179 *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 180 *> ( v1 1 ) ( 1 v2 v2 v2 ) 181 *> ( v1 v2 1 ) ( 1 v3 v3 ) 182 *> ( v1 v2 v3 ) 183 *> ( v1 v2 v3 ) 184 *> 185 *> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 186 *> 187 *> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 188 *> ( v1 v2 v3 ) ( v2 v2 v2 1 ) 189 *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 190 *> ( 1 v3 ) 191 *> ( 1 ) 192 *> \endverbatim 193 *> 194 * ===================================================================== 195 SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, 196 $ T, LDT, C, LDC, WORK, LDWORK ) 197 * 198 * -- LAPACK auxiliary routine (version 3.4.0) -- 199 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 201 * November 2011 202 * 203 * .. Scalar Arguments .. 204 CHARACTER DIRECT, SIDE, STOREV, TRANS 205 INTEGER K, LDC, LDT, LDV, LDWORK, M, N 206 * .. 207 * .. Array Arguments .. 208 COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), 209 $ WORK( LDWORK, * ) 210 * .. 211 * 212 * ===================================================================== 213 * 214 * .. Parameters .. 215 COMPLEX ONE 216 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 217 * .. 218 * .. Local Scalars .. 219 CHARACTER TRANST 220 INTEGER I, J, LASTV, LASTC 221 * .. 222 * .. External Functions .. 223 LOGICAL LSAME 224 INTEGER ILACLR, ILACLC 225 EXTERNAL LSAME, ILACLR, ILACLC 226 * .. 227 * .. External Subroutines .. 228 EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM 229 * .. 230 * .. Intrinsic Functions .. 231 INTRINSIC CONJG 232 * .. 233 * .. Executable Statements .. 234 * 235 * Quick return if possible 236 * 237 IF( M.LE.0 .OR. N.LE.0 ) 238 $ RETURN 239 * 240 IF( LSAME( TRANS, 'N' ) ) THEN 241 TRANST = 'C' 242 ELSE 243 TRANST = 'N' 244 END IF 245 * 246 IF( LSAME( STOREV, 'C' ) ) THEN 247 * 248 IF( LSAME( DIRECT, 'F' ) ) THEN 249 * 250 * Let V = ( V1 ) (first K rows) 251 * ( V2 ) 252 * where V1 is unit lower triangular. 253 * 254 IF( LSAME( SIDE, 'L' ) ) THEN 255 * 256 * Form H * C or H**H * C where C = ( C1 ) 257 * ( C2 ) 258 * 259 LASTV = MAX( K, ILACLR( M, K, V, LDV ) ) 260 LASTC = ILACLC( LASTV, N, C, LDC ) 261 * 262 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) 263 * 264 * W := C1**H 265 * 266 DO 10 J = 1, K 267 CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 268 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 269 10 CONTINUE 270 * 271 * W := W * V1 272 * 273 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 274 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 275 IF( LASTV.GT.K ) THEN 276 * 277 * W := W + C2**H *V2 278 * 279 CALL CGEMM( 'Conjugate transpose', 'No transpose', 280 $ LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC, 281 $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK ) 282 END IF 283 * 284 * W := W * T**H or W * T 285 * 286 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 287 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 288 * 289 * C := C - V * W**H 290 * 291 IF( M.GT.K ) THEN 292 * 293 * C2 := C2 - V2 * W**H 294 * 295 CALL CGEMM( 'No transpose', 'Conjugate transpose', 296 $ LASTV-K, LASTC, K, -ONE, V( K+1, 1 ), LDV, 297 $ WORK, LDWORK, ONE, C( K+1, 1 ), LDC ) 298 END IF 299 * 300 * W := W * V1**H 301 * 302 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 303 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 304 * 305 * C1 := C1 - W**H 306 * 307 DO 30 J = 1, K 308 DO 20 I = 1, LASTC 309 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) ) 310 20 CONTINUE 311 30 CONTINUE 312 * 313 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 314 * 315 * Form C * H or C * H**H where C = ( C1 C2 ) 316 * 317 LASTV = MAX( K, ILACLR( N, K, V, LDV ) ) 318 LASTC = ILACLR( M, LASTV, C, LDC ) 319 * 320 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 321 * 322 * W := C1 323 * 324 DO 40 J = 1, K 325 CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 326 40 CONTINUE 327 * 328 * W := W * V1 329 * 330 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 331 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 332 IF( LASTV.GT.K ) THEN 333 * 334 * W := W + C2 * V2 335 * 336 CALL CGEMM( 'No transpose', 'No transpose', 337 $ LASTC, K, LASTV-K, 338 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, 339 $ ONE, WORK, LDWORK ) 340 END IF 341 * 342 * W := W * T or W * T**H 343 * 344 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 345 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 346 * 347 * C := C - W * V**H 348 * 349 IF( LASTV.GT.K ) THEN 350 * 351 * C2 := C2 - W * V2**H 352 * 353 CALL CGEMM( 'No transpose', 'Conjugate transpose', 354 $ LASTC, LASTV-K, K, 355 $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, 356 $ ONE, C( 1, K+1 ), LDC ) 357 END IF 358 * 359 * W := W * V1**H 360 * 361 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 362 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 363 * 364 * C1 := C1 - W 365 * 366 DO 60 J = 1, K 367 DO 50 I = 1, LASTC 368 C( I, J ) = C( I, J ) - WORK( I, J ) 369 50 CONTINUE 370 60 CONTINUE 371 END IF 372 * 373 ELSE 374 * 375 * Let V = ( V1 ) 376 * ( V2 ) (last K rows) 377 * where V2 is unit upper triangular. 378 * 379 IF( LSAME( SIDE, 'L' ) ) THEN 380 * 381 * Form H * C or H**H * C where C = ( C1 ) 382 * ( C2 ) 383 * 384 LASTV = MAX( K, ILACLR( M, K, V, LDV ) ) 385 LASTC = ILACLC( LASTV, N, C, LDC ) 386 * 387 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) 388 * 389 * W := C2**H 390 * 391 DO 70 J = 1, K 392 CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 393 $ WORK( 1, J ), 1 ) 394 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 395 70 CONTINUE 396 * 397 * W := W * V2 398 * 399 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 400 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 401 $ WORK, LDWORK ) 402 IF( LASTV.GT.K ) THEN 403 * 404 * W := W + C1**H*V1 405 * 406 CALL CGEMM( 'Conjugate transpose', 'No transpose', 407 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 408 $ ONE, WORK, LDWORK ) 409 END IF 410 * 411 * W := W * T**H or W * T 412 * 413 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 414 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 415 * 416 * C := C - V * W**H 417 * 418 IF( LASTV.GT.K ) THEN 419 * 420 * C1 := C1 - V1 * W**H 421 * 422 CALL CGEMM( 'No transpose', 'Conjugate transpose', 423 $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 424 $ ONE, C, LDC ) 425 END IF 426 * 427 * W := W * V2**H 428 * 429 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 430 $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 431 $ WORK, LDWORK ) 432 * 433 * C2 := C2 - W**H 434 * 435 DO 90 J = 1, K 436 DO 80 I = 1, LASTC 437 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - 438 $ CONJG( WORK( I, J ) ) 439 80 CONTINUE 440 90 CONTINUE 441 * 442 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 443 * 444 * Form C * H or C * H**H where C = ( C1 C2 ) 445 * 446 LASTV = MAX( K, ILACLR( N, K, V, LDV ) ) 447 LASTC = ILACLR( M, LASTV, C, LDC ) 448 * 449 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 450 * 451 * W := C2 452 * 453 DO 100 J = 1, K 454 CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1, 455 $ WORK( 1, J ), 1 ) 456 100 CONTINUE 457 * 458 * W := W * V2 459 * 460 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 461 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 462 $ WORK, LDWORK ) 463 IF( LASTV.GT.K ) THEN 464 * 465 * W := W + C1 * V1 466 * 467 CALL CGEMM( 'No transpose', 'No transpose', 468 $ LASTC, K, LASTV-K, 469 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 470 END IF 471 * 472 * W := W * T or W * T**H 473 * 474 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 475 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 476 * 477 * C := C - W * V**H 478 * 479 IF( LASTV.GT.K ) THEN 480 * 481 * C1 := C1 - W * V1**H 482 * 483 CALL CGEMM( 'No transpose', 'Conjugate transpose', 484 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 485 $ ONE, C, LDC ) 486 END IF 487 * 488 * W := W * V2**H 489 * 490 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 491 $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 492 $ WORK, LDWORK ) 493 * 494 * C2 := C2 - W 495 * 496 DO 120 J = 1, K 497 DO 110 I = 1, LASTC 498 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) 499 $ - WORK( I, J ) 500 110 CONTINUE 501 120 CONTINUE 502 END IF 503 END IF 504 * 505 ELSE IF( LSAME( STOREV, 'R' ) ) THEN 506 * 507 IF( LSAME( DIRECT, 'F' ) ) THEN 508 * 509 * Let V = ( V1 V2 ) (V1: first K columns) 510 * where V1 is unit upper triangular. 511 * 512 IF( LSAME( SIDE, 'L' ) ) THEN 513 * 514 * Form H * C or H**H * C where C = ( C1 ) 515 * ( C2 ) 516 * 517 LASTV = MAX( K, ILACLC( K, M, V, LDV ) ) 518 LASTC = ILACLC( LASTV, N, C, LDC ) 519 * 520 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) 521 * 522 * W := C1**H 523 * 524 DO 130 J = 1, K 525 CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 526 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 527 130 CONTINUE 528 * 529 * W := W * V1**H 530 * 531 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 532 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 533 IF( LASTV.GT.K ) THEN 534 * 535 * W := W + C2**H*V2**H 536 * 537 CALL CGEMM( 'Conjugate transpose', 538 $ 'Conjugate transpose', LASTC, K, LASTV-K, 539 $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, 540 $ ONE, WORK, LDWORK ) 541 END IF 542 * 543 * W := W * T**H or W * T 544 * 545 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 546 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 547 * 548 * C := C - V**H * W**H 549 * 550 IF( LASTV.GT.K ) THEN 551 * 552 * C2 := C2 - V2**H * W**H 553 * 554 CALL CGEMM( 'Conjugate transpose', 555 $ 'Conjugate transpose', LASTV-K, LASTC, K, 556 $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, 557 $ ONE, C( K+1, 1 ), LDC ) 558 END IF 559 * 560 * W := W * V1 561 * 562 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 563 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 564 * 565 * C1 := C1 - W**H 566 * 567 DO 150 J = 1, K 568 DO 140 I = 1, LASTC 569 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) ) 570 140 CONTINUE 571 150 CONTINUE 572 * 573 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 574 * 575 * Form C * H or C * H**H where C = ( C1 C2 ) 576 * 577 LASTV = MAX( K, ILACLC( K, N, V, LDV ) ) 578 LASTC = ILACLR( M, LASTV, C, LDC ) 579 * 580 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) 581 * 582 * W := C1 583 * 584 DO 160 J = 1, K 585 CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 586 160 CONTINUE 587 * 588 * W := W * V1**H 589 * 590 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 591 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 592 IF( LASTV.GT.K ) THEN 593 * 594 * W := W + C2 * V2**H 595 * 596 CALL CGEMM( 'No transpose', 'Conjugate transpose', 597 $ LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC, 598 $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK ) 599 END IF 600 * 601 * W := W * T or W * T**H 602 * 603 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 604 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 605 * 606 * C := C - W * V 607 * 608 IF( LASTV.GT.K ) THEN 609 * 610 * C2 := C2 - W * V2 611 * 612 CALL CGEMM( 'No transpose', 'No transpose', 613 $ LASTC, LASTV-K, K, 614 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, 615 $ ONE, C( 1, K+1 ), LDC ) 616 END IF 617 * 618 * W := W * V1 619 * 620 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 621 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 622 * 623 * C1 := C1 - W 624 * 625 DO 180 J = 1, K 626 DO 170 I = 1, LASTC 627 C( I, J ) = C( I, J ) - WORK( I, J ) 628 170 CONTINUE 629 180 CONTINUE 630 * 631 END IF 632 * 633 ELSE 634 * 635 * Let V = ( V1 V2 ) (V2: last K columns) 636 * where V2 is unit lower triangular. 637 * 638 IF( LSAME( SIDE, 'L' ) ) THEN 639 * 640 * Form H * C or H**H * C where C = ( C1 ) 641 * ( C2 ) 642 * 643 LASTV = MAX( K, ILACLC( K, M, V, LDV ) ) 644 LASTC = ILACLC( LASTV, N, C, LDC ) 645 * 646 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) 647 * 648 * W := C2**H 649 * 650 DO 190 J = 1, K 651 CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 652 $ WORK( 1, J ), 1 ) 653 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 654 190 CONTINUE 655 * 656 * W := W * V2**H 657 * 658 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 659 $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 660 $ WORK, LDWORK ) 661 IF( LASTV.GT.K ) THEN 662 * 663 * W := W + C1**H * V1**H 664 * 665 CALL CGEMM( 'Conjugate transpose', 666 $ 'Conjugate transpose', LASTC, K, LASTV-K, 667 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 668 END IF 669 * 670 * W := W * T**H or W * T 671 * 672 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 673 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 674 * 675 * C := C - V**H * W**H 676 * 677 IF( LASTV.GT.K ) THEN 678 * 679 * C1 := C1 - V1**H * W**H 680 * 681 CALL CGEMM( 'Conjugate transpose', 682 $ 'Conjugate transpose', LASTV-K, LASTC, K, 683 $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC ) 684 END IF 685 * 686 * W := W * V2 687 * 688 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 689 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 690 $ WORK, LDWORK ) 691 * 692 * C2 := C2 - W**H 693 * 694 DO 210 J = 1, K 695 DO 200 I = 1, LASTC 696 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - 697 $ CONJG( WORK( I, J ) ) 698 200 CONTINUE 699 210 CONTINUE 700 * 701 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 702 * 703 * Form C * H or C * H**H where C = ( C1 C2 ) 704 * 705 LASTV = MAX( K, ILACLC( K, N, V, LDV ) ) 706 LASTC = ILACLR( M, LASTV, C, LDC ) 707 * 708 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) 709 * 710 * W := C2 711 * 712 DO 220 J = 1, K 713 CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1, 714 $ WORK( 1, J ), 1 ) 715 220 CONTINUE 716 * 717 * W := W * V2**H 718 * 719 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 720 $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 721 $ WORK, LDWORK ) 722 IF( LASTV.GT.K ) THEN 723 * 724 * W := W + C1 * V1**H 725 * 726 CALL CGEMM( 'No transpose', 'Conjugate transpose', 727 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE, 728 $ WORK, LDWORK ) 729 END IF 730 * 731 * W := W * T or W * T**H 732 * 733 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 734 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 735 * 736 * C := C - W * V 737 * 738 IF( LASTV.GT.K ) THEN 739 * 740 * C1 := C1 - W * V1 741 * 742 CALL CGEMM( 'No transpose', 'No transpose', 743 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 744 $ ONE, C, LDC ) 745 END IF 746 * 747 * W := W * V2 748 * 749 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 750 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 751 $ WORK, LDWORK ) 752 * 753 * C1 := C1 - W 754 * 755 DO 240 J = 1, K 756 DO 230 I = 1, LASTC 757 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) 758 $ - WORK( I, J ) 759 230 CONTINUE 760 240 CONTINUE 761 * 762 END IF 763 * 764 END IF 765 END IF 766 * 767 RETURN 768 * 769 * End of CLARFB 770 * 771 END 772
