QR factorization for triangular-pentagonal matrices

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

16 files changed, 5518 insertions, 0 deletions

diff --git a/SRC/ctpmqrt.f b/SRC/ctpmqrt.f
new file mode 100644
index 00000000..1bcd7b7f
--- /dev/null
+++ b/SRC/ctpmqrt.f
@@ -0,0 +1,264 @@
+      SUBROUTINE CTPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT,
+     $                    A, LDA, B, LDB, WORK, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 SIDE, TRANS
+      INTEGER   INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT
+*     ..
+*     .. Array Arguments ..
+      COMPLEX   V( LDV, * ), A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CTPMQRT applies a complex orthogonal matrix Q obtained from a 
+*  "triangular-pentagonal" complex block reflector H to a general
+*  complex matrix C, which consists of two blocks A and B.
+*
+*  Arguments
+*  =========
+*
+*  SIDE    (input) CHARACTER*1
+*          = 'L': apply Q or Q**H from the Left;
+*          = 'R': apply Q or Q**H from the Right.
+*
+*  TRANS   (input) CHARACTER*1
+*          = 'N':  No transpose, apply Q;
+*          = 'C':  Transpose, apply Q**H.
+*
+*  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 number of elementary reflectors whose product defines
+*          the matrix Q.
+*
+*  L       (input) INTEGER
+*          The order of the trapezoidal part of V.  
+*          K >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size used for the storage of T.  K >= NB >= 1.
+*          This must be the same value of NB used to generate T
+*          in CTPQRT.
+*
+*  V       (input) COMPLEX array, dimension (LDA,K)
+*          The i-th column must contain the vector which defines the
+*          elementary reflector H(i), for i = 1,2,...,k, as returned by
+*          CTPQRT in B.  See Further Details.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.
+*          If SIDE = 'L', LDA >= max(1,M);
+*          if SIDE = 'R', LDA >= max(1,N).
+*
+*  T       (input) COMPLEX array, dimension (LDT,K)
+*          The upper triangular factors of the block reflectors
+*          as returned by CTPQRT, stored as a NB-by-K matrix.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  A       (input/output) COMPLEX 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 
+*          Q*C or Q**H*C or C*Q or C*Q**H.  See Further 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) COMPLEX array, dimension (LDB,N)
+*          On entry, the M-by-N matrix B.
+*          On exit, B is overwritten by the corresponding block of
+*          Q*C or Q**H*C or C*Q or C*Q**H.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B. 
+*          LDB >= max(1,M).
+*
+*  WORK    (workspace/output) COMPLEX array.  The dimension of WORK is
+*           N*NB if SIDE = 'L', or  M*NB if SIDE = 'R'.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The columns of the pentagonal matrix V contain the elementary reflectors
+*  H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a 
+*  trapezoidal block V2:
+*
+*        V = [V1]
+*            [V2].
+*
+*  The size of the trapezoidal block V2 is determined by the parameter L, 
+*  where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L
+*  rows of a K-by-K upper triangular matrix.  If L=K, V2 is upper triangular;
+*  if L=0, there is no trapezoidal block, hence V = V1 is rectangular.
+*
+*  If SIDE = 'L':  C = [A]  where A is K-by-N,  B is M-by-N and V is M-by-K. 
+*                      [B]   
+*  
+*  If SIDE = 'R':  C = [A B]  where A is M-by-K, B is M-by-N and V is N-by-K.
+*
+*  The complex orthogonal matrix Q is formed from V and T.
+*
+*  If TRANS='N' and SIDE='L', C is on exit replaced with Q * C.
+*
+*  If TRANS='C' and SIDE='L', C is on exit replaced with Q**H * C.
+*
+*  If TRANS='N' and SIDE='R', C is on exit replaced with C * Q.
+*
+*  If TRANS='C' and SIDE='R', C is on exit replaced with C * Q**H.
+*
+*  =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
+      INTEGER            I, IB, MB, LB, KF, Q
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, CLARFB
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     .. Test the input arguments ..
+*
+      INFO   = 0
+      LEFT   = LSAME( SIDE,  'L' )
+      RIGHT  = LSAME( SIDE,  'R' )
+      TRAN   = LSAME( TRANS, 'C' )
+      NOTRAN = LSAME( TRANS, 'N' )
+*      
+      IF( LEFT ) THEN
+         Q = M
+      ELSE IF ( RIGHT ) THEN
+         Q = N
+      END IF
+      IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
+         INFO = -1
+      ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
+         INFO = -2
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( K.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
+         INFO = -6         
+      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+         INFO = -7
+      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CTPMQRT', -INFO )
+         RETURN
+      END IF
+*
+*     .. Quick return if possible ..
+*
+      IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN
+*
+      IF( LEFT .AND. TRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF
+            CALL CTPRFB( 'L', 'C', 'F', 'C', MB, N, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*         
+      ELSE IF( RIGHT .AND. NOTRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL CTPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      ELSE IF( LEFT .AND. NOTRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )  
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF                   
+            CALL CTPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*
+      ELSE IF( RIGHT .AND. TRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )         
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL CTPRFB( 'R', 'C', 'F', 'C', M, MB, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of CTPMQRT
+*
+      END
diff --git a/SRC/ctpqrt.f b/SRC/ctpqrt.f
new file mode 100644
index 00000000..ada345b4
--- /dev/null
+++ b/SRC/ctpqrt.f
@@ -0,0 +1,185 @@
+      SUBROUTINE CTPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK,
+     $                   INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 ..
+      INTEGER INFO, LDA, LDB, LDT, N, M, L, NB
+*     ..
+*     .. Array Arguments ..
+      COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CTPQRT computes a blocked QR factorization of a complex 
+*  "triangular-pentagonal" matrix C, which is composed of a 
+*  triangular block A and pentagonal block B, using the compact 
+*  WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of the
+*          triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size to be used in the blocked QR.  N >= NB >= 1.
+*
+*  A       (input/output) COMPLEX array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) COMPLEX array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) COMPLEX array, dimension (LDT,N)
+*          The upper triangular block reflectors stored in compact form
+*          as a sequence of upper triangular blocks.  See Further Details.
+*          
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  WORK    (workspace) COMPLEX array, dimension (NB*N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*  
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*
+*  The number of blocks is B = ceiling(N/NB), where each
+*  block is of order NB except for the last block, which is of order 
+*  IB = N - (B-1)*NB.  For each of the B blocks, a upper triangular block
+*  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB 
+*  for the last block) T's are stored in the NB-by-N matrix T as
+*
+*               T = [T1 T2 ... TB].
+*
+* =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      INTEGER    I, IB, LB, MB, IINFO
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL   CTPQRT2, CTPRFB, XERBLA
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+         INFO = -4
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CTPQRT', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( M.EQ.0 .OR. N.EQ.0 ) RETURN
+*
+      DO I = 1, N, NB
+*     
+*     Compute the QR factorization of the current block
+*
+         IB = MIN( N-I+1, NB )
+         MB = MIN( M-L+I+IB-1, M )
+         IF( I.GE.L ) THEN
+            LB = 0
+         ELSE
+            LB = MB-M+L-I+1
+         END IF
+*
+         CALL CTPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB, 
+     $                 T(1, I ), LDT, IINFO )
+*
+*     Update by applying H**H to B(:,I+IB:N) from the left
+*
+         IF( I+IB.LE.N ) THEN
+            CALL CTPRFB( 'L', 'C', 'F', 'C', MB, N-I-IB+1, IB, LB,
+     $                    B( 1, I ), LDB, T( 1, I ), LDT, 
+     $                    A( I, I+IB ), LDA, B( 1, I+IB ), LDB, 
+     $                    WORK, IB )
+         END IF
+      END DO
+      RETURN
+*     
+*     End of CTPQRT
+*
+      END
diff --git a/SRC/ctpqrt2.f b/SRC/ctpqrt2.f
new file mode 100644
index 00000000..33ea7b92
--- /dev/null
+++ b/SRC/ctpqrt2.f
@@ -0,0 +1,224 @@
+      SUBROUTINE CTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK 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 ..
+      INTEGER   INFO, LDA, LDB, LDT, N, M, L
+*     ..
+*     .. Array Arguments ..
+      COMPLEX   A( LDA, * ), B( LDB, * ), T( LDT, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CTPQRT2 computes a QR factorization of a complex "triangular-pentagonal"
+*  matrix C, which is composed of a triangular block A and pentagonal block B, 
+*  using the compact WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The total number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of
+*          the triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.  
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  A       (input/output) COMPLEX array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) COMPLEX array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) COMPLEX array, dimension (LDT,N)
+*          The N-by-N upper triangular factor T of the block reflector.
+*          See Further Details.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= max(1,N)
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*  The (M+N)-by-(M+N) block reflector H is then given by
+*
+*               H = I - W * T * W**H
+*
+*  where W**H is the conjugate transpose of W and T is the upper triangular
+*  factor of the block reflector.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX  ONE, ZERO
+      PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER   I, J, P, MP, NP
+      COMPLEX   ALPHA
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL  CLARFG, CGEMV, CGERC, CTRMV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -7
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CTPQRT2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
+*      
+      DO I = 1, N
+*
+*        Generate elementary reflector H(I) to annihilate B(:,I)
+*
+         P = M-L+MIN( L, I )
+         CALL CLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
+         IF( I.LT.N ) THEN
+*
+*           W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)]
+*
+            DO J = 1, N-I
+               T( J, N ) = CONJG(A( I, I+J ))
+            END DO
+            CALL CGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB, 
+     $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
+*
+*           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H
+*
+            ALPHA = -CONJG(T( I, 1 ))            
+            DO J = 1, N-I
+               A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N ))
+            END DO
+            CALL CGERC( P, N-I, ALPHA, B( 1, I ), 1, 
+     $           T( 1, N ), 1, B( 1, I+1 ), LDB )
+         END IF
+      END DO
+*
+      DO I = 2, N
+*
+*        T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I))
+*
+         ALPHA = -T( I, 1 )
+
+         DO J = 1, I-1
+            T( J, I ) = ZERO
+         END DO
+         P = MIN( I-1, L )
+         MP = MIN( M-L+1, M )
+         NP = MIN( P+1, N )
+*
+*        Triangular part of B2
+*
+         DO J = 1, P
+            T( J, I ) = ALPHA*B( M-L+J, I )
+         END DO
+         CALL CTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB,
+     $               T( 1, I ), 1 )
+*
+*        Rectangular part of B2
+*
+         CALL CGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB, 
+     $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
+*
+*        B1
+*
+         CALL CGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, 
+     $               ONE, T( 1, I ), 1 )         
+*
+*        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
+*
+         CALL CTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
+*
+*        T(I,I) = tau(I)
+*
+         T( I, I ) = T( I, 1 )
+         T( I, 1 ) = ZERO
+      END DO
+   
+*
+*     End of CTPQRT2
+*
+      END
diff --git a/SRC/ctprfb.f b/SRC/ctprfb.f
new file mode 100644
index 00000000..9c7980c0
--- /dev/null
+++ b/SRC/ctprfb.f
@@ -0,0 +1,708 @@
+      SUBROUTINE CTPRFB( 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 ..
+      COMPLEX   A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          V( LDV, * ), WORK( LDWORK, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CTPRFB applies a complex "triangular-pentagonal" block reflector H or its 
+*  conjugate transpose H**H to a complex 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) COMPLEX 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) COMPLEX 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) COMPLEX 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) COMPLEX 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) COMPLEX 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 ..
+      COMPLEX   ONE, ZERO
+      PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,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  CCOPY, CGEMM, CLACGV, CTRMM
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC CONJG
+*     ..
+*     .. 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 CTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV,
+     $               WORK, LDWORK )         
+         CALL CGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL CGEMM( 'C', '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 CTRMM( '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 CGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL CGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ),  LDB )    
+         CALL CTRMM( '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 CTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
+     $               WORK, LDWORK )
+         CALL CGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, 
+     $               V, LDV, ONE, WORK, LDWORK )
+         CALL CGEMM( '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 CTRMM( '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 CGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL CGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL CTRMM( 'R', 'U', 'C', '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 CTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL CGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, 
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL CGEMM( 'C', '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 CTRMM( '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 CGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, 
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL CGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
+     $               WORK, LDWORK, ONE, B,  LDB )
+         CALL CTRMM( '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 CTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL CGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, 
+     $               V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL CGEMM( '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 CTRMM( '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 CGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
+         CALL CGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL CTRMM( 'R', 'L', 'C', '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 CTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
+     $               WORK, LDB )
+         CALL CGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL CGEMM( '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 CTRMM( '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 CGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL CGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, 
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL CTRMM( 'L', 'L', 'C', '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 CTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV,
+     $               WORK, LDWORK )
+         CALL CGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV,
+     $               ONE, WORK, LDWORK )
+         CALL CGEMM( 'N', 'C', 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 CTRMM( '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 CGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V, LDV, ONE, B, LDB )
+         CALL CGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )   
+         CALL CTRMM( '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 CTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL CGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL CGEMM( '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 CTRMM( '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 CGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL CGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV, 
+     $               WORK, LDWORK, ONE, B, LDB )
+         CALL CTRMM( 'L', 'U', 'C', '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 CTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL CGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB,
+     $               V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL CGEMM( 'N', 'C', 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 CTRMM( '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 CGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL CGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,    
+     $               V, LDV, ONE, B, LDB )
+         CALL CTRMM( '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 CTPRFB
+*
+      END
diff --git a/SRC/dtpmqrt.f b/SRC/dtpmqrt.f
new file mode 100644
index 00000000..73b257a4
--- /dev/null
+++ b/SRC/dtpmqrt.f
@@ -0,0 +1,264 @@
+      SUBROUTINE DTPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT,
+     $                    A, LDA, B, LDB, WORK, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 SIDE, TRANS
+      INTEGER   INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   V( LDV, * ), A( LDA, * ), B( LDB, * ), 
+     $                   T( LDT, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DTPMQRT applies a real orthogonal matrix Q obtained from a 
+*  "triangular-pentagonal" real block reflector H to a general
+*  real matrix C, which consists of two blocks A and B.
+*
+*  Arguments
+*  =========
+*
+*  SIDE    (input) CHARACTER*1
+*          = 'L': apply Q or Q**T from the Left;
+*          = 'R': apply Q or Q**T from the Right.
+*
+*  TRANS   (input) CHARACTER*1
+*          = 'N':  No transpose, apply Q;
+*          = 'C':  Transpose, apply Q**T.
+*
+*  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 number of elementary reflectors whose product defines
+*          the matrix Q.
+*
+*  L       (input) INTEGER
+*          The order of the trapezoidal part of V.  
+*          K >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size used for the storage of T.  K >= NB >= 1.
+*          This must be the same value of NB used to generate T
+*          in CTPQRT.
+*
+*  V       (input) DOUBLE PRECISION array, dimension (LDA,K)
+*          The i-th column must contain the vector which defines the
+*          elementary reflector H(i), for i = 1,2,...,k, as returned by
+*          CTPQRT in B.  See Further Details.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.
+*          If SIDE = 'L', LDA >= max(1,M);
+*          if SIDE = 'R', LDA >= max(1,N).
+*
+*  T       (input) DOUBLE PRECISION array, dimension (LDT,K)
+*          The upper triangular factors of the block reflectors
+*          as returned by CTPQRT, stored as a NB-by-K matrix.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  A       (input/output) DOUBLE PRECISION 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 
+*          Q*C or Q**T*C or C*Q or C*Q**T.  See Further 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) DOUBLE PRECISION array, dimension (LDB,N)
+*          On entry, the M-by-N matrix B.
+*          On exit, B is overwritten by the corresponding block of
+*          Q*C or Q**T*C or C*Q or C*Q**T.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B. 
+*          LDB >= max(1,M).
+*
+*  WORK    (workspace/output) DOUBLE PRECISION array.  The dimension of WORK is
+*           N*NB if SIDE = 'L', or  M*NB if SIDE = 'R'.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The columns of the pentagonal matrix V contain the elementary reflectors
+*  H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a 
+*  trapezoidal block V2:
+*
+*        V = [V1]
+*            [V2].
+*
+*  The size of the trapezoidal block V2 is determined by the parameter L, 
+*  where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L
+*  rows of a K-by-K upper triangular matrix.  If L=K, V2 is upper triangular;
+*  if L=0, there is no trapezoidal block, hence V = V1 is rectangular.
+*
+*  If SIDE = 'L':  C = [A]  where A is K-by-N,  B is M-by-N and V is M-by-K. 
+*                      [B]   
+*  
+*  If SIDE = 'R':  C = [A B]  where A is M-by-K, B is M-by-N and V is N-by-K.
+*
+*  The real orthogonal matrix Q is formed from V and T.
+*
+*  If TRANS='N' and SIDE='L', C is on exit replaced with Q * C.
+*
+*  If TRANS='T' and SIDE='L', C is on exit replaced with Q**T * C.
+*
+*  If TRANS='N' and SIDE='R', C is on exit replaced with C * Q.
+*
+*  If TRANS='T' and SIDE='R', C is on exit replaced with C * Q**T.
+*
+*  =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
+      INTEGER            I, IB, MB, LB, KF, Q
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, DLARFB
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     .. Test the input arguments ..
+*
+      INFO   = 0
+      LEFT   = LSAME( SIDE,  'L' )
+      RIGHT  = LSAME( SIDE,  'R' )
+      TRAN   = LSAME( TRANS, 'T' )
+      NOTRAN = LSAME( TRANS, 'N' )
+*      
+      IF( LEFT ) THEN
+         Q = M
+      ELSE IF ( RIGHT ) THEN
+         Q = N
+      END IF
+      IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
+         INFO = -1
+      ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
+         INFO = -2
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( K.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
+         INFO = -6         
+      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+         INFO = -7
+      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DTPMQRT', -INFO )
+         RETURN
+      END IF
+*
+*     .. Quick return if possible ..
+*
+      IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN
+*
+      IF( LEFT .AND. TRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF
+            CALL DTPRFB( 'L', 'T', 'F', 'C', MB, N, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*         
+      ELSE IF( RIGHT .AND. NOTRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL DTPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      ELSE IF( LEFT .AND. NOTRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )  
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF                   
+            CALL DTPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*
+      ELSE IF( RIGHT .AND. TRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )         
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL DTPRFB( 'R', 'T', 'F', 'C', M, MB, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of DTPMQRT
+*
+      END
diff --git a/SRC/dtpqrt.f b/SRC/dtpqrt.f
new file mode 100644
index 00000000..bf1aabb6
--- /dev/null
+++ b/SRC/dtpqrt.f
@@ -0,0 +1,185 @@
+      SUBROUTINE DTPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK,
+     $                   INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 ..
+      INTEGER INFO, LDA, LDB, LDT, N, M, L, NB
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DTPQRT computes a blocked QR factorization of a real 
+*  "triangular-pentagonal" matrix C, which is composed of a 
+*  triangular block A and pentagonal block B, using the compact 
+*  WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of the
+*          triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size to be used in the blocked QR.  N >= NB >= 1.
+*
+*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) DOUBLE PRECISION array, dimension (LDT,N)
+*          The upper triangular block reflectors stored in compact form
+*          as a sequence of upper triangular blocks.  See Further Details.
+*          
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  WORK    (workspace) DOUBLE PRECISION array, dimension (NB*N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*  
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*
+*  The number of blocks is B = ceiling(N/NB), where each
+*  block is of order NB except for the last block, which is of order 
+*  IB = N - (B-1)*NB.  For each of the B blocks, a upper triangular block
+*  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB 
+*  for the last block) T's are stored in the NB-by-N matrix T as
+*
+*               T = [T1 T2 ... TB].
+*
+* =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      INTEGER    I, IB, LB, MB, IINFO
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL   DTPQRT2, DTPRFB, XERBLA
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+         INFO = -4
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DTPQRT', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( M.EQ.0 .OR. N.EQ.0 ) RETURN
+*
+      DO I = 1, N, NB
+*     
+*     Compute the QR factorization of the current block
+*
+         IB = MIN( N-I+1, NB )
+         MB = MIN( M-L+I+IB-1, M )
+         IF( I.GE.L ) THEN
+            LB = 0
+         ELSE
+            LB = MB-M+L-I+1
+         END IF
+*
+         CALL DTPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB, 
+     $                 T(1, I ), LDT, IINFO )
+*
+*     Update by applying H**T to B(:,I+IB:N) from the left
+*
+         IF( I+IB.LE.N ) THEN
+            CALL DTPRFB( 'L', 'T', 'F', 'C', MB, N-I-IB+1, IB, LB,
+     $                    B( 1, I ), LDB, T( 1, I ), LDT, 
+     $                    A( I, I+IB ), LDA, B( 1, I+IB ), LDB, 
+     $                    WORK, IB )
+         END IF
+      END DO
+      RETURN
+*     
+*     End of DTPQRT
+*
+      END
diff --git a/SRC/dtpqrt2.f b/SRC/dtpqrt2.f
new file mode 100644
index 00000000..a6fd36b3
--- /dev/null
+++ b/SRC/dtpqrt2.f
@@ -0,0 +1,224 @@
+      SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK 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 ..
+      INTEGER   INFO, LDA, LDB, LDT, N, M, L
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), T( LDT, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DTPQRT2 computes a QR factorization of a real "triangular-pentagonal"
+*  matrix C, which is composed of a triangular block A and pentagonal block B, 
+*  using the compact WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The total number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of
+*          the triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.  
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) DOUBLE PRECISION array, dimension (LDT,N)
+*          The N-by-N upper triangular factor T of the block reflector.
+*          See Further Details.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= max(1,N)
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*  The (M+N)-by-(M+N) block reflector H is then given by
+*
+*               H = I - W * T * W**T
+*
+*  where W^H is the conjugate transpose of W and T is the upper triangular
+*  factor of the block reflector.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION  ONE, ZERO
+      PARAMETER( ONE = 1.0, ZERO = 0.0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER   I, J, P, MP, NP
+      DOUBLE PRECISION   ALPHA
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL  DLARFG, DGEMV, DGER, DTRMV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -7
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DTPQRT2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
+*      
+      DO I = 1, N
+*
+*        Generate elementary reflector H(I) to annihilate B(:,I)
+*
+         P = M-L+MIN( L, I )
+         CALL DLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
+         IF( I.LT.N ) THEN
+*
+*           W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)]
+*
+            DO J = 1, N-I
+               T( J, N ) = (A( I, I+J ))
+            END DO
+            CALL DGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB, 
+     $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
+*
+*           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H
+*
+            ALPHA = -(T( I, 1 ))            
+            DO J = 1, N-I
+               A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N ))
+            END DO
+            CALL DGER( P, N-I, ALPHA, B( 1, I ), 1, 
+     $           T( 1, N ), 1, B( 1, I+1 ), LDB )
+         END IF
+      END DO
+*
+      DO I = 2, N
+*
+*        T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I))
+*
+         ALPHA = -T( I, 1 )
+
+         DO J = 1, I-1
+            T( J, I ) = ZERO
+         END DO
+         P = MIN( I-1, L )
+         MP = MIN( M-L+1, M )
+         NP = MIN( P+1, N )
+*
+*        Triangular part of B2
+*
+         DO J = 1, P
+            T( J, I ) = ALPHA*B( M-L+J, I )
+         END DO
+         CALL DTRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB,
+     $               T( 1, I ), 1 )
+*
+*        Rectangular part of B2
+*
+         CALL DGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB, 
+     $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
+*
+*        B1
+*
+         CALL DGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, 
+     $               ONE, T( 1, I ), 1 )         
+*
+*        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
+*
+         CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
+*
+*        T(I,I) = tau(I)
+*
+         T( I, I ) = T( I, 1 )
+         T( I, 1 ) = ZERO
+      END DO
+   
+*
+*     End of DTPQRT2
+*
+      END
diff --git a/SRC/dtprfb.f b/SRC/dtprfb.f
new file mode 100644
index 00000000..c7b4d58a
--- /dev/null
+++ b/SRC/dtprfb.f
@@ -0,0 +1,705 @@
+      SUBROUTINE DTPRFB( 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 ..
+      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          V( LDV, * ), WORK( LDWORK, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DTPRFB applies a real "triangular-pentagonal" block reflector H or its 
+*  transpose H**T 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**T from the Left
+*          = 'R': apply H or H**T from the Right
+*
+*  TRANS   (input) CHARACTER*1
+*          = 'N': apply H (No transpose)
+*          = 'T': apply H**T (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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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**T*C or C*H or C*H**T.  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) DOUBLE PRECISION 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**T*C or C*H or C*H**T.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B. 
+*          LDB >= max(1,M).
+*
+*  WORK    (workspace) DOUBLE PRECISION 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 ..
+      DOUBLE PRECISION   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  DGEMM, DTRMM
+*     ..
+*     .. 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**T C  where  C = [ A ]  (K-by-N)
+*                                          [ B ]  (M-by-N)
+*
+*        H = I - W T W**T          or  H**T = I - W T**T W**T
+*
+*        A = A -   T (A + V**T B)  or  A = A -   T**T (A + V**T B)
+*        B = B - V T (A + V**T B)  or  B = B - V T**T (A + V**T 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 DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV,
+     $               WORK, LDWORK )         
+         CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ),  LDB )    
+         CALL DTRMM( '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**T  where  C = [ A B ] (A is M-by-K, B is M-by-N)
+*
+*        H = I - W T W**T          or  H**T = I - W T**T W**T
+*
+*        A = A - (A + B V) T      or  A = A - (A + B V) T**T
+*        B = B - (A + B V) T V**T  or  B = B - (A + B V) T**T V**T
+*
+* ---------------------------------------------------------------------------
+*
+         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 DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
+     $               WORK, LDWORK )
+         CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, 
+     $               V, LDV, ONE, WORK, LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL DTRMM( '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**T C  where  C = [ B ]  (M-by-N)
+*                                          [ A ]  (K-by-N)
+*
+*        H = I - W T W**T          or  H**T = I - W T**T W**T
+*
+*        A = A -   T (A + V**T B)  or  A = A -   T**T (A + V**T B)
+*        B = B - V T (A + V**T B)  or  B = B - V T**T (A + V**T 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 DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, 
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, 
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
+     $               WORK, LDWORK, ONE, B,  LDB )
+         CALL DTRMM( '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**T  where  C = [ B A ] (B is M-by-N, A is M-by-K)
+*
+*        H = I - W T W**T          or  H**T = I - W T**T W**T
+*
+*        A = A - (A + B V) T      or  A = A - (A + B V) T**T
+*        B = B - (A + B V) T V**T  or  B = B - (A + B V) T**T V**T
+*
+* ---------------------------------------------------------------------------
+*
+         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 DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, 
+     $               V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
+         CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL DTRMM( '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**T C  where  C = [ A ]  (K-by-N)
+*                                          [ B ]  (M-by-N)
+*
+*        H = I - W**T T W          or  H**T = I - W**T T**T W
+*
+*        A = A -     T (A + V B)  or  A = A -     T**T (A + V B)
+*        B = B - V**T T (A + V B)  or  B = B - V**T T**T (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 DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
+     $               WORK, LDB )
+         CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, 
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL DTRMM( '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**T  where  C = [ A B ] (A is M-by-K, B is M-by-N)
+*
+*        H = I - W**T T W            or  H**T = I - W**T T**T W
+*
+*        A = A - (A + B V**T) T      or  A = A - (A + B V**T) T**T
+*        B = B - (A + B V**T) T V    or  B = B - (A + B V**T) T**T 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 DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV,
+     $               WORK, LDWORK )
+         CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV,
+     $               ONE, WORK, LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V, LDV, ONE, B, LDB )
+         CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )   
+         CALL DTRMM( '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**T C  where  C = [ B ]  (M-by-N)
+*                                          [ A ]  (K-by-N)
+*
+*        H = I - W**T T W          or  H**T = I - W**T T**T W
+*
+*        A = A -     T (A + V B)  or  A = A -     T**T (A + V B)
+*        B = B - V**T T (A + V B)  or  B = B - V**T T**T (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 DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, 
+     $               WORK, LDWORK, ONE, B, LDB )
+         CALL DTRMM( '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**T  where  C = [ B A ] (A is M-by-K, B is M-by-N)
+*
+*        H = I - W**T T W            or  H**T = I - W**T T**T W
+*
+*        A = A - (A + B V**T) T      or  A = A - (A + B V**T) T**T
+*        B = B - (A + B V**T) T V    or  B = B - (A + B V**T) T**T 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 DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB,
+     $               V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL DGEMM( '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 DTRMM( '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 DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,    
+     $               V, LDV, ONE, B, LDB )
+         CALL DTRMM( '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 DTPRFB
+*
+      END
diff --git a/SRC/stpmqrt.f b/SRC/stpmqrt.f
new file mode 100644
index 00000000..c6c2c858
--- /dev/null
+++ b/SRC/stpmqrt.f
@@ -0,0 +1,264 @@
+      SUBROUTINE STPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT,
+     $                    A, LDA, B, LDB, WORK, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 SIDE, TRANS
+      INTEGER   INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT
+*     ..
+*     .. Array Arguments ..
+      REAL   V( LDV, * ), A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  STPMQRT applies a real orthogonal matrix Q obtained from a 
+*  "triangular-pentagonal" real block reflector H to a general
+*  real matrix C, which consists of two blocks A and B.
+*
+*  Arguments
+*  =========
+*
+*  SIDE    (input) CHARACTER*1
+*          = 'L': apply Q or Q^H from the Left;
+*          = 'R': apply Q or Q^H from the Right.
+*
+*  TRANS   (input) CHARACTER*1
+*          = 'N':  No transpose, apply Q;
+*          = 'C':  Transpose, apply Q^H.
+*
+*  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 number of elementary reflectors whose product defines
+*          the matrix Q.
+*
+*  L       (input) INTEGER
+*          The order of the trapezoidal part of V.  
+*          K >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size used for the storage of T.  K >= NB >= 1.
+*          This must be the same value of NB used to generate T
+*          in CTPQRT.
+*
+*  V       (input) REAL array, dimension (LDA,K)
+*          The i-th column must contain the vector which defines the
+*          elementary reflector H(i), for i = 1,2,...,k, as returned by
+*          CTPQRT in B.  See Further Details.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.
+*          If SIDE = 'L', LDA >= max(1,M);
+*          if SIDE = 'R', LDA >= max(1,N).
+*
+*  T       (input) REAL array, dimension (LDT,K)
+*          The upper triangular factors of the block reflectors
+*          as returned by CTPQRT, stored as a NB-by-K matrix.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  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 
+*          Q*C or Q^H*C or C*Q or C*Q^H.  See Further 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
+*          Q*C or Q^H*C or C*Q or C*Q^H.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B. 
+*          LDB >= max(1,M).
+*
+*  WORK    (workspace/output) REAL array.  The dimension of WORK is
+*           N*NB if SIDE = 'L', or  M*NB if SIDE = 'R'.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The columns of the pentagonal matrix V contain the elementary reflectors
+*  H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a 
+*  trapezoidal block V2:
+*
+*        V = [V1]
+*            [V2].
+*
+*  The size of the trapezoidal block V2 is determined by the parameter L, 
+*  where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L
+*  rows of a K-by-K upper triangular matrix.  If L=K, V2 is upper triangular;
+*  if L=0, there is no trapezoidal block, hence V = V1 is rectangular.
+*
+*  If SIDE = 'L':  C = [A]  where A is K-by-N,  B is M-by-N and V is M-by-K. 
+*                      [B]   
+*  
+*  If SIDE = 'R':  C = [A B]  where A is M-by-K, B is M-by-N and V is N-by-K.
+*
+*  The real orthogonal matrix Q is formed from V and T.
+*
+*  If TRANS='N' and SIDE='L', C is on exit replaced with Q * C.
+*
+*  If TRANS='C' and SIDE='L', C is on exit replaced with Q^H * C.
+*
+*  If TRANS='N' and SIDE='R', C is on exit replaced with C * Q.
+*
+*  If TRANS='C' and SIDE='R', C is on exit replaced with C * Q^H.
+*
+*  =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
+      INTEGER            I, IB, MB, LB, KF, Q
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, SLARFB
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     .. Test the input arguments ..
+*
+      INFO   = 0
+      LEFT   = LSAME( SIDE,  'L' )
+      RIGHT  = LSAME( SIDE,  'R' )
+      TRAN   = LSAME( TRANS, 'T' )
+      NOTRAN = LSAME( TRANS, 'N' )
+*      
+      IF( LEFT ) THEN
+         Q = M
+      ELSE IF ( RIGHT ) THEN
+         Q = N
+      END IF
+      IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
+         INFO = -1
+      ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
+         INFO = -2
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( K.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
+         INFO = -6         
+      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+         INFO = -7
+      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'STPMQRT', -INFO )
+         RETURN
+      END IF
+*
+*     .. Quick return if possible ..
+*
+      IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN
+*
+      IF( LEFT .AND. TRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF
+            CALL STPRFB( 'L', 'T', 'F', 'C', MB, N, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*         
+      ELSE IF( RIGHT .AND. NOTRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL STPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      ELSE IF( LEFT .AND. NOTRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )  
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF                   
+            CALL STPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*
+      ELSE IF( RIGHT .AND. TRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )         
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL STPRFB( 'R', 'T', 'F', 'C', M, MB, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of STPMQRT
+*
+      END
diff --git a/SRC/stpqrt.f b/SRC/stpqrt.f
new file mode 100644
index 00000000..07b9a9a8
--- /dev/null
+++ b/SRC/stpqrt.f
@@ -0,0 +1,185 @@
+      SUBROUTINE STPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK,
+     $                   INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 ..
+      INTEGER INFO, LDA, LDB, LDT, N, M, L, NB
+*     ..
+*     .. Array Arguments ..
+      REAL A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  STPQRT computes a blocked QR factorization of a real 
+*  "triangular-pentagonal" matrix C, which is composed of a 
+*  triangular block A and pentagonal block B, using the compact 
+*  WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of the
+*          triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size to be used in the blocked QR.  N >= NB >= 1.
+*
+*  A       (input/output) REAL array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) REAL array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) REAL array, dimension (LDT,N)
+*          The upper triangular block reflectors stored in compact form
+*          as a sequence of upper triangular blocks.  See Further Details.
+*          
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  WORK    (workspace) REAL array, dimension (NB*N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*  
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*
+*  The number of blocks is B = ceiling(N/NB), where each
+*  block is of order NB except for the last block, which is of order 
+*  IB = N - (B-1)*NB.  For each of the B blocks, a upper triangular block
+*  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB 
+*  for the last block) T's are stored in the NB-by-N matrix T as
+*
+*               T = [T1 T2 ... TB].
+*
+* =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      INTEGER    I, IB, LB, MB, IINFO
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL   STPQRT2, STPRFB, XERBLA
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+         INFO = -4
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'STPQRT', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( M.EQ.0 .OR. N.EQ.0 ) RETURN
+*
+      DO I = 1, N, NB
+*     
+*     Compute the QR factorization of the current block
+*
+         IB = MIN( N-I+1, NB )
+         MB = MIN( M-L+I+IB-1, M )
+         IF( I.GE.L ) THEN
+            LB = 0
+         ELSE
+            LB = MB-M+L-I+1
+         END IF
+*
+         CALL STPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB, 
+     $                 T(1, I ), LDT, IINFO )
+*
+*     Update by applying H^H to B(:,I+IB:N) from the left
+*
+         IF( I+IB.LE.N ) THEN
+            CALL STPRFB( 'L', 'T', 'F', 'C', MB, N-I-IB+1, IB, LB,
+     $                    B( 1, I ), LDB, T( 1, I ), LDT, 
+     $                    A( I, I+IB ), LDA, B( 1, I+IB ), LDB, 
+     $                    WORK, IB )
+         END IF
+      END DO
+      RETURN
+*     
+*     End of STPQRT
+*
+      END
diff --git a/SRC/stpqrt2.f b/SRC/stpqrt2.f
new file mode 100644
index 00000000..edbed209
--- /dev/null
+++ b/SRC/stpqrt2.f
@@ -0,0 +1,224 @@
+      SUBROUTINE STPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK 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 ..
+      INTEGER   INFO, LDA, LDB, LDT, N, M, L
+*     ..
+*     .. Array Arguments ..
+      REAL   A( LDA, * ), B( LDB, * ), T( LDT, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  STPQRT2 computes a QR factorization of a real "triangular-pentagonal"
+*  matrix C, which is composed of a triangular block A and pentagonal block B, 
+*  using the compact WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The total number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of
+*          the triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.  
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  A       (input/output) REAL array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) REAL array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) REAL array, dimension (LDT,N)
+*          The N-by-N upper triangular factor T of the block reflector.
+*          See Further Details.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= max(1,N)
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*  The (M+N)-by-(M+N) block reflector H is then given by
+*
+*               H = I - W * T * W^H
+*
+*  where W^H is the conjugate transpose of W and T is the upper triangular
+*  factor of the block reflector.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL  ONE, ZERO
+      PARAMETER( ONE = 1.0, ZERO = 0.0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER   I, J, P, MP, NP
+      REAL   ALPHA
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL  SLARFG, SGEMV, SGER, STRMV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -7
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'STPQRT2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
+*      
+      DO I = 1, N
+*
+*        Generate elementary reflector H(I) to annihilate B(:,I)
+*
+         P = M-L+MIN( L, I )
+         CALL SLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
+         IF( I.LT.N ) THEN
+*
+*           W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)]
+*
+            DO J = 1, N-I
+               T( J, N ) = (A( I, I+J ))
+            END DO
+            CALL SGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB, 
+     $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
+*
+*           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H
+*
+            ALPHA = -(T( I, 1 ))            
+            DO J = 1, N-I
+               A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N ))
+            END DO
+            CALL SGER( P, N-I, ALPHA, B( 1, I ), 1, 
+     $           T( 1, N ), 1, B( 1, I+1 ), LDB )
+         END IF
+      END DO
+*
+      DO I = 2, N
+*
+*        T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I))
+*
+         ALPHA = -T( I, 1 )
+
+         DO J = 1, I-1
+            T( J, I ) = ZERO
+         END DO
+         P = MIN( I-1, L )
+         MP = MIN( M-L+1, M )
+         NP = MIN( P+1, N )
+*
+*        Triangular part of B2
+*
+         DO J = 1, P
+            T( J, I ) = ALPHA*B( M-L+J, I )
+         END DO
+         CALL STRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB,
+     $               T( 1, I ), 1 )
+*
+*        Rectangular part of B2
+*
+         CALL SGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB, 
+     $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
+*
+*        B1
+*
+         CALL SGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, 
+     $               ONE, T( 1, I ), 1 )         
+*
+*        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
+*
+         CALL STRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
+*
+*        T(I,I) = tau(I)
+*
+         T( I, I ) = T( I, 1 )
+         T( I, 1 ) = ZERO
+      END DO
+   
+*
+*     End of STPQRT2
+*
+      END
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
diff --git a/SRC/ztpmqrt.f b/SRC/ztpmqrt.f
new file mode 100644
index 00000000..ffb0ce25
--- /dev/null
+++ b/SRC/ztpmqrt.f
@@ -0,0 +1,264 @@
+      SUBROUTINE ZTPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT,
+     $                    A, LDA, B, LDB, WORK, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 SIDE, TRANS
+      INTEGER   INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16   V( LDV, * ), A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZTPMQRT applies a complex orthogonal matrix Q obtained from a 
+*  "triangular-pentagonal" complex block reflector H to a general
+*  complex matrix C, which consists of two blocks A and B.
+*
+*  Arguments
+*  =========
+*
+*  SIDE    (input) CHARACTER*1
+*          = 'L': apply Q or Q**H from the Left;
+*          = 'R': apply Q or Q**H from the Right.
+*
+*  TRANS   (input) CHARACTER*1
+*          = 'N':  No transpose, apply Q;
+*          = 'C':  Transpose, apply Q**H.
+*
+*  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 number of elementary reflectors whose product defines
+*          the matrix Q.
+*
+*  L       (input) INTEGER
+*          The order of the trapezoidal part of V.  
+*          K >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size used for the storage of T.  K >= NB >= 1.
+*          This must be the same value of NB used to generate T
+*          in CTPQRT.
+*
+*  V       (input) COMPLEX*16 array, dimension (LDA,K)
+*          The i-th column must contain the vector which defines the
+*          elementary reflector H(i), for i = 1,2,...,k, as returned by
+*          CTPQRT in B.  See Further Details.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.
+*          If SIDE = 'L', LDA >= max(1,M);
+*          if SIDE = 'R', LDA >= max(1,N).
+*
+*  T       (input) COMPLEX*16 array, dimension (LDT,K)
+*          The upper triangular factors of the block reflectors
+*          as returned by CTPQRT, stored as a NB-by-K matrix.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  A       (input/output) COMPLEX*16 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 
+*          Q*C or Q**H*C or C*Q or C*Q**H.  See Further 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) COMPLEX*16 array, dimension (LDB,N)
+*          On entry, the M-by-N matrix B.
+*          On exit, B is overwritten by the corresponding block of
+*          Q*C or Q**H*C or C*Q or C*Q**H.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B. 
+*          LDB >= max(1,M).
+*
+*  WORK    (workspace/output) COMPLEX*16 array.  The dimension of WORK is
+*           N*NB if SIDE = 'L', or  M*NB if SIDE = 'R'.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The columns of the pentagonal matrix V contain the elementary reflectors
+*  H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a 
+*  trapezoidal block V2:
+*
+*        V = [V1]
+*            [V2].
+*
+*  The size of the trapezoidal block V2 is determined by the parameter L, 
+*  where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L
+*  rows of a K-by-K upper triangular matrix.  If L=K, V2 is upper triangular;
+*  if L=0, there is no trapezoidal block, hence V = V1 is rectangular.
+*
+*  If SIDE = 'L':  C = [A]  where A is K-by-N,  B is M-by-N and V is M-by-K. 
+*                      [B]   
+*  
+*  If SIDE = 'R':  C = [A B]  where A is M-by-K, B is M-by-N and V is N-by-K.
+*
+*  The complex orthogonal matrix Q is formed from V and T.
+*
+*  If TRANS='N' and SIDE='L', C is on exit replaced with Q * C.
+*
+*  If TRANS='C' and SIDE='L', C is on exit replaced with Q**H * C.
+*
+*  If TRANS='N' and SIDE='R', C is on exit replaced with C * Q.
+*
+*  If TRANS='C' and SIDE='R', C is on exit replaced with C * Q**H.
+*
+*  =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
+      INTEGER            I, IB, MB, LB, KF, Q
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, ZLARFB
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     .. Test the input arguments ..
+*
+      INFO   = 0
+      LEFT   = LSAME( SIDE,  'L' )
+      RIGHT  = LSAME( SIDE,  'R' )
+      TRAN   = LSAME( TRANS, 'C' )
+      NOTRAN = LSAME( TRANS, 'N' )
+*      
+      IF( LEFT ) THEN
+         Q = M
+      ELSE IF ( RIGHT ) THEN
+         Q = N
+      END IF
+      IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
+         INFO = -1
+      ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
+         INFO = -2
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( K.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
+         INFO = -6         
+      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+         INFO = -7
+      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -12
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZTPMQRT', -INFO )
+         RETURN
+      END IF
+*
+*     .. Quick return if possible ..
+*
+      IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN
+*
+      IF( LEFT .AND. TRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF
+            CALL ZTPRFB( 'L', 'C', 'F', 'C', MB, N, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*         
+      ELSE IF( RIGHT .AND. NOTRAN ) THEN
+*
+         DO I = 1, K, NB
+            IB = MIN( NB, K-I+1 )
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL ZTPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, 
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      ELSE IF( LEFT .AND. NOTRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )  
+            MB = MIN( M-L+I+IB-1, M )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-M+L-I+1
+            END IF                   
+            CALL ZTPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( I, 1 ), LDA, B, LDB, WORK, IB )
+         END DO
+*
+      ELSE IF( RIGHT .AND. TRAN ) THEN
+*
+         KF = ((K-1)/NB)*NB+1
+         DO I = KF, 1, -NB
+            IB = MIN( NB, K-I+1 )         
+            MB = MIN( N-L+I+IB-1, N )
+            IF( I.GE.L ) THEN
+               LB = 0
+            ELSE
+               LB = MB-N+L-I+1
+            END IF
+            CALL ZTPRFB( 'R', 'C', 'F', 'C', M, MB, IB, LB,
+     $                   V( 1, I ), LDV, T( 1, I ), LDT, 
+     $                   A( 1, I ), LDA, B, LDB, WORK, M )
+         END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of ZTPMQRT
+*
+      END
diff --git a/SRC/ztpqrt.f b/SRC/ztpqrt.f
new file mode 100644
index 00000000..d0570407
--- /dev/null
+++ b/SRC/ztpqrt.f
@@ -0,0 +1,185 @@
+      SUBROUTINE ZTPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK,
+     $                   INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK routine (version 3.?) --
+*  -- 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 ..
+      INTEGER INFO, LDA, LDB, LDT, N, M, L, NB
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZTPQRT computes a blocked QR factorization of a complex 
+*  "triangular-pentagonal" matrix C, which is composed of a 
+*  triangular block A and pentagonal block B, using the compact 
+*  WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of the
+*          triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  NB      (input) INTEGER
+*          The block size to be used in the blocked QR.  N >= NB >= 1.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) COMPLEX*16 array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) COMPLEX*16 array, dimension (LDT,N)
+*          The upper triangular block reflectors stored in compact form
+*          as a sequence of upper triangular blocks.  See Further Details.
+*          
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= NB.
+*
+*  WORK    (workspace) COMPLEX*16 array, dimension (NB*N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*  
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*
+*  The number of blocks is B = ceiling(N/NB), where each
+*  block is of order NB except for the last block, which is of order 
+*  IB = N - (B-1)*NB.  For each of the B blocks, a upper triangular block
+*  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB 
+*  for the last block) T's are stored in the NB-by-N matrix T as
+*
+*               T = [T1 T2 ... TB].
+*
+* =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      INTEGER    I, IB, LB, MB, IINFO
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL   ZTPQRT2, ZTPRFB, XERBLA
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+         INFO = -4
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.NB ) THEN
+         INFO = -10
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZTPQRT', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( M.EQ.0 .OR. N.EQ.0 ) RETURN
+*
+      DO I = 1, N, NB
+*     
+*     Compute the QR factorization of the current block
+*
+         IB = MIN( N-I+1, NB )
+         MB = MIN( M-L+I+IB-1, M )
+         IF( I.GE.L ) THEN
+            LB = 0
+         ELSE
+            LB = MB-M+L-I+1
+         END IF
+*
+         CALL ZTPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB, 
+     $                 T(1, I ), LDT, IINFO )
+*
+*     Update by applying H**H to B(:,I+IB:N) from the left
+*
+         IF( I+IB.LE.N ) THEN
+            CALL ZTPRFB( 'L', 'C', 'F', 'C', MB, N-I-IB+1, IB, LB,
+     $                    B( 1, I ), LDB, T( 1, I ), LDT, 
+     $                    A( I, I+IB ), LDA, B( 1, I+IB ), LDB, 
+     $                    WORK, IB )
+         END IF
+      END DO
+      RETURN
+*     
+*     End of ZTPQRT
+*
+      END
diff --git a/SRC/ztpqrt2.f b/SRC/ztpqrt2.f
new file mode 100644
index 00000000..0d087fd6
--- /dev/null
+++ b/SRC/ztpqrt2.f
@@ -0,0 +1,224 @@
+      SUBROUTINE ZTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
+      IMPLICIT NONE
+*
+*  -- LAPACK 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 ..
+      INTEGER   INFO, LDA, LDB, LDT, N, M, L
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16   A( LDA, * ), B( LDB, * ), T( LDT, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZTPQRT2 computes a QR factorization of a complex "triangular-pentagonal"
+*  matrix C, which is composed of a triangular block A and pentagonal block B, 
+*  using the compact WY representation for Q.
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The total number of rows of the matrix B.  
+*          M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix B, and the order of
+*          the triangular matrix A.
+*          N >= 0.
+*
+*  L       (input) INTEGER
+*          The number of rows of the upper trapezoidal part of B.  
+*          MIN(M,N) >= L >= 0.  See Further Details.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the upper triangular N-by-N matrix A.
+*          On exit, the elements on and above the diagonal of the array
+*          contain the upper triangular matrix R.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  B       (input/output) COMPLEX*16 array, dimension (LDB,N)
+*          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
+*          are rectangular, and the last L rows are upper trapezoidal.
+*          On exit, B contains the pentagonal matrix V.  See Further Details.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,M).
+*
+*  T       (output) COMPLEX*16 array, dimension (LDT,N)
+*          The N-by-N upper triangular factor T of the block reflector.
+*          See Further Details.
+*
+*  LDT     (input) INTEGER
+*          The leading dimension of the array T.  LDT >= max(1,N)
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*
+*  Further Details
+*  ===============
+*
+*  The input matrix C is a (N+M)-by-N matrix  
+*
+*               C = [ A ]
+*                   [ B ]        
+*
+*  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
+*  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
+*  upper trapezoidal matrix B2:
+*
+*               B = [ B1 ]  <- (M-L)-by-N rectangular
+*                   [ B2 ]  <-     L-by-N upper trapezoidal.
+*
+*  The upper trapezoidal matrix B2 consists of the first L rows of a
+*  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
+*  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
+*
+*  The matrix W stores the elementary reflectors H(i) in the i-th column
+*  below the diagonal (of A) in the (N+M)-by-N input matrix C
+*
+*               C = [ A ]  <- upper triangular N-by-N
+*                   [ B ]  <- M-by-N pentagonal
+*
+*  so that W can be represented as
+*
+*               W = [ I ]  <- identity, N-by-N
+*                   [ V ]  <- M-by-N, same form as B.
+*
+*  Thus, all of information needed for W is contained on exit in B, which
+*  we call V above.  Note that V has the same form as B; that is, 
+*
+*               V = [ V1 ] <- (M-L)-by-N rectangular
+*                   [ V2 ] <-     L-by-N upper trapezoidal.
+*
+*  The columns of V represent the vectors which define the H(i)'s.  
+*  The (M+N)-by-(M+N) block reflector H is then given by
+*
+*               H = I - W * T * W**H
+*
+*  where W**H is the conjugate transpose of W and T is the upper triangular
+*  factor of the block reflector.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16  ONE, ZERO
+      PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER   I, J, P, MP, NP
+      COMPLEX*16   ALPHA
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL  ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+      IF( N.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+         INFO = -7
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZTPQRT2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
+*      
+      DO I = 1, N
+*
+*        Generate elementary reflector H(I) to annihilate B(:,I)
+*
+         P = M-L+MIN( L, I )
+         CALL ZLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
+         IF( I.LT.N ) THEN
+*
+*           W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)]
+*
+            DO J = 1, N-I
+               T( J, N ) = CONJG(A( I, I+J ))
+            END DO
+            CALL ZGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB, 
+     $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
+*
+*           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H
+*
+            ALPHA = -CONJG(T( I, 1 ))            
+            DO J = 1, N-I
+               A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N ))
+            END DO
+            CALL ZGERC( P, N-I, ALPHA, B( 1, I ), 1, 
+     $           T( 1, N ), 1, B( 1, I+1 ), LDB )
+         END IF
+      END DO
+*
+      DO I = 2, N
+*
+*        T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I))
+*
+         ALPHA = -T( I, 1 )
+
+         DO J = 1, I-1
+            T( J, I ) = ZERO
+         END DO
+         P = MIN( I-1, L )
+         MP = MIN( M-L+1, M )
+         NP = MIN( P+1, N )
+*
+*        Triangular part of B2
+*
+         DO J = 1, P
+            T( J, I ) = ALPHA*B( M-L+J, I )
+         END DO
+         CALL ZTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB,
+     $               T( 1, I ), 1 )
+*
+*        Rectangular part of B2
+*
+         CALL ZGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB, 
+     $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
+*
+*        B1
+*
+         CALL ZGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, 
+     $               ONE, T( 1, I ), 1 )         
+*
+*        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
+*
+         CALL ZTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
+*
+*        T(I,I) = tau(I)
+*
+         T( I, I ) = T( I, 1 )
+         T( I, 1 ) = ZERO
+      END DO
+   
+*
+*     End of ZTPQRT2
+*
+      END
diff --git a/SRC/ztprfb.f b/SRC/ztprfb.f
new file mode 100644
index 00000000..c97c35d1
--- /dev/null
+++ b/SRC/ztprfb.f
@@ -0,0 +1,708 @@
+      SUBROUTINE ZTPRFB( 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 ..
+      COMPLEX*16   A( LDA, * ), B( LDB, * ), T( LDT, * ), 
+     $          V( LDV, * ), WORK( LDWORK, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZTPRFB applies a complex "triangular-pentagonal" block reflector H or its 
+*  conjugate transpose H**H to a complex 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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 ..
+      COMPLEX*16   ONE, ZERO
+      PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,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  ZCOPY, ZGEMM, ZLACGV, ZTRMM
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC CONJG
+*     ..
+*     .. 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 ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV,
+     $               WORK, LDWORK )         
+         CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL ZGEMM( 'C', '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 ZTRMM( '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 ZGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ),  LDB )    
+         CALL ZTRMM( '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 ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
+     $               WORK, LDWORK )
+         CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, 
+     $               V, LDV, ONE, WORK, LDWORK )
+         CALL ZGEMM( '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 ZTRMM( '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 ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL ZTRMM( 'R', 'U', 'C', '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 ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, 
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL ZGEMM( 'C', '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 ZTRMM( '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 ZGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, 
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
+     $               WORK, LDWORK, ONE, B,  LDB )
+         CALL ZTRMM( '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 ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, 
+     $               V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL ZGEMM( '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 ZTRMM( '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 ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
+     $               V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
+         CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK,
+     $               V, LDV, ONE, B, LDB )
+         CALL ZTRMM( 'R', 'L', 'C', '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 ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
+     $               WORK, LDB )
+         CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, 
+     $               ONE, WORK, LDWORK )
+         CALL ZGEMM( '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 ZTRMM( '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 ZGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
+     $               ONE, B, LDB )
+         CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, 
+     $               WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL ZTRMM( 'L', 'L', 'C', '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 ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV,
+     $               WORK, LDWORK )
+         CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV,
+     $               ONE, WORK, LDWORK )
+         CALL ZGEMM( 'N', 'C', 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 ZTRMM( '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 ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V, LDV, ONE, B, LDB )
+         CALL ZGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
+     $               V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )   
+         CALL ZTRMM( '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 ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
+     $               WORK( KP, 1 ), LDWORK )
+         CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
+     $               B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
+         CALL ZGEMM( '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 ZTRMM( '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 ZGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
+     $               WORK, LDWORK, ONE, B( MP, 1 ), LDB )
+         CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV, 
+     $               WORK, LDWORK, ONE, B, LDB )
+         CALL ZTRMM( 'L', 'U', 'C', '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 ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV,
+     $               WORK( 1, KP ), LDWORK )
+         CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB,
+     $               V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
+         CALL ZGEMM( 'N', 'C', 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 ZTRMM( '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 ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 
+     $               V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
+         CALL ZGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,    
+     $               V, LDV, ONE, B, LDB )
+         CALL ZTRMM( '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 ZTPRFB
+*
+      END