src/LAPACK/dlarfb.f

   1       SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
   2      $                   T, LDT, C, LDC, WORK, LDWORK )
   3 *
   4 *  -- LAPACK auxiliary routine (version 3.1) --
   5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   6 *     November 2006
   7 *
   8 *     .. Scalar Arguments ..
   9       CHARACTER          DIRECT, SIDE, STOREV, TRANS
  10       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
  11 *     ..
  12 *     .. Array Arguments ..
  13       DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
  14      $                   WORK( LDWORK, * )
  15 *     ..
  16 *
  17 *  Purpose
  18 *  =======
  19 *
  20 *  DLARFB applies a real block reflector H or its transpose H' to a
  21 *  real m by n matrix C, from either the left or the right.
  22 *
  23 *  Arguments
  24 *  =========
  25 *
  26 *  SIDE    (input) CHARACTER*1
  27 *          = 'L': apply H or H' from the Left
  28 *          = 'R': apply H or H' from the Right
  29 *
  30 *  TRANS   (input) CHARACTER*1
  31 *          = 'N': apply H (No transpose)
  32 *          = 'T': apply H' (Transpose)
  33 *
  34 *  DIRECT  (input) CHARACTER*1
  35 *          Indicates how H is formed from a product of elementary
  36 *          reflectors
  37 *          = 'F': H = H(1) H(2) . . . H(k) (Forward)
  38 *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
  39 *
  40 *  STOREV  (input) CHARACTER*1
  41 *          Indicates how the vectors which define the elementary
  42 *          reflectors are stored:
  43 *          = 'C': Columnwise
  44 *          = 'R': Rowwise
  45 *
  46 *  M       (input) INTEGER
  47 *          The number of rows of the matrix C.
  48 *
  49 *  N       (input) INTEGER
  50 *          The number of columns of the matrix C.
  51 *
  52 *  K       (input) INTEGER
  53 *          The order of the matrix T (= the number of elementary
  54 *          reflectors whose product defines the block reflector).
  55 *
  56 *  V       (input) DOUBLE PRECISION array, dimension
  57 *                                (LDV,K) if STOREV = 'C'
  58 *                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
  59 *                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
  60 *          The matrix V. See further details.
  61 *
  62 *  LDV     (input) INTEGER
  63 *          The leading dimension of the array V.
  64 *          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
  65 *          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
  66 *          if STOREV = 'R', LDV >= K.
  67 *
  68 *  T       (input) DOUBLE PRECISION array, dimension (LDT,K)
  69 *          The triangular k by k matrix T in the representation of the
  70 *          block reflector.
  71 *
  72 *  LDT     (input) INTEGER
  73 *          The leading dimension of the array T. LDT >= K.
  74 *
  75 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
  76 *          On entry, the m by n matrix C.
  77 *          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
  78 *
  79 *  LDC     (input) INTEGER
  80 *          The leading dimension of the array C. LDA >= max(1,M).
  81 *
  82 *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
  83 *
  84 *  LDWORK  (input) INTEGER
  85 *          The leading dimension of the array WORK.
  86 *          If SIDE = 'L', LDWORK >= max(1,N);
  87 *          if SIDE = 'R', LDWORK >= max(1,M).
  88 *
  89 *  =====================================================================
  90 *
  91 *     .. Parameters ..
  92       DOUBLE PRECISION   ONE
  93       PARAMETER          ( ONE = 1.0D+0 )
  94 *     ..
  95 *     .. Local Scalars ..
  96       CHARACTER          TRANST
  97       INTEGER            I, J
  98 *     ..
  99 *     .. External Functions ..
 100       LOGICAL            LSAME
 101       EXTERNAL           LSAME
 102 *     ..
 103 *     .. External Subroutines ..
 104       EXTERNAL           DCOPY, DGEMM, DTRMM
 105 *     ..
 106 *     .. Executable Statements ..
 107 *
 108 *     Quick return if possible
 109 *
 110       IF( M.LE.0 .OR. N.LE.0 )
 111      $   RETURN
 112 *
 113       IF( LSAME( TRANS, 'N' ) ) THEN
 114          TRANST = 'T'
 115       ELSE
 116          TRANST = 'N'
 117       END IF
 118 *
 119       IF( LSAME( STOREV, 'C' ) ) THEN
 120 *
 121          IF( LSAME( DIRECT, 'F' ) ) THEN
 122 *
 123 *           Let  V =  ( V1 )    (first K rows)
 124 *                     ( V2 )
 125 *           where  V1  is unit lower triangular.
 126 *
 127             IF( LSAME( SIDE, 'L' ) ) THEN
 128 *
 129 *              Form  H * C  or  H' * C  where  C = ( C1 )
 130 *                                                  ( C2 )
 131 *
 132 *              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
 133 *
 134 *              W := C1'
 135 *
 136                DO 10 J = 1, K
 137                   CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 138    10          CONTINUE
 139 *
 140 *              W := W * V1
 141 *
 142                CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
 143      $                     K, ONE, V, LDV, WORK, LDWORK )
 144                IF( M.GT.K ) THEN
 145 *
 146 *                 W := W + C2'*V2
 147 *
 148                   CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K,
 149      $                        ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
 150      $                        ONE, WORK, LDWORK )
 151                END IF
 152 *
 153 *              W := W * T'  or  W * T
 154 *
 155                CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
 156      $                     ONE, T, LDT, WORK, LDWORK )
 157 *
 158 *              C := C - V * W'
 159 *
 160                IF( M.GT.K ) THEN
 161 *
 162 *                 C2 := C2 - V2 * W'
 163 *
 164                   CALL DGEMM( 'No transpose', 'Transpose', M-K, N, K,
 165      $                        -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
 166      $                        C( K+1, 1 ), LDC )
 167                END IF
 168 *
 169 *              W := W * V1'
 170 *
 171                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K,
 172      $                     ONE, V, LDV, WORK, LDWORK )
 173 *
 174 *              C1 := C1 - W'
 175 *
 176                DO 30 J = 1, K
 177                   DO 20 I = 1, N
 178                      C( J, I ) = C( J, I ) - WORK( I, J )
 179    20             CONTINUE
 180    30          CONTINUE
 181 *
 182             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 183 *
 184 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 185 *
 186 *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
 187 *
 188 *              W := C1
 189 *
 190                DO 40 J = 1, K
 191                   CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 192    40          CONTINUE
 193 *
 194 *              W := W * V1
 195 *
 196                CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
 197      $                     K, ONE, V, LDV, WORK, LDWORK )
 198                IF( N.GT.K ) THEN
 199 *
 200 *                 W := W + C2 * V2
 201 *
 202                   CALL DGEMM( 'No transpose', 'No transpose', M, K, N-K,
 203      $                        ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
 204      $                        ONE, WORK, LDWORK )
 205                END IF
 206 *
 207 *              W := W * T  or  W * T'
 208 *
 209                CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
 210      $                     ONE, T, LDT, WORK, LDWORK )
 211 *
 212 *              C := C - W * V'
 213 *
 214                IF( N.GT.K ) THEN
 215 *
 216 *                 C2 := C2 - W * V2'
 217 *
 218                   CALL DGEMM( 'No transpose', 'Transpose', M, N-K, K,
 219      $                        -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
 220      $                        C( 1, K+1 ), LDC )
 221                END IF
 222 *
 223 *              W := W * V1'
 224 *
 225                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K,
 226      $                     ONE, V, LDV, WORK, LDWORK )
 227 *
 228 *              C1 := C1 - W
 229 *
 230                DO 60 J = 1, K
 231                   DO 50 I = 1, M
 232                      C( I, J ) = C( I, J ) - WORK( I, J )
 233    50             CONTINUE
 234    60          CONTINUE
 235             END IF
 236 *
 237          ELSE
 238 *
 239 *           Let  V =  ( V1 )
 240 *                     ( V2 )    (last K rows)
 241 *           where  V2  is unit upper triangular.
 242 *
 243             IF( LSAME( SIDE, 'L' ) ) THEN
 244 *
 245 *              Form  H * C  or  H' * C  where  C = ( C1 )
 246 *                                                  ( C2 )
 247 *
 248 *              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
 249 *
 250 *              W := C2'
 251 *
 252                DO 70 J = 1, K
 253                   CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
 254    70          CONTINUE
 255 *
 256 *              W := W * V2
 257 *
 258                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
 259      $                     K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
 260                IF( M.GT.K ) THEN
 261 *
 262 *                 W := W + C1'*V1
 263 *
 264                   CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K,
 265      $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
 266                END IF
 267 *
 268 *              W := W * T'  or  W * T
 269 *
 270                CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
 271      $                     ONE, T, LDT, WORK, LDWORK )
 272 *
 273 *              C := C - V * W'
 274 *
 275                IF( M.GT.K ) THEN
 276 *
 277 *                 C1 := C1 - V1 * W'
 278 *
 279                   CALL DGEMM( 'No transpose', 'Transpose', M-K, N, K,
 280      $                        -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC )
 281                END IF
 282 *
 283 *              W := W * V2'
 284 *
 285                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K,
 286      $                     ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
 287 *
 288 *              C2 := C2 - W'
 289 *
 290                DO 90 J = 1, K
 291                   DO 80 I = 1, N
 292                      C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J )
 293    80             CONTINUE
 294    90          CONTINUE
 295 *
 296             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 297 *
 298 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 299 *
 300 *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
 301 *
 302 *              W := C2
 303 *
 304                DO 100 J = 1, K
 305                   CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
 306   100          CONTINUE
 307 *
 308 *              W := W * V2
 309 *
 310                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
 311      $                     K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
 312                IF( N.GT.K ) THEN
 313 *
 314 *                 W := W + C1 * V1
 315 *
 316                   CALL DGEMM( 'No transpose', 'No transpose', M, K, N-K,
 317      $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
 318                END IF
 319 *
 320 *              W := W * T  or  W * T'
 321 *
 322                CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
 323      $                     ONE, T, LDT, WORK, LDWORK )
 324 *
 325 *              C := C - W * V'
 326 *
 327                IF( N.GT.K ) THEN
 328 *
 329 *                 C1 := C1 - W * V1'
 330 *
 331                   CALL DGEMM( 'No transpose', 'Transpose', M, N-K, K,
 332      $                        -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
 333                END IF
 334 *
 335 *              W := W * V2'
 336 *
 337                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K,
 338      $                     ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
 339 *
 340 *              C2 := C2 - W
 341 *
 342                DO 120 J = 1, K
 343                   DO 110 I = 1, M
 344                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
 345   110             CONTINUE
 346   120          CONTINUE
 347             END IF
 348          END IF
 349 *
 350       ELSE IF( LSAME( STOREV, 'R' ) ) THEN
 351 *
 352          IF( LSAME( DIRECT, 'F' ) ) THEN
 353 *
 354 *           Let  V =  ( V1  V2 )    (V1: first K columns)
 355 *           where  V1  is unit upper triangular.
 356 *
 357             IF( LSAME( SIDE, 'L' ) ) THEN
 358 *
 359 *              Form  H * C  or  H' * C  where  C = ( C1 )
 360 *                                                  ( C2 )
 361 *
 362 *              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
 363 *
 364 *              W := C1'
 365 *
 366                DO 130 J = 1, K
 367                   CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 368   130          CONTINUE
 369 *
 370 *              W := W * V1'
 371 *
 372                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K,
 373      $                     ONE, V, LDV, WORK, LDWORK )
 374                IF( M.GT.K ) THEN
 375 *
 376 *                 W := W + C2'*V2'
 377 *
 378                   CALL DGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE,
 379      $                        C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
 380      $                        WORK, LDWORK )
 381                END IF
 382 *
 383 *              W := W * T'  or  W * T
 384 *
 385                CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
 386      $                     ONE, T, LDT, WORK, LDWORK )
 387 *
 388 *              C := C - V' * W'
 389 *
 390                IF( M.GT.K ) THEN
 391 *
 392 *                 C2 := C2 - V2' * W'
 393 *
 394                   CALL DGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE,
 395      $                        V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
 396      $                        C( K+1, 1 ), LDC )
 397                END IF
 398 *
 399 *              W := W * V1
 400 *
 401                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
 402      $                     K, ONE, V, LDV, WORK, LDWORK )
 403 *
 404 *              C1 := C1 - W'
 405 *
 406                DO 150 J = 1, K
 407                   DO 140 I = 1, N
 408                      C( J, I ) = C( J, I ) - WORK( I, J )
 409   140             CONTINUE
 410   150          CONTINUE
 411 *
 412             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 413 *
 414 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 415 *
 416 *              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
 417 *
 418 *              W := C1
 419 *
 420                DO 160 J = 1, K
 421                   CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 422   160          CONTINUE
 423 *
 424 *              W := W * V1'
 425 *
 426                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K,
 427      $                     ONE, V, LDV, WORK, LDWORK )
 428                IF( N.GT.K ) THEN
 429 *
 430 *                 W := W + C2 * V2'
 431 *
 432                   CALL DGEMM( 'No transpose', 'Transpose', M, K, N-K,
 433      $                        ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
 434      $                        ONE, WORK, LDWORK )
 435                END IF
 436 *
 437 *              W := W * T  or  W * T'
 438 *
 439                CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
 440      $                     ONE, T, LDT, WORK, LDWORK )
 441 *
 442 *              C := C - W * V
 443 *
 444                IF( N.GT.K ) THEN
 445 *
 446 *                 C2 := C2 - W * V2
 447 *
 448                   CALL DGEMM( 'No transpose', 'No transpose', M, N-K, K,
 449      $                        -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
 450      $                        C( 1, K+1 ), LDC )
 451                END IF
 452 *
 453 *              W := W * V1
 454 *
 455                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
 456      $                     K, ONE, V, LDV, WORK, LDWORK )
 457 *
 458 *              C1 := C1 - W
 459 *
 460                DO 180 J = 1, K
 461                   DO 170 I = 1, M
 462                      C( I, J ) = C( I, J ) - WORK( I, J )
 463   170             CONTINUE
 464   180          CONTINUE
 465 *
 466             END IF
 467 *
 468          ELSE
 469 *
 470 *           Let  V =  ( V1  V2 )    (V2: last K columns)
 471 *           where  V2  is unit lower triangular.
 472 *
 473             IF( LSAME( SIDE, 'L' ) ) THEN
 474 *
 475 *              Form  H * C  or  H' * C  where  C = ( C1 )
 476 *                                                  ( C2 )
 477 *
 478 *              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
 479 *
 480 *              W := C2'
 481 *
 482                DO 190 J = 1, K
 483                   CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
 484   190          CONTINUE
 485 *
 486 *              W := W * V2'
 487 *
 488                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K,
 489      $                     ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
 490                IF( M.GT.K ) THEN
 491 *
 492 *                 W := W + C1'*V1'
 493 *
 494                   CALL DGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE,
 495      $                        C, LDC, V, LDV, ONE, WORK, LDWORK )
 496                END IF
 497 *
 498 *              W := W * T'  or  W * T
 499 *
 500                CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
 501      $                     ONE, T, LDT, WORK, LDWORK )
 502 *
 503 *              C := C - V' * W'
 504 *
 505                IF( M.GT.K ) THEN
 506 *
 507 *                 C1 := C1 - V1' * W'
 508 *
 509                   CALL DGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE,
 510      $                        V, LDV, WORK, LDWORK, ONE, C, LDC )
 511                END IF
 512 *
 513 *              W := W * V2
 514 *
 515                CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
 516      $                     K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
 517 *
 518 *              C2 := C2 - W'
 519 *
 520                DO 210 J = 1, K
 521                   DO 200 I = 1, N
 522                      C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J )
 523   200             CONTINUE
 524   210          CONTINUE
 525 *
 526             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 527 *
 528 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 529 *
 530 *              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
 531 *
 532 *              W := C2
 533 *
 534                DO 220 J = 1, K
 535                   CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
 536   220          CONTINUE
 537 *
 538 *              W := W * V2'
 539 *
 540                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K,
 541      $                     ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
 542                IF( N.GT.K ) THEN
 543 *
 544 *                 W := W + C1 * V1'
 545 *
 546                   CALL DGEMM( 'No transpose', 'Transpose', M, K, N-K,
 547      $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
 548                END IF
 549 *
 550 *              W := W * T  or  W * T'
 551 *
 552                CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
 553      $                     ONE, T, LDT, WORK, LDWORK )
 554 *
 555 *              C := C - W * V
 556 *
 557                IF( N.GT.K ) THEN
 558 *
 559 *                 C1 := C1 - W * V1
 560 *
 561                   CALL DGEMM( 'No transpose', 'No transpose', M, N-K, K,
 562      $                        -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
 563                END IF
 564 *
 565 *              W := W * V2
 566 *
 567                CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
 568      $                     K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
 569 *
 570 *              C1 := C1 - W
 571 *
 572                DO 240 J = 1, K
 573                   DO 230 I = 1, M
 574                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
 575   230             CONTINUE
 576   240          CONTINUE
 577 *
 578             END IF
 579 *
 580          END IF
 581       END IF
 582 *
 583       RETURN
 584 *
 585 *     End of DLARFB
 586 *
 587       END