SRC/slarfb.f

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