QR factorization for triangular-pentagonal matrices

(generalization of triangle-over-square and triangle-over-triangle)

1 files changed, 705 insertions, 0 deletions

diff --git a/SRC/stprfb.f b/SRC/stprfb.f
new file mode 100644
index 00000000..234b8790
--- /dev/null
+++ b/SRC/stprfb.f
@@ -0,0 +1,705 @@
+      SUBROUTINE STPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, 
+     $                   V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
+      IMPLICIT NONE
+*
+*  -- LAPACK auxiliary routine (version 3.x) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*  -- July 2011                                                      --
+*
+*     .. Scalar Arguments ..
+      CHARACTER DIRECT, SIDE, STOREV, TRANS
+      INTEGER   K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
+*     ..
+*     .. Array Arguments ..
+      REAL   A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          V( LDV, * ), WORK( LDWORK, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  STPRFB applies a real "triangular-pentagonal" block reflector H or its 
+*  conjugate transpose H^H to a real matrix C, which is composed of two 
+*  blocks A and B, either from the left or right.
+*  
+*  Arguments
+*  =========
+*
+*  SIDE    (input) CHARACTER*1
+*          = 'L': apply H or H^H from the Left
+*          = 'R': apply H or H^H from the Right
+*
+*  TRANS   (input) CHARACTER*1
+*          = 'N': apply H (No transpose)
+*          = 'C': apply H^H (Conjugate transpose)
+*
+*  DIRECT  (input) CHARACTER*1
+*          Indicates how H is formed from a product of elementary
+*          reflectors
+*          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*
+*  STOREV  (input) CHARACTER*1
+*          Indicates how the vectors which define the elementary
+*          reflectors are stored:
+*          = 'C': Columns
+*          = 'R': Rows
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B.  
+*          N >= 0.
+*
+*  K       (input) INTEGER
+*          The order of the matrix T, i.e. the number of elementary
+*          reflectors whose product defines the block reflector.  
+*          K >= 0.
+*
+*  L       (input) INTEGER
+*          The order of the trapezoidal part of V.  
+*          K >= L >= 0.  See Further Details.
+*
+*  V       (input) REAL array, dimension
+*                                (LDV,K) if STOREV = 'C'
+*                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
+*                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
+*          The pentagonal matrix V, which contains the elementary reflectors
+*          H(1), H(2), ..., H(K).  See Further Details.
+*
+*  LDV     (input) INTEGER
+*          The leading dimension of the array V.
+*          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
+*          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
+*          if STOREV = 'R', LDV >= K.
+*
+*  T       (input) REAL array, dimension (LDT,K)
+*          The triangular K-by-K matrix T in the representation of the
+*          block reflector.  
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T. 
+*          LDT >= K.
+*
+*  A       (input/output) REAL array, dimension 
+*          (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
+*          On entry, the K-by-N or M-by-K matrix A.
+*          On exit, A is overwritten by the corresponding block of 
+*          H*C or H^H*C or C*H or C*H^H.  See Futher Details.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A. 
+*          If SIDE = 'L', LDC >= max(1,K);
+*          If SIDE = 'R', LDC >= max(1,M). 
+*
+*  B       (input/output) REAL array, dimension (LDB,N)
+*          On entry, the M-by-N matrix B.
+*          On exit, B is overwritten by the corresponding block of
+*          H*C or H^H*C or C*H or C*H^H.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B. 
+*          LDB >= max(1,M).
+*
+*  WORK    (workspace) REAL array, dimension 
+*          (LDWORK,N) if SIDE = 'L',
+*          (LDWORK,K) if SIDE = 'R'.
+*
+*  LDWORK  (input) INTEGER
+*          The leading dimension of the array WORK.
+*          If SIDE = 'L', LDWORK >= K; 
+*          if SIDE = 'R', LDWORK >= M.
+*
+*  Further Details
+*  ===============
+*
+*  The matrix C is a composite matrix formed from blocks A and B.
+*  The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, 
+*  and if SIDE = 'L', A is of size K-by-N.
+*
+*  If SIDE = 'R' and DIRECT = 'F', C = [A B].
+*
+*  If SIDE = 'L' and DIRECT = 'F', C = [A] 
+*                                      [B].
+*
+*  If SIDE = 'R' and DIRECT = 'B', C = [B A].
+*
+*  If SIDE = 'L' and DIRECT = 'B', C = [B]
+*                                      [A]. 
+*
+*  The pentagonal matrix V is composed of a rectangular block V1 and a 
+*  trapezoidal block V2.  The size of the trapezoidal block is determined by 
+*  the parameter L, where 0<=L<=K.  If L=K, the V2 block of V is triangular;
+*  if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
+*
+*  If DIRECT = 'F' and STOREV = 'C':  V = [V1]
+*                                         [V2]
+*     - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
+*
+*  If DIRECT = 'F' and STOREV = 'R':  V = [V1 V2]
+*
+*     - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
+*
+*  If DIRECT = 'B' and STOREV = 'C':  V = [V2]
+*                                         [V1]
+*     - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
+*
+*  If DIRECT = 'B' and STOREV = 'R':  V = [V2 V1]
+*    
+*     - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
+*
+*  If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
+*
+*  If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
+*
+*  If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
+*
+*  If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
+*
+*  ==========================================================================
+*
+*     .. Parameters ..
+      REAL   ONE, ZERO
+      PARAMETER ( ONE = 1.0, ZERO = 0.0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER   I, J, MP, NP, KP
+      LOGICAL   LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
+*     ..
+*     .. External Functions ..
+      LOGICAL   LSAME
+      EXTERNAL  LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL  SCOPY, SGEMM, SLACGV, STRMM
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
+*
+      IF( LSAME( STOREV, 'C' ) ) THEN
+         COLUMN = .TRUE.
+         ROW = .FALSE.
+      ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
+         COLUMN = .FALSE.
+         ROW = .TRUE.
+      ELSE
+         COLUMN = .FALSE.
+         ROW = .FALSE.
+      END IF
+*
+      IF( LSAME( SIDE, 'L' ) ) THEN
+         LEFT = .TRUE.
+         RIGHT = .FALSE.
+      ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+         LEFT = .FALSE.
+         RIGHT = .TRUE.
+      ELSE
+         LEFT = .FALSE.
+         RIGHT = .FALSE.
+      END IF
+*
+      IF( LSAME( DIRECT, 'F' ) ) THEN
+         FORWARD = .TRUE.
+         BACKWARD = .FALSE.
+      ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
+         FORWARD = .FALSE.
+         BACKWARD = .TRUE.
+      ELSE
+         FORWARD = .FALSE.
+         BACKWARD = .FALSE.
+      END IF
+*
+* ---------------------------------------------------------------------------
+*      
+      IF( COLUMN .AND. FORWARD .AND. LEFT  ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ I ]    (K-by-K)
+*                  [ V ]    (M-by-K)
+*
+*        Form  H C  or  H^H C  where  C = [ A ]  (K-by-N)
+*                                         [ B ]  (M-by-N)
+*
+*        H = I - W T W^H          or  H^H = I - W T^H W^H
+*
+*        A = A -   T (A + V^H B)  or  A = A -   T^H (A + V^H B)
+*        B = B - V T (A + V^H B)  or  B = B - V T^H (A + V^H B) 
+*
+* ---------------------------------------------------------------------------
+*
+         MP = MIN( M-L+1, M )
+         KP = MIN( L+1, K )
+*         
+         DO J = 1, N
+            DO I = 1, L
+               WORK( I, J ) = B( M-L+I, J )
+            END DO
+         END DO
+         CALL STRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV,
+     $               WORK, LDWORK )         
+         CALL SGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL SGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, 
+     $               B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
+*     
+         DO J = 1, N
+            DO I = 1, K
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, 
+     $               WORK, LDWORK )
+*     
+         DO J = 1, N
+            DO I = 1, K
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL SGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ),  LDB )    
+         CALL STRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
+     $               WORK, LDWORK )
+         DO J = 1, N
+            DO I = 1, L
+               B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ I ]    (K-by-K)
+*                  [ V ]    (N-by-K)
+*
+*        Form  C H or  C H^H  where  C = [ A B ] (A is M-by-K, B is M-by-N)
+*
+*        H = I - W T W^H          or  H^H = I - W T^H W^H
+*
+*        A = A - (A + B V) T      or  A = A - (A + B V) T^H
+*        B = B - (A + B V) T V^H  or  B = B - (A + B V) T^H V^H
+*
+* ---------------------------------------------------------------------------
+*
+         NP = MIN( N-L+1, N )
+         KP = MIN( L+1, K )
+*         
+         DO J = 1, L
+            DO I = 1, M
+               WORK( I, J ) = B( I, N-L+J )
+            END DO
+         END DO
+         CALL STRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
+     $               WORK, LDWORK )
+         CALL SGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, 
+     $               V, LDV, ONE, WORK, LDWORK )
+         CALL SGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, 
+     $               V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
+*     
+         DO J = 1, K
+            DO I = 1, M
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, 
+     $               WORK, LDWORK )
+*     
+         DO J = 1, K
+            DO I = 1, M
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL SGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL STRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV,
+     $               WORK, LDWORK )
+         DO J = 1, L
+            DO I = 1, M
+               B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ V ]    (M-by-K)
+*                  [ I ]    (K-by-K)
+*
+*        Form  H C  or  H^H C  where  C = [ B ]  (M-by-N)
+*                                         [ A ]  (K-by-N)
+*
+*        H = I - W T W^H          or  H^H = I - W T^H W^H
+*
+*        A = A -   T (A + V^H B)  or  A = A -   T^H (A + V^H B)
+*        B = B - V T (A + V^H B)  or  B = B - V T^H (A + V^H B) 
+*
+* ---------------------------------------------------------------------------
+*
+         MP = MIN( L+1, M )
+         KP = MIN( K-L+1, K )
+*
+         DO J = 1, N
+            DO I = 1, L
+               WORK( K-L+I, J ) = B( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL SGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, 
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL SGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV,
+     $               B, LDB, ZERO, WORK, LDWORK )         
+*
+         DO J = 1, N
+            DO I = 1, K
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, 
+     $               WORK, LDWORK )
+*     
+         DO J = 1, N
+            DO I = 1, K
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, 
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL SGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
+     $               WORK, LDWORK, ONE, B,  LDB )
+         CALL STRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         DO J = 1, N
+            DO I = 1, L
+               B( I, J ) = B( I, J ) - WORK( K-L+I, J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ V ]    (N-by-K)
+*                  [ I ]    (K-by-K)
+*
+*        Form  C H  or  C H^H  where  C = [ B A ] (B is M-by-N, A is M-by-K)
+*
+*        H = I - W T W^H          or  H^H = I - W T^H W^H
+*
+*        A = A - (A + B V) T      or  A = A - (A + B V) T^H
+*        B = B - (A + B V) T V^H  or  B = B - (A + B V) T^H V^H
+*
+* ---------------------------------------------------------------------------
+*
+         NP = MIN( L+1, N )
+         KP = MIN( K-L+1, K )
+*         
+         DO J = 1, L
+            DO I = 1, M
+               WORK( I, K-L+J ) = B( I, J )
+            END DO
+         END DO
+         CALL STRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL SGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, 
+     $               V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL SGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, 
+     $               V, LDV, ZERO, WORK, LDWORK )
+*     
+         DO J = 1, K
+            DO I = 1, M
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, 
+     $               WORK, LDWORK )
+*     
+         DO J = 1, K
+            DO I = 1, M
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
+         CALL SGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL STRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         DO J = 1, L
+            DO I = 1, M
+               B( I, J ) = B( I, J ) - WORK( I, K-L+J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ I V ] ( I is K-by-K, V is K-by-M )
+*
+*        Form  H C  or  H^H C  where  C = [ A ]  (K-by-N)
+*                                         [ B ]  (M-by-N)
+*
+*        H = I - W^H T W          or  H^H = I - W^H T^H W
+*
+*        A = A -     T (A + V B)  or  A = A -     T^H (A + V B)
+*        B = B - V^H T (A + V B)  or  B = B - V^H T^H (A + V B) 
+*
+* ---------------------------------------------------------------------------
+*
+         MP = MIN( M-L+1, M )
+         KP = MIN( L+1, K )
+*
+         DO J = 1, N
+            DO I = 1, L
+               WORK( I, J ) = B( M-L+I, J )
+            END DO
+         END DO 
+         CALL STRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
+     $               WORK, LDB )
+         CALL SGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL SGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, 
+     $               B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
+*
+         DO J = 1, N
+            DO I = 1, K
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, 
+     $               WORK, LDWORK )
+*
+         DO J = 1, N
+            DO I = 1, K
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL SGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, 
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL STRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV,
+     $               WORK, LDWORK )
+         DO J = 1, N
+            DO I = 1, L
+               B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ I V ] ( I is K-by-K, V is K-by-N )
+*
+*        Form  C H  or  C H^H  where  C = [ A B ] (A is M-by-K, B is M-by-N)
+*
+*        H = I - W^H T W            or  H^H = I - W^H T^H W
+*
+*        A = A - (A + B V^H) T      or  A = A - (A + B V^H) T^H
+*        B = B - (A + B V^H) T V    or  B = B - (A + B V^H) T^H V
+*
+* ---------------------------------------------------------------------------
+*
+         NP = MIN( N-L+1, N )
+         KP = MIN( L+1, K )
+*
+         DO J = 1, L
+            DO I = 1, M
+               WORK( I, J ) = B( I, N-L+J )
+            END DO
+         END DO
+         CALL STRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV,
+     $               WORK, LDWORK )
+         CALL SGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV,
+     $               ONE, WORK, LDWORK )
+         CALL SGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, 
+     $               V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
+*
+         DO J = 1, K
+            DO I = 1, M
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, 
+     $               WORK, LDWORK )
+*
+         DO J = 1, K
+            DO I = 1, M
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V, LDV, ONE, B, LDB )
+         CALL SGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )   
+         CALL STRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
+     $               WORK, LDWORK )
+         DO J = 1, L
+            DO I = 1, M
+               B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ V I ] ( I is K-by-K, V is K-by-M )
+*
+*        Form  H C  or  H^H C  where  C = [ B ]  (M-by-N)
+*                                         [ A ]  (K-by-N)
+*
+*        H = I - W^H T W          or  H^H = I - W^H T^H W
+*
+*        A = A -     T (A + V B)  or  A = A -     T^H (A + V B)
+*        B = B - V^H T (A + V B)  or  B = B - V^H T^H (A + V B) 
+*
+* ---------------------------------------------------------------------------
+*
+         MP = MIN( L+1, M )
+         KP = MIN( K-L+1, K )
+*
+         DO J = 1, N
+            DO I = 1, L
+               WORK( K-L+I, J ) = B( I, J )
+            END DO
+         END DO
+         CALL STRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL SGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL SGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
+     $               ZERO, WORK, LDWORK )
+*
+         DO J = 1, N
+            DO I = 1, K
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
+     $               WORK, LDWORK )
+*
+         DO J = 1, N
+            DO I = 1, K
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL SGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, 
+     $               WORK, LDWORK, ONE, B, LDB )
+         CALL STRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV,
+     $               WORK( KP, 1 ), LDWORK )     
+         DO J = 1, N
+            DO I = 1, L
+               B( I, J ) = B( I, J ) - WORK( K-L+I, J )
+            END DO
+         END DO
+*
+* ---------------------------------------------------------------------------
+*      
+      ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
+*
+* ---------------------------------------------------------------------------
+*
+*        Let  W =  [ V I ] ( I is K-by-K, V is K-by-N )
+*
+*        Form  C H  or  C H^H  where  C = [ B A ] (A is M-by-K, B is M-by-N)
+*
+*        H = I - W^H T W            or  H^H = I - W^H T^H W
+*
+*        A = A - (A + B V^H) T      or  A = A - (A + B V^H) T^H
+*        B = B - (A + B V^H) T V    or  B = B - (A + B V^H) T^H V
+*
+* ---------------------------------------------------------------------------
+*
+         NP = MIN( L+1, N )
+         KP = MIN( K-L+1, K )
+*
+         DO J = 1, L
+            DO I = 1, M
+               WORK( I, K-L+J ) = B( I, J )
+            END DO
+         END DO
+         CALL STRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL SGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB,
+     $               V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL SGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV,
+     $               ZERO, WORK, LDWORK )                     
+*
+         DO J = 1, K
+            DO I = 1, M
+               WORK( I, J ) = WORK( I, J ) + A( I, J )
+            END DO
+         END DO
+*
+         CALL STRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,         
+     $               WORK, LDWORK )
+*
+         DO J = 1, K
+            DO I = 1, M
+               A( I, J ) = A( I, J ) - WORK( I, J )
+            END DO
+         END DO
+*
+         CALL SGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL SGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,    
+     $               V, LDV, ONE, B, LDB )
+         CALL STRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         DO J = 1, L
+            DO I = 1, M
+               B( I, J ) = B( I, J ) - WORK( I, K-L+J )
+            END DO
+         END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of STPRFB
+*
+      END