1 *> \brief \b DLARFB 2 * 3 * =========== DOCUMENTATION =========== 4 * 5 * Online html documentation available at 6 * http://www.netlib.org/lapack/explore-html/ 7 * 8 *> \htmlonly 9 *> Download DLARFB + dependencies 10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb.f"> 11 *> [TGZ]</a> 12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb.f"> 13 *> [ZIP]</a> 14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb.f"> 15 *> [TXT]</a> 16 *> \endhtmlonly 17 * 18 * Definition: 19 * =========== 20 * 21 * SUBROUTINE DLARFB( 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 * DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), 30 * $ WORK( LDWORK, * ) 31 * .. 32 * 33 * 34 *> \par Purpose: 35 * ============= 36 *> 37 *> \verbatim 38 *> 39 *> DLARFB applies a real block reflector H or its transpose H**T to a 40 *> real 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**T from the Left 50 *> = 'R': apply H or H**T from the Right 51 *> \endverbatim 52 *> 53 *> \param[in] TRANS 54 *> \verbatim 55 *> TRANS is CHARACTER*1 56 *> = 'N': apply H (No transpose) 57 *> = 'T': apply H**T (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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDC,N) 131 *> On entry, the m by n matrix C. 132 *> On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T. 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 DOUBLE PRECISION 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 doubleOTHERauxiliary 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 DLARFB( 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 DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), 209 $ WORK( LDWORK, * ) 210 * .. 211 * 212 * ===================================================================== 213 * 214 * .. Parameters .. 215 DOUBLE PRECISION ONE 216 PARAMETER ( ONE = 1.0D+0 ) 217 * .. 218 * .. Local Scalars .. 219 CHARACTER TRANST 220 INTEGER I, J, LASTV, LASTC 221 * .. 222 * .. External Functions .. 223 LOGICAL LSAME 224 INTEGER ILADLR, ILADLC 225 EXTERNAL LSAME, ILADLR, ILADLC 226 * .. 227 * .. External Subroutines .. 228 EXTERNAL DCOPY, DGEMM, DTRMM 229 * .. 230 * .. Executable Statements .. 231 * 232 * Quick return if possible 233 * 234 IF( M.LE.0 .OR. N.LE.0 ) 235 $ RETURN 236 * 237 IF( LSAME( TRANS, 'N' ) ) THEN 238 TRANST = 'T' 239 ELSE 240 TRANST = 'N' 241 END IF 242 * 243 IF( LSAME( STOREV, 'C' ) ) THEN 244 * 245 IF( LSAME( DIRECT, 'F' ) ) THEN 246 * 247 * Let V = ( V1 ) (first K rows) 248 * ( V2 ) 249 * where V1 is unit lower triangular. 250 * 251 IF( LSAME( SIDE, 'L' ) ) THEN 252 * 253 * Form H * C or H**T * C where C = ( C1 ) 254 * ( C2 ) 255 * 256 LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) 257 LASTC = ILADLC( LASTV, N, C, LDC ) 258 * 259 * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) 260 * 261 * W := C1**T 262 * 263 DO 10 J = 1, K 264 CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 265 10 CONTINUE 266 * 267 * W := W * V1 268 * 269 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 270 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 271 IF( LASTV.GT.K ) THEN 272 * 273 * W := W + C2**T *V2 274 * 275 CALL DGEMM( 'Transpose', 'No transpose', 276 $ LASTC, K, LASTV-K, 277 $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV, 278 $ ONE, WORK, LDWORK ) 279 END IF 280 * 281 * W := W * T**T or W * T 282 * 283 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 284 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 285 * 286 * C := C - V * W**T 287 * 288 IF( LASTV.GT.K ) THEN 289 * 290 * C2 := C2 - V2 * W**T 291 * 292 CALL DGEMM( 'No transpose', 'Transpose', 293 $ LASTV-K, LASTC, K, 294 $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE, 295 $ C( K+1, 1 ), LDC ) 296 END IF 297 * 298 * W := W * V1**T 299 * 300 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 301 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 302 * 303 * C1 := C1 - W**T 304 * 305 DO 30 J = 1, K 306 DO 20 I = 1, LASTC 307 C( J, I ) = C( J, I ) - WORK( I, J ) 308 20 CONTINUE 309 30 CONTINUE 310 * 311 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 312 * 313 * Form C * H or C * H**T where C = ( C1 C2 ) 314 * 315 LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) 316 LASTC = ILADLR( M, LASTV, C, LDC ) 317 * 318 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 319 * 320 * W := C1 321 * 322 DO 40 J = 1, K 323 CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 324 40 CONTINUE 325 * 326 * W := W * V1 327 * 328 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 329 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 330 IF( LASTV.GT.K ) THEN 331 * 332 * W := W + C2 * V2 333 * 334 CALL DGEMM( 'No transpose', 'No transpose', 335 $ LASTC, K, LASTV-K, 336 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, 337 $ ONE, WORK, LDWORK ) 338 END IF 339 * 340 * W := W * T or W * T**T 341 * 342 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 343 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 344 * 345 * C := C - W * V**T 346 * 347 IF( LASTV.GT.K ) THEN 348 * 349 * C2 := C2 - W * V2**T 350 * 351 CALL DGEMM( 'No transpose', 'Transpose', 352 $ LASTC, LASTV-K, K, 353 $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE, 354 $ C( 1, K+1 ), LDC ) 355 END IF 356 * 357 * W := W * V1**T 358 * 359 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 360 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 361 * 362 * C1 := C1 - W 363 * 364 DO 60 J = 1, K 365 DO 50 I = 1, LASTC 366 C( I, J ) = C( I, J ) - WORK( I, J ) 367 50 CONTINUE 368 60 CONTINUE 369 END IF 370 * 371 ELSE 372 * 373 * Let V = ( V1 ) 374 * ( V2 ) (last K rows) 375 * where V2 is unit upper triangular. 376 * 377 IF( LSAME( SIDE, 'L' ) ) THEN 378 * 379 * Form H * C or H**T * C where C = ( C1 ) 380 * ( C2 ) 381 * 382 LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) 383 LASTC = ILADLC( LASTV, N, C, LDC ) 384 * 385 * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) 386 * 387 * W := C2**T 388 * 389 DO 70 J = 1, K 390 CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 391 $ WORK( 1, J ), 1 ) 392 70 CONTINUE 393 * 394 * W := W * V2 395 * 396 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 397 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 398 $ WORK, LDWORK ) 399 IF( LASTV.GT.K ) THEN 400 * 401 * W := W + C1**T*V1 402 * 403 CALL DGEMM( 'Transpose', 'No transpose', 404 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 405 $ ONE, WORK, LDWORK ) 406 END IF 407 * 408 * W := W * T**T or W * T 409 * 410 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 411 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 412 * 413 * C := C - V * W**T 414 * 415 IF( LASTV.GT.K ) THEN 416 * 417 * C1 := C1 - V1 * W**T 418 * 419 CALL DGEMM( 'No transpose', 'Transpose', 420 $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 421 $ ONE, C, LDC ) 422 END IF 423 * 424 * W := W * V2**T 425 * 426 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 427 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 428 $ WORK, LDWORK ) 429 * 430 * C2 := C2 - W**T 431 * 432 DO 90 J = 1, K 433 DO 80 I = 1, LASTC 434 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) 435 80 CONTINUE 436 90 CONTINUE 437 * 438 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 439 * 440 * Form C * H or C * H**T where C = ( C1 C2 ) 441 * 442 LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) 443 LASTC = ILADLR( M, LASTV, C, LDC ) 444 * 445 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 446 * 447 * W := C2 448 * 449 DO 100 J = 1, K 450 CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) 451 100 CONTINUE 452 * 453 * W := W * V2 454 * 455 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 456 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 457 $ WORK, LDWORK ) 458 IF( LASTV.GT.K ) THEN 459 * 460 * W := W + C1 * V1 461 * 462 CALL DGEMM( 'No transpose', 'No transpose', 463 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 464 $ ONE, WORK, LDWORK ) 465 END IF 466 * 467 * W := W * T or W * T**T 468 * 469 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 470 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 471 * 472 * C := C - W * V**T 473 * 474 IF( LASTV.GT.K ) THEN 475 * 476 * C1 := C1 - W * V1**T 477 * 478 CALL DGEMM( 'No transpose', 'Transpose', 479 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 480 $ ONE, C, LDC ) 481 END IF 482 * 483 * W := W * V2**T 484 * 485 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 486 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 487 $ WORK, LDWORK ) 488 * 489 * C2 := C2 - W 490 * 491 DO 120 J = 1, K 492 DO 110 I = 1, LASTC 493 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) 494 110 CONTINUE 495 120 CONTINUE 496 END IF 497 END IF 498 * 499 ELSE IF( LSAME( STOREV, 'R' ) ) THEN 500 * 501 IF( LSAME( DIRECT, 'F' ) ) THEN 502 * 503 * Let V = ( V1 V2 ) (V1: first K columns) 504 * where V1 is unit upper triangular. 505 * 506 IF( LSAME( SIDE, 'L' ) ) THEN 507 * 508 * Form H * C or H**T * C where C = ( C1 ) 509 * ( C2 ) 510 * 511 LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) 512 LASTC = ILADLC( LASTV, N, C, LDC ) 513 * 514 * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) 515 * 516 * W := C1**T 517 * 518 DO 130 J = 1, K 519 CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 520 130 CONTINUE 521 * 522 * W := W * V1**T 523 * 524 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 525 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 526 IF( LASTV.GT.K ) THEN 527 * 528 * W := W + C2**T*V2**T 529 * 530 CALL DGEMM( 'Transpose', 'Transpose', 531 $ LASTC, K, LASTV-K, 532 $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, 533 $ ONE, WORK, LDWORK ) 534 END IF 535 * 536 * W := W * T**T or W * T 537 * 538 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 539 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 540 * 541 * C := C - V**T * W**T 542 * 543 IF( LASTV.GT.K ) THEN 544 * 545 * C2 := C2 - V2**T * W**T 546 * 547 CALL DGEMM( 'Transpose', 'Transpose', 548 $ LASTV-K, LASTC, K, 549 $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, 550 $ ONE, C( K+1, 1 ), LDC ) 551 END IF 552 * 553 * W := W * V1 554 * 555 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 556 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 557 * 558 * C1 := C1 - W**T 559 * 560 DO 150 J = 1, K 561 DO 140 I = 1, LASTC 562 C( J, I ) = C( J, I ) - WORK( I, J ) 563 140 CONTINUE 564 150 CONTINUE 565 * 566 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 567 * 568 * Form C * H or C * H**T where C = ( C1 C2 ) 569 * 570 LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) 571 LASTC = ILADLR( M, LASTV, C, LDC ) 572 * 573 * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) 574 * 575 * W := C1 576 * 577 DO 160 J = 1, K 578 CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 579 160 CONTINUE 580 * 581 * W := W * V1**T 582 * 583 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 584 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 585 IF( LASTV.GT.K ) THEN 586 * 587 * W := W + C2 * V2**T 588 * 589 CALL DGEMM( 'No transpose', 'Transpose', 590 $ LASTC, K, LASTV-K, 591 $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV, 592 $ ONE, WORK, LDWORK ) 593 END IF 594 * 595 * W := W * T or W * T**T 596 * 597 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 598 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 599 * 600 * C := C - W * V 601 * 602 IF( LASTV.GT.K ) THEN 603 * 604 * C2 := C2 - W * V2 605 * 606 CALL DGEMM( 'No transpose', 'No transpose', 607 $ LASTC, LASTV-K, K, 608 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, 609 $ ONE, C( 1, K+1 ), LDC ) 610 END IF 611 * 612 * W := W * V1 613 * 614 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 615 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 616 * 617 * C1 := C1 - W 618 * 619 DO 180 J = 1, K 620 DO 170 I = 1, LASTC 621 C( I, J ) = C( I, J ) - WORK( I, J ) 622 170 CONTINUE 623 180 CONTINUE 624 * 625 END IF 626 * 627 ELSE 628 * 629 * Let V = ( V1 V2 ) (V2: last K columns) 630 * where V2 is unit lower triangular. 631 * 632 IF( LSAME( SIDE, 'L' ) ) THEN 633 * 634 * Form H * C or H**T * C where C = ( C1 ) 635 * ( C2 ) 636 * 637 LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) 638 LASTC = ILADLC( LASTV, N, C, LDC ) 639 * 640 * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) 641 * 642 * W := C2**T 643 * 644 DO 190 J = 1, K 645 CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 646 $ WORK( 1, J ), 1 ) 647 190 CONTINUE 648 * 649 * W := W * V2**T 650 * 651 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 652 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 653 $ WORK, LDWORK ) 654 IF( LASTV.GT.K ) THEN 655 * 656 * W := W + C1**T * V1**T 657 * 658 CALL DGEMM( 'Transpose', 'Transpose', 659 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 660 $ ONE, WORK, LDWORK ) 661 END IF 662 * 663 * W := W * T**T or W * T 664 * 665 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 666 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 667 * 668 * C := C - V**T * W**T 669 * 670 IF( LASTV.GT.K ) THEN 671 * 672 * C1 := C1 - V1**T * W**T 673 * 674 CALL DGEMM( 'Transpose', 'Transpose', 675 $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 676 $ ONE, C, LDC ) 677 END IF 678 * 679 * W := W * V2 680 * 681 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 682 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 683 $ WORK, LDWORK ) 684 * 685 * C2 := C2 - W**T 686 * 687 DO 210 J = 1, K 688 DO 200 I = 1, LASTC 689 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) 690 200 CONTINUE 691 210 CONTINUE 692 * 693 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 694 * 695 * Form C * H or C * H**T where C = ( C1 C2 ) 696 * 697 LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) 698 LASTC = ILADLR( M, LASTV, C, LDC ) 699 * 700 * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) 701 * 702 * W := C2 703 * 704 DO 220 J = 1, K 705 CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1, 706 $ WORK( 1, J ), 1 ) 707 220 CONTINUE 708 * 709 * W := W * V2**T 710 * 711 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 712 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 713 $ WORK, LDWORK ) 714 IF( LASTV.GT.K ) THEN 715 * 716 * W := W + C1 * V1**T 717 * 718 CALL DGEMM( 'No transpose', 'Transpose', 719 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 720 $ ONE, WORK, LDWORK ) 721 END IF 722 * 723 * W := W * T or W * T**T 724 * 725 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 726 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 727 * 728 * C := C - W * V 729 * 730 IF( LASTV.GT.K ) THEN 731 * 732 * C1 := C1 - W * V1 733 * 734 CALL DGEMM( 'No transpose', 'No transpose', 735 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 736 $ ONE, C, LDC ) 737 END IF 738 * 739 * W := W * V2 740 * 741 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 742 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 743 $ WORK, LDWORK ) 744 * 745 * C1 := C1 - W 746 * 747 DO 240 J = 1, K 748 DO 230 I = 1, LASTC 749 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) 750 230 CONTINUE 751 240 CONTINUE 752 * 753 END IF 754 * 755 END IF 756 END IF 757 * 758 RETURN 759 * 760 * End of DLARFB 761 * 762 END 763
