SRC/stprfb.f

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