Changes in TS QR API

20 files changed, 1924 insertions, 1550 deletions

diff --git a/SRC/cgelq.f b/SRC/cgelq.f
index c6c962d7..0abd2d72 100644
--- a/SRC/cgelq.f
+++ b/SRC/cgelq.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE CGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-*                          INFO)
+*       SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+*                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX              A( LDA, * ), WORK1( * ), WORK2( * )
+*       COMPLEX  A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> CGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using CLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise CGELQT:
-*> A = L * Q .
+*> CGELQ computes a LQ factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,8 +42,8 @@
 *>          On exit, the elements on and below the diagonal of the array
 *>          contain the M-by-min(M,N) lower trapezoidal matrix L
 *>          (L is lower triangular if M <= N);
-*>          the elements above the diagonal are the rows of
-*>          blocked V representing Q (see Further Details).
+*>          the elements above the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -56,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is COMPLEX array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    CLASWLQ or CGELQT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): horizontal block size
-*>          WORK1(5): vertical block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           CLASWLQ or CGELQT
+*>          T is COMPLEX array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2),  and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
-*>
+*>          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -114,26 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLASWLQ or DGELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*>  the LQ factorization.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE CGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -142,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      COMPLEX           A( LDA, * ), WORK1( * ), WORK2( * )
+      COMPLEX  A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -173,7 +196,18 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
@@ -184,33 +218,31 @@
         MB = 1
         NB = N
       END IF
-      IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
-      IF( NB.GT.N.OR.NB.LE.M) NB = N
-      MINLW1 = M + 5
-      IF ((NB.GT.M).AND.(N.GT.M)) THEN
-        IF(MOD(N-M, NB-M).EQ.0) THEN
-          NBLCKS = (N-M)/(NB-M)
+      IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+      IF( NB.GT.N .OR. NB.LE.M ) NB = N
+      MINTSZ = M + 5
+      IF ( NB.GT.M .AND. N.GT.M ) THEN
+        IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+          NBLCKS = ( N - M ) / ( NB - M )
         ELSE
-          NBLCKS = (N-M)/(NB-M) + 1
+          NBLCKS = ( N - M ) / ( NB - M ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
-*     Determine if the workspace size satisfies minimum size
+*
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
-     $    .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+     $    .AND. ( .NOT.LQUERY) ) THEN
+        IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             MB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             NB = N
         END IF
-        IF (LWORK2.LT.MB*M) THEN
+        IF ( LWORK.LT.MB*M ) THEN
             LMINWS = .TRUE.
             MB = 1
         END IF
@@ -222,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
-     $   .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = MB*M*NBLCKS+5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = MB * M
-        WORK2(2) = M
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = MB*M*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, MB*M )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGELQ', -INFO )
         RETURN
       ELSE IF (LQUERY) THEN
-       RETURN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The LQ Decomposition
 *
-      IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
-        CALL CGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+      IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+        CALL CGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
       ELSE
-         CALL CLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
-     $                    LWORK2, INFO)
+        CALL CLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+     $                    LWORK, INFO)
       END IF
+      WORK(1) = MAX( 1, MB*M )
       RETURN
 *
 *     End of CGELQ
diff --git a/SRC/cgemlq.f b/SRC/cgemlq.f
index 1a551ca3..03dae76d 100644
--- a/SRC/cgemlq.f
+++ b/SRC/cgemlq.f
@@ -2,17 +2,16 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      COMPLEX        A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      COMPLEX           A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
@@ -21,12 +20,12 @@
 *>     CGEMLQ overwrites the general real M-by-N matrix C with
 *>
 *>
-*>                    SIDE = 'L'     SIDE = 'R'
-*>    TRANS = 'N':      Q * C          C * Q
-*>    TRANS = 'T':      Q**T * C       C * Q**T
-*>    where Q is a complex orthogonal matrix defined as the product
-*>    of blocked elementary reflectors computed by short wide LQ
-*>    factorization (DGELQ)
+*>                      SIDE = 'L'     SIDE = 'R'
+*>      TRANS = 'N':      Q * C          C * Q
+*>      TRANS = 'C':      Q**H * C       C * Q**H
+*>      where Q is a complex unitary matrix defined as the product
+*>      of blocked elementary reflectors computed by short wide
+*>      LQ factorization (CGELQ)
 *> \endverbatim
 *
 *  Arguments:
@@ -62,12 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is COMPLEX array, dimension (LDA,K)
-*>          The i-th row must contain the vector which defines the blocked
-*>          elementary reflector H(i), for i = 1,2,...,k, as returned by
-*>          DLASWLQ in the first k rows of its array argument A.
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,41 +75,42 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is COMPLEX array, dimension (MAX(1,LWORK1)) is
-*>          returned by GEQR.
+*>          T is COMPLEX array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *>          C is COMPLEX array, dimension (LDC,N)
 *>          On entry, the M-by-N matrix C.
 *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
 *> \param[in] LDC
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
+*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -128,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LASWQR or GELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*>  the LQ factorization.
+*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LAMSWLQ or GEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $      C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -157,10 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      COMPLEX           A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+      COMPLEX           A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -174,7 +180,7 @@
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     .. External Subroutines ..
-      EXTERNAL           CLAMSWLQ, CGEMLQT, XERBLA
+      EXTERNAL           ZLAMSWLQ, ZGEMLQT, XERBLA
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -182,26 +188,27 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'C' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
-      IF (LEFT) THEN
+      MB = INT(T(2))
+      NB = INT(T(3))
+      IF ( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
         LW = M * MB
         MN = N
       END IF
-      IF ((NB.GT.K).AND.(MN.GT.K)) THEN
-        IF(MOD(MN-K, NB-K).EQ.0) THEN
-          NBLCKS = (MN-K)/(NB-K)
+*
+      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+          NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
-          NBLCKS = (MN-K)/(NB-K) + 1
+          NBLCKS = ( MN - K ) / ( NB - K ) + 1
         END IF
       ELSE
         NBLCKS = 1
@@ -220,40 +227,43 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+     $         .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
-         INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+        INFO = -11
+      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = REAL( LW )
       END IF
+*
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGEMLQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
-     $   (NB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
         CALL CGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
-     $        WORK1(6), MB, C, LDC, WORK2, INFO)
+     $        T(4), MB, C, LDC, WORK, INFO)
       ELSE
-        CALL CLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $        MB, C, LDC, WORK2, LWORK2, INFO )
+        CALL CLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $    MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = REAL ( LW )
+*
       RETURN
 *
 *     End of CGEMLQ
diff --git a/SRC/cgemqr.f b/SRC/cgemqr.f
index 51d38b85..a5d6263f 100644
--- a/SRC/cgemqr.f
+++ b/SRC/cgemqr.f
@@ -2,17 +2,16 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                     LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      COMPLEX        A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      COMPLEX           A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
@@ -23,8 +22,8 @@
 *>
 *>                      SIDE = 'L'     SIDE = 'R'
 *>      TRANS = 'N':      Q * C          C * Q
-*>      TRANS = 'T':      Q**T * C       C * Q**T
-*>      where Q is a complex orthogonal matrix defined as the product
+*>      TRANS = 'T':      Q**H * C       C * Q**H
+*>      where Q is a complex unitary matrix defined as the product
 *>      of blocked elementary reflectors computed by tall skinny
 *>      QR factorization (CGEQR)
 *> \endverbatim
@@ -62,13 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is COMPLEX array, dimension (LDA,K)
-*>          The i-th column must contain the vector which defines the
-*>          blockedelementary reflector H(i), for i = 1,2,...,k, as
-*>          returned by DGETSQR in the first k columns of
-*>          its array argument A.
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -79,16 +75,16 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is COMPLEX array, dimension (MAX(1,LWORK1)) as
-*>          it is returned by GEQR.
+*>          T is COMPLEX array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
@@ -100,21 +96,21 @@
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
+*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -130,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LATSQR or GEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
+*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LATMSQR or GEMQRT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $        C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -159,11 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      COMPLEX           A( LDA, * ), WORK1( * ), C(LDC, * ),
-     $               WORK2( * )
+      COMPLEX        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -185,14 +188,14 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'C' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
+      MB = INT(T(2))
+      NB = INT(T(3))
       IF(LEFT) THEN
         LW = N * NB
         MN = M
@@ -201,12 +204,12 @@
         MN = N
       END IF
 *
-      IF ((MB.GT.K).AND.(MN.GT.K)) THEN
-          IF(MOD(MN-K, MB-K).EQ.0) THEN
-             NBLCKS = (MN-K)/(MB-K)
-          ELSE
-             NBLCKS = (MN-K)/(MB-K) + 1
-          END IF
+      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+           NBLCKS = ( MN - K ) / ( MB - K )
+        ELSE
+           NBLCKS = ( MN - K ) / ( MB - K ) + 1
+        END IF
       ELSE
         NBLCKS = 1
       END IF
@@ -224,24 +227,25 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+     $        .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
          INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
 *     Determine the block size if it is tall skinny or short and wide
 *
-      IF( INFO.EQ.0)  THEN
-          WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGEMQR', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
        RETURN
       END IF
 *
@@ -251,16 +255,17 @@
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
-     $   (MB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
         CALL CGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        WORK1(6), NB, C, LDC, WORK2, INFO)
+     $        T(4), NB, C, LDC, WORK, INFO )
       ELSE
-        CALL CLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $      NB, C, LDC, WORK2, LWORK2, INFO )
+        CALL CLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $      NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = LW
+*
       RETURN
 *
 *     End of CGEMQR
diff --git a/SRC/cgeqr.f b/SRC/cgeqr.f
index 330fda5c..1743960f 100644
--- a/SRC/cgeqr.f
+++ b/SRC/cgeqr.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE CGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+*       SUBROUTINE CGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
 *                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX           A( LDA, * ), WORK1( * ), WORK2( * )
+*       COMPLEX           A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> CGEQR computes a QR factorization of an M-by-N matrix A,
-*> using CLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise CGEQRT:
-*> A = Q * R .
+*> CGEQR computes a QR factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,7 +42,8 @@
 *>          On exit, the elements on and above the diagonal of the array
 *>          contain the min(M,N)-by-N upper trapezoidal matrix R
 *>          (R is upper triangular if M >= N);
-*>          the elements below the diagonal represent Q (see Further Details).
+*>          the elements below the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -55,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is COMPLEX array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    CLATSQR or CGEQRT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): row block size
-*>          WORK1(5): column block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           CLATSQR or CGEQRT
+*>          T is COMPLEX array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2), and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
+*>          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -113,26 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLATSQR or DGEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE CGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE CGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      COMPLEX           A( LDA, * ), WORK1( * ), WORK2( * )
+      COMPLEX           A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -172,7 +196,18 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
@@ -183,34 +218,31 @@
         MB = M
         NB = 1
       END IF
-      IF( MB.GT.M.OR.MB.LE.N) MB = M
-      IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
-      MINLW1 = N + 5
-      IF ((MB.GT.N).AND.(M.GT.N)) THEN
-        IF(MOD(M-N, MB-N).EQ.0) THEN
-          NBLCKS = (M-N)/(MB-N)
+      IF( MB.GT.M .OR. MB.LE.N ) MB = M
+      IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+      MINTSZ = N + 5
+      IF ( MB.GT.N .AND. M.GT.N ) THEN
+        IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+          NBLCKS = ( M - N ) / ( MB - N )
         ELSE
-          NBLCKS = (M-N)/(MB-N) + 1
+          NBLCKS = ( M - N ) / ( MB - N ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
-     $    .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             NB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             MB = M
         END IF
-        IF (LWORK2.LT.NB*N) THEN
+        IF ( LWORK.LT.NB*N ) THEN
             LMINWS = .TRUE.
             NB = 1
         END IF
@@ -222,45 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
-     $   .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = NB * N * NBLCKS + 5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = NB * N
-        WORK2(2) = N
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = NB*N*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, NB*N )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGEQR', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
-       RETURN
+      ELSE IF ( LQUERY ) THEN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The QR Decomposition
 *
-      IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
-         CALL CGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
-         RETURN
+      IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+        CALL CGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
       ELSE
-         CALL CLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
-     $                    LWORK2, INFO)
+        CALL CLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+     $                    LWORK, INFO )
       END IF
+      WORK(1) = MAX( 1, NB*N )
       RETURN
 *
 *     End of CGEQR
diff --git a/SRC/cgetsls.f b/SRC/cgetsls.f
index af5bd2cb..c2d5368f 100644
--- a/SRC/cgetsls.f
+++ b/SRC/cgetsls.f
@@ -1,16 +1,15 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-*     $                   , WORK, LWORK, INFO )
-
+*       SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+*     $                     WORK, LWORK, INFO )
 *
 *       .. Scalar Arguments ..
 *       CHARACTER          TRANS
 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX   A( LDA, * ), B( LDB, * ), WORK( * )
+*       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
 *       ..
 *
 *
@@ -19,10 +18,11 @@
 *>
 *> \verbatim
 *>
-*> CGETSLS solves overdetermined or underdetermined real linear systems
-*> involving an M-by-N matrix A, or its transpose, using a tall skinny
-*> QR or short wide LQfactorization of A.  It is assumed that A has
-*> full rank.
+*> CGETSLS solves overdetermined or underdetermined complex linear systems
+*> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
+*> factorization of A.  It is assumed that A has full rank.
+*>
+*>
 *>
 *> The following options are provided:
 *>
@@ -80,10 +80,8 @@
 *>          A is COMPLEX array, dimension (LDA,N)
 *>          On entry, the M-by-N matrix A.
 *>          On exit,
-*>            if M >= N, A is overwritten by details of its QR
-*>                       factorization as returned by DGEQRF;
-*>            if M <  N, A is overwritten by details of its LQ
-*>                       factorization as returned by DGELQF.
+*>          A is overwritten by details of its QR or LQ
+*>          factorization as returned by CGEQR or CGEQL.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -97,21 +95,17 @@
 *>          B is COMPLEX array, dimension (LDB,NRHS)
 *>          On entry, the matrix B of right hand side vectors, stored
 *>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
-*>          if TRANS = 'T'.
+*>          if TRANS = 'C'.
 *>          On exit, if INFO = 0, B is overwritten by the solution
 *>          vectors, stored columnwise:
 *>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
-*>          squares solution vectors; the residual sum of squares for the
-*>          solution in each column is given by the sum of squares of
-*>          elements N+1 to M in that column;
+*>          squares solution vectors.
 *>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
 *>          minimum norm solution vectors;
-*>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
+*>          if TRANS = 'C' and m >= n, rows 1 to M of B contain the
 *>          minimum norm solution vectors;
-*>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
-*>          least squares solution vectors; the residual sum of squares
-*>          for the solution in each column is given by the sum of
-*>          squares of elements M+1 to N in that column.
+*>          if TRANS = 'C' and m < n, rows 1 to M of B contain the
+*>          least squares solution vectors.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -122,23 +116,21 @@
 *>
 *> \param[out] WORK
 *> \verbatim
-*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
-*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*>          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
 *>          The dimension of the array WORK.
-*>          LWORK >= max( 1, MN + max( MN, NRHS ) ).
-*>          For optimal performance,
-*>          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
-*>          where MN = min(M,N) and NB is the optimum block size.
-*>
-*>          If LWORK = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK array, returns
-*>          this value as the first entry of the WORK array, and no error
-*>          message related to LWORK is issued by XERBLA.
+*>          If LWORK = -1 or -2, then a workspace query is assumed.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -162,9 +154,11 @@
 *
 *> \date November 2011
 *
+*> \ingroup complexGEsolve
+*
 *  =====================================================================
-      SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-     $                   , WORK, LWORK, INFO )
+      SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+     $                     WORK, LWORK, INFO )
 *
 *  -- LAPACK driver routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -173,7 +167,7 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER          TRANS
-      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, MB
+      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
 *     ..
 *     .. Array Arguments ..
       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
@@ -190,10 +184,11 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY, TRAN
-      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
-     $                   SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2,
-     $                   INFO2, NB
+      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+     $                   SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2,
+     $                   INFO2
       REAL               ANRM, BIGNUM, BNRM, SMLNUM
+      COMPLEX            WRK1(5), WRK2
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -206,7 +201,7 @@
      $                   CTRTRS, XERBLA, CGELQ, CGEMLQ
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          REAL, MAX, MIN
+      INTRINSIC          REAL, MAX, MIN, INT
 *     ..
 *     .. Executable Statements ..
 *
@@ -218,7 +213,7 @@
       MNK   = MAX(MINMN,NRHS)
       TRAN  = LSAME( TRANS, 'C' )
 *
-      LQUERY = ( LWORK.EQ.-1 )
+      LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
       IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
      $    LSAME( TRANS, 'C' ) ) ) THEN
          INFO = -1
@@ -234,53 +229,58 @@
          INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
+      IF( INFO.EQ.0 )  THEN
 *
 *     Determine the block size and minimum LWORK
 *
-      IF ( M.GE.N ) THEN
-        CALL CGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        MB = INT(WORK(4))
-        NB = INT(WORK(5))
-        LW = INT(WORK(6))
-        CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
-      ELSE
-        CALL CGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        MB = INT(WORK(4))
-        NB = INT(WORK(5))
-        LW = INT(WORK(6))
-        CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
-      END IF
+       IF ( M.GE.N ) THEN
+         CALL CGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1 , INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL CGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZM, B, LDB, WRK2, -2, INFO2 )
+         LWM = MAX( LWM, INT(WRK2) )
+       ELSE
+         CALL CGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1, INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL CGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -2, INFO2 )
+         LWM  = MAX( LWM, INT(WRK2) )
+       END IF
 *
-       IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+       IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. ( .NOT.LQUERY ) ) THEN
           INFO=-10
        END IF
       END IF
 *
       IF( INFO.NE.0 ) THEN
-        CALL XERBLA( 'DGETSLS', -INFO )
-        WORK( 1 ) = REAL( WSIZEO )
-        WORK( 2 ) = REAL( WSIZEM )
+        CALL XERBLA( 'CGETSLS', -INFO )
+        WORK( 1 ) = REAL( TSZ + LW )
         RETURN
-      ELSE IF (LQUERY) THEN
-        WORK( 1 ) = REAL( WSIZEO )
-        WORK( 2 ) = REAL( WSIZEM )
+      END IF
+      IF ( LQUERY ) THEN
+        WORK( 1 ) = REAL( TSZ + LW )
         RETURN
       END IF
-      IF(LWORK.LT.WSIZEO) THEN
-        LW1=INT(WORK(3))
-        LW2=MAX(LW,INT(WORK(6)))
+      IF( LWORK.LT.( TSZ + LW ) ) THEN
+        LW1 = TSZM
+        LW2 = LWM
       ELSE
-        LW1=INT(WORK(2))
-        LW2=MAX(LW,INT(WORK(6)))
+        LW1 = TSZ
+        LW2 = LW
       END IF
 *
 *     Quick return if possible
@@ -347,9 +347,9 @@
 *
 *        compute QR factorization of A
 *
-        CALL CGEQR( M, N, A, LDA, WORK(LW2+1), LW1
-     $    , WORK(1), LW2, INFO )
-        IF (.NOT.TRAN) THEN
+        CALL CGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
+     $              WORK(1), LW2, INFO )
+        IF ( .NOT.TRAN ) THEN
 *
 *           Least-Squares Problem min || A * X - B ||
 *
@@ -401,8 +401,8 @@
 *
 *        Compute LQ factorization of A
 *
-        CALL CGELQ( M, N, A, LDA, WORK(LW2+1), LW1
-     $    , WORK(1), LW2, INFO )
+         CALL CGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+     $              WORK(1), LW2, INFO )
 *
 *        workspace at least M, optimally M*NB.
 *
@@ -482,10 +482,9 @@
       END IF
 *
    50 CONTINUE
-       WORK( 1 ) = REAL( WSIZEO )
-       WORK( 2 ) = REAL( WSIZEM )
+      WORK( 1 ) = REAL( TSZ + LW )
       RETURN
 *
-*     End of CGETSLS
+*     End of ZGETSLS
 *
       END
diff --git a/SRC/dgelq.f b/SRC/dgelq.f
index d73f7454..59d9fa91 100644
--- a/SRC/dgelq.f
+++ b/SRC/dgelq.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE DGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-*                          INFO)
+*       SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+*                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION  A( LDA, * ), WORK1( * ), WORK2( * )
+*       DOUBLE PRECISION  A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> DGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using DLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise DGELQT:
-*> A = L * Q .
+*> DGELQ computes a LQ factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,8 +42,8 @@
 *>          On exit, the elements on and below the diagonal of the array
 *>          contain the M-by-min(M,N) lower trapezoidal matrix L
 *>          (L is lower triangular if M <= N);
-*>          the elements above the diagonal are the rows of
-*>          blocked V representing Q (see Further Details).
+*>          the elements above the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -56,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    DLASWLQ or DGELQT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): horizontal block size
-*>          WORK1(5): vertical block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           DLASWLQ or DGELQT
+*>          T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2),  and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
-*>
+*>          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -114,27 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLASWLQ or DGELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*>  the LQ factorization.
 *> \endverbatim
 *>
-*>
 *  =====================================================================
-      SUBROUTINE DGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -143,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION  A( LDA, * ), WORK1( * ), WORK2( * )
+      DOUBLE PRECISION  A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -174,45 +196,53 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
       IF ( MIN(M,N).GT.0 ) THEN
-        MB = ILAENV( 1, 'DGELQ ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'DGELQ ', ' ', M, N, 2, -1)
+        MB = ILAENV( 1, 'DGELQ ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'DGELQ ', ' ', M, N, 2, -1 )
       ELSE
         MB = 1
         NB = N
       END IF
-      IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
-      IF( NB.GT.N.OR.NB.LE.M) NB = N
-      MINLW1 = M + 5
-      IF ((NB.GT.M).AND.(N.GT.M)) THEN
-        IF(MOD(N-M, NB-M).EQ.0) THEN
-          NBLCKS = (N-M)/(NB-M)
+      IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+      IF( NB.GT.N .OR. NB.LE.M ) NB = N
+      MINTSZ = M + 5
+      IF ( NB.GT.M .AND. N.GT.M ) THEN
+        IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+          NBLCKS = ( N - M ) / ( NB - M )
         ELSE
-          NBLCKS = (N-M)/(NB-M) + 1
+          NBLCKS = ( N - M ) / ( NB - M ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
-     $    .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+     $    .AND. ( .NOT.LQUERY) ) THEN
+        IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             MB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             NB = N
         END IF
-        IF (LWORK2.LT.MB*M) THEN
+        IF ( LWORK.LT.MB*M ) THEN
             LMINWS = .TRUE.
             MB = 1
         END IF
@@ -224,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
-     $   .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = MB*M*NBLCKS+5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = MB * M
-        WORK2(2) = M
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = MB*M*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, MB*M )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGELQ', -INFO )
         RETURN
       ELSE IF (LQUERY) THEN
-       RETURN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The LQ Decomposition
 *
-      IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
-        CALL DGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+      IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+        CALL DGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
       ELSE
-        CALL DLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
-     $                    LWORK2, INFO)
+        CALL DLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+     $                    LWORK, INFO)
       END IF
+      WORK(1) = MAX( 1, MB*M )
       RETURN
 *
 *     End of DGELQ
diff --git a/SRC/dgemlq.f b/SRC/dgemlq.f
index 7bdf97a1..17c4de5c 100644
--- a/SRC/dgemlq.f
+++ b/SRC/dgemlq.f
@@ -2,17 +2,16 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      DOUBLE        A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      DOUBLE PRECISION  A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
@@ -62,12 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is DOUBLE PRECISION array, dimension (LDA,K)
-*>          The i-th row must contain the vector which defines the blocked
-*>          elementary reflector H(i), for i = 1,2,...,k, as returned by
-*>          DLASWLQ in the first k rows of its array argument A.
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,41 +75,42 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1)) is
-*>          returned by GEQR.
+*>          T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *>          C is DOUBLE PRECISION array, dimension (LDC,N)
 *>          On entry, the M-by-N matrix C.
 *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
 *> \param[in] LDC
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
+*>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -128,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LASWQR or GELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*>  the LQ factorization.
+*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LAMSWLQ or GEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $      C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -157,10 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+      DOUBLE PRECISION  A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -174,7 +180,7 @@
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     .. External Subroutines ..
-      EXTERNAL           DTPMLQT, DGEMLQT, XERBLA
+      EXTERNAL           DLAMSWLQ, DGEMLQT, XERBLA
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -182,26 +188,27 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = (LWORK2.LT.0)
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'T' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
-      IF (LEFT) THEN
+      MB = INT(T(2))
+      NB = INT(T(3))
+      IF ( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
         LW = M * MB
         MN = N
       END IF
-      IF ((NB.GT.K).AND.(MN.GT.K)) THEN
-        IF(MOD(MN-K, NB-K).EQ.0) THEN
-          NBLCKS = (MN-K)/(NB-K)
+*
+      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+          NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
-          NBLCKS = (MN-K)/(NB-K) + 1
+          NBLCKS = ( MN - K ) / ( NB - K ) + 1
         END IF
       ELSE
         NBLCKS = 1
@@ -220,40 +227,42 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+     $         .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
-         INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+        INFO = -11
+      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = LW
       END IF
+*
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGEMLQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
-     $   (NB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
         CALL DGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
-     $        WORK1(6), MB, C, LDC, WORK2, INFO)
+     $        T(4), MB, C, LDC, WORK, INFO)
       ELSE
-        CALL DLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $    MB, C, LDC, WORK2, LWORK2, INFO )
+        CALL DLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $    MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = LW
 *
       RETURN
 *
diff --git a/SRC/dgemqr.f b/SRC/dgemqr.f
index f47e6a87..ee02d98e 100644
--- a/SRC/dgemqr.f
+++ b/SRC/dgemqr.f
@@ -2,20 +2,16 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                     LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      DOUBLE PRECISION  A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
-*       ..
-*
-*
+*      DOUBLE PRECISION  A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
@@ -65,13 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is DOUBLE PRECISION array, dimension (LDA,K)
-*>          The i-th column must contain the vector which defines the
-*>          blockedelementary reflector H(i), for i = 1,2,...,k, as
-*>          returned by DGETSQR in the first k columns of
-*>          its array argument A.
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -82,16 +75,16 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1)) as
-*>          it is returned by GEQR.
+*>          T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
@@ -103,21 +96,21 @@
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
+*>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -133,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLATSQR or DGEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
+*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LATMSQR or GEMQRT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $        C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -162,11 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A( LDA, * ), WORK1( * ), C(LDC, * ),
-     $               WORK2( * )
+      DOUBLE PRECISION  A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -180,7 +180,7 @@
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     .. External Subroutines ..
-      EXTERNAL           DGEMQRT, DTPMQRT, XERBLA
+      EXTERNAL           DGEMQRT, DLAMTSQR, XERBLA
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,14 +188,14 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'T' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
+      MB = INT(T(2))
+      NB = INT(T(3))
       IF(LEFT) THEN
         LW = N * NB
         MN = M
@@ -204,11 +204,11 @@
         MN = N
       END IF
 *
-      IF ((MB.GT.K).AND.(MN.GT.K)) THEN
-          IF(MOD(MN-K, MB-K).EQ.0) THEN
-             NBLCKS = (MN-K)/(MB-K)
+      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+          IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+             NBLCKS = ( MN - K ) / ( MB - K )
           ELSE
-             NBLCKS = (MN-K)/(MB-K) + 1
+             NBLCKS = ( MN - K ) / ( MB - K ) + 1
           END IF
       ELSE
         NBLCKS = 1
@@ -227,43 +227,44 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+     $        .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ).AND.MIN(M,N,K).NE.0 ) THEN
+      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
          INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
 *     Determine the block size if it is tall skinny or short and wide
 *
-      IF( INFO.EQ.0)  THEN
-         WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGEMQR', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
        RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
-     $   (MB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
         CALL DGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        WORK1(6), NB, C, LDC, WORK2, INFO)
+     $        T(4), NB, C, LDC, WORK, INFO )
       ELSE
-        CALL DLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $      NB, C, LDC, WORK2, LWORK2, INFO )
+        CALL DLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $      NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = LW
 *
       RETURN
 *
diff --git a/SRC/dgeqr.f b/SRC/dgeqr.f
index da0fc4ad..93f02bd8 100644
--- a/SRC/dgeqr.f
+++ b/SRC/dgeqr.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE DGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+*       SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
 *                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION  A( LDA, * ), WORK1( * ), WORK2( * )
+*       DOUBLE PRECISION  A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> DGEQR computes a QR factorization of an M-by-N matrix A,
-*> using DLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise DGEQRT:
-*> A = Q * R .
+*> DGEQR computes a QR factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,7 +42,8 @@
 *>          On exit, the elements on and above the diagonal of the array
 *>          contain the min(M,N)-by-N upper trapezoidal matrix R
 *>          (R is upper triangular if M >= N);
-*>          the elements below the diagonal represent Q (see Further Details).
+*>          the elements below the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -55,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    DLATSQR or DGEQRT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): row block size
-*>          WORK1(5): column block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           DLATSQR or DGEQRT
+*>          T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2), and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
+*>          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -113,26 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLATSQR or DGEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE DGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION  A( LDA, * ), WORK1( * ), WORK2( * )
+      DOUBLE PRECISION  A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -172,45 +196,53 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
       IF ( MIN(M,N).GT.0 ) THEN
-        MB = ILAENV( 1, 'DGEQR ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'DGEQR ', ' ', M, N, 2, -1)
+        MB = ILAENV( 1, 'DGEQR ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'DGEQR ', ' ', M, N, 2, -1 )
       ELSE
         MB = M
         NB = 1
       END IF
-      IF( MB.GT.M.OR.MB.LE.N) MB = M
-      IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
-      MINLW1 = N + 5
-      IF ((MB.GT.N).AND.(M.GT.N)) THEN
-        IF(MOD(M-N, MB-N).EQ.0) THEN
-          NBLCKS = (M-N)/(MB-N)
+      IF( MB.GT.M .OR. MB.LE.N ) MB = M
+      IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+      MINTSZ = N + 5
+      IF ( MB.GT.N .AND. M.GT.N ) THEN
+        IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+          NBLCKS = ( M - N ) / ( MB - N )
         ELSE
-          NBLCKS = (M-N)/(MB-N) + 1
+          NBLCKS = ( M - N ) / ( MB - N ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
-     $    .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             NB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             MB = M
         END IF
-        IF (LWORK2.LT.NB*N) THEN
+        IF ( LWORK.LT.NB*N ) THEN
             LMINWS = .TRUE.
             NB = 1
         END IF
@@ -222,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
-     $   .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
-
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = NB * N * NBLCKS + 5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = NB * N
-        WORK2(2) = N
+*
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = NB*N*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, NB*N )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGEQR', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
-       RETURN
+      ELSE IF ( LQUERY ) THEN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The QR Decomposition
 *
-      IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
-         CALL DGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
+      IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+        CALL DGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
       ELSE
-         CALL DLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
-     $                    LWORK2, INFO)
+        CALL DLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+     $                    LWORK, INFO )
       END IF
+      WORK(1) = MAX( 1, NB*N )
       RETURN
 *
 *     End of DGEQR
diff --git a/SRC/dgetsls.f b/SRC/dgetsls.f
index f2a7b3a7..0a94694e 100644
--- a/SRC/dgetsls.f
+++ b/SRC/dgetsls.f
@@ -29,16 +29,21 @@
 *> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
 *>    an overdetermined system, i.e., solve the least squares problem
 *>                 minimize || B - A*X ||.
-
+*>
 *> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
 *>    an underdetermined system A * X = B.
-
+*>
 *> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
 *>    an undetermined system A**T * X = B.
-
+*>
 *> 4. If TRANS = 'T' and m < n:  find the least squares solution of
 *>    an overdetermined system, i.e., solve the least squares problem
 *>                 minimize || B - A**T * X ||.
+*>
+*> Several right hand side vectors b and solution vectors x can be
+*> handled in a single call; they are stored as the columns of the
+*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
+*> matrix X.
 *> \endverbatim
 *
 *  Arguments:
@@ -76,7 +81,7 @@
 *>          On entry, the M-by-N matrix A.
 *>          On exit,
 *>          A is overwritten by details of its QR or LQ
-*>          factorization as returned by DGETSQR.
+*>          factorization as returned by DGEQR or DGEQL.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -111,18 +116,21 @@
 *>
 *> \param[out] WORK
 *> \verbatim
-*>          WORK is DOUBLE PRECISION array, dimension (MAX(12,LWORK))
-*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK,
-*>          and WORK(2) returns the minimum LWORK.
+*>          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
 *>          The dimension of the array WORK.
-*>          IF LWORK=-1, workspace query is assumed, and
-*>          WORK(1) returns the optimal LWORK,
-*>          and WORK(2) returns the minimum LWORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -174,10 +182,10 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY, TRAN
-      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
-     $                   SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2,
+      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+     $                   SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
      $                   INFO2
-      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMLNUM
+      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMLNUM, WRK1(5), WRK2
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -190,7 +198,7 @@
      $                   DTRTRS, XERBLA, DGELQ, DGEMLQ
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          DBLE, MAX, MIN
+      INTRINSIC          DBLE, MAX, MIN, INT
 *     ..
 *     .. Executable Statements ..
 *
@@ -199,10 +207,10 @@
       INFO=0
       MINMN = MIN( M, N )
       MAXMN = MAX( M, N )
-      MNK   = MAX(MINMN,NRHS)
+      MNK   = MAX( MINMN, NRHS )
       TRAN  = LSAME( TRANS, 'T' )
 *
-      LQUERY = ( LWORK.EQ.-1 )
+      LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
       IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
      $    LSAME( TRANS, 'T' ) ) ) THEN
          INFO = -1
@@ -218,49 +226,58 @@
          INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
+      IF( INFO.EQ.0 )  THEN
 *
 *     Determine the block size and minimum LWORK
 *
        IF ( M.GE.N ) THEN
-        CALL DGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        LW = INT(WORK(6))
-        CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+         CALL DGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1 , INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL DGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZM, B, LDB, WRK2, -2, INFO2 )
+         LWM = MAX( LWM, INT(WRK2) )
        ELSE
-        CALL DGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        LW = INT(WORK(6))
-        CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+         CALL DGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1, INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL DGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -2, INFO2 )
+         LWM  = MAX( LWM, INT(WRK2) )
        END IF
 *
-       IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+       IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
           INFO=-10
        END IF
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGETSLS', -INFO )
-        WORK( 1 ) = DBLE( WSIZEO )
-        WORK( 2 ) = DBLE( WSIZEM )
+        WORK( 1 ) = DBLE( TSZ + LW )
         RETURN
-      ELSE IF (LQUERY) THEN
-        WORK( 1 ) = DBLE( WSIZEO )
-        WORK( 2 ) = DBLE( WSIZEM )
+      END IF
+      IF ( LQUERY ) THEN
+        WORK( 1 ) = DBLE( TSZ + LW )
         RETURN
       END IF
-      IF(LWORK.LT.WSIZEO) THEN
-        LW1=INT(WORK(3))
-        LW2=MAX(LW,INT(WORK(6)))
+      IF( LWORK.LT.( TSZ + LW ) ) THEN
+        LW1 = TSZM
+        LW2 = LWM
       ELSE
-        LW1=INT(WORK(2))
-        LW2=MAX(LW,INT(WORK(6)))
+        LW1 = TSZ
+        LW2 = LW
       END IF
 *
 *     Quick return if possible
@@ -323,13 +340,13 @@
          IBSCL = 2
       END IF
 *
-      IF ( M.GE.N) THEN
+      IF ( M.GE.N ) THEN
 *
 *        compute QR factorization of A
 *
         CALL DGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
      $              WORK(1), LW2, INFO )
-        IF (.NOT.TRAN) THEN
+        IF ( .NOT.TRAN ) THEN
 *
 *           Least-Squares Problem min || A * X - B ||
 *
@@ -381,7 +398,7 @@
 *
 *        Compute LQ factorization of A
 *
-        CALL DGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+         CALL DGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
      $              WORK(1), LW2, INFO )
 *
 *        workspace at least M, optimally M*NB.
@@ -454,16 +471,15 @@
      $                INFO )
       END IF
       IF( IBSCL.EQ.1 ) THEN
-         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+        CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
       ELSE IF( IBSCL.EQ.2 ) THEN
-         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+        CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
      $                INFO )
       END IF
 *
    50 CONTINUE
-       WORK( 1 ) = DBLE( WSIZEO )
-       WORK( 2 ) = DBLE( WSIZEM )
+      WORK( 1 ) = DBLE( TSZ + LW )
       RETURN
 *
 *     End of DGETSLS
diff --git a/SRC/sgelq.f b/SRC/sgelq.f
index 8a759834..adc606d9 100644
--- a/SRC/sgelq.f
+++ b/SRC/sgelq.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE SGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-*                          INFO)
+*       SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+*                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       REAL              A( LDA, * ), WORK1( * ), WORK2( * )
+*       REAL  A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> SGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using SLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise SGELQT:
-*> A = L * Q .
+*> SGELQ computes a LQ factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,8 +42,8 @@
 *>          On exit, the elements on and below the diagonal of the array
 *>          contain the M-by-min(M,N) lower trapezoidal matrix L
 *>          (L is lower triangular if M <= N);
-*>          the elements above the diagonal are the rows of
-*>          blocked V representing Q (see Further Details).
+*>          the elements above the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -56,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is REAL array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    SLASWLQ or SGELQT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): horizontal block size
-*>          WORK1(5): vertical block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           SLASWLQ or SGELQT
+*>          T is REAL array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2),  and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) REAL array, dimension (MAX(1,LWORK2))
-*>
+*>          (workspace) REAL array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -114,27 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLASWLQ or DGELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*>  the LQ factorization.
 *> \endverbatim
 *>
-*>
 *  =====================================================================
-      SUBROUTINE SGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -143,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      REAL              A( LDA, * ), WORK1( * ), WORK2( * )
+      REAL              A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -174,7 +196,18 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
@@ -185,34 +218,31 @@
         MB = 1
         NB = N
       END IF
-      IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
-      IF( NB.GT.N.OR.NB.LE.M) NB = N
-      MINLW1 = M + 5
-      IF ((NB.GT.M).AND.(N.GT.M)) THEN
-        IF(MOD(N-M, NB-M).EQ.0) THEN
-          NBLCKS = (N-M)/(NB-M)
+      IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+      IF( NB.GT.N .OR. NB.LE.M ) NB = N
+      MINTSZ = M + 5
+      IF ( NB.GT.M .AND. N.GT.M ) THEN
+        IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+          NBLCKS = ( N - M ) / ( NB - M )
         ELSE
-          NBLCKS = (N-M)/(NB-M) + 1
+          NBLCKS = ( N - M ) / ( NB - M ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
-     $    .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+     $    .AND. ( .NOT.LQUERY) ) THEN
+        IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             MB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             NB = N
         END IF
-        IF (LWORK2.LT.MB*M) THEN
+        IF ( LWORK.LT.MB*M ) THEN
             LMINWS = .TRUE.
             MB = 1
         END IF
@@ -224,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
-     $   .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = MB*M*NBLCKS+5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = MB * M
-        WORK2(2) = M
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = MB*M*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, MB*M )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGELQ', -INFO )
         RETURN
       ELSE IF (LQUERY) THEN
-       RETURN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The LQ Decomposition
 *
-      IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
-        CALL SGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+      IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+        CALL SGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
       ELSE
-         CALL SLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
-     $                    LWORK2, INFO)
+        CALL SLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+     $                    LWORK, INFO)
       END IF
+      WORK(1) = MAX( 1, MB*M )
       RETURN
 *
 *     End of SGELQ
diff --git a/SRC/sgemlq.f b/SRC/sgemlq.f
index 14a37a4d..a9cd54bd 100644
--- a/SRC/sgemlq.f
+++ b/SRC/sgemlq.f
@@ -2,23 +2,22 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      REAL              A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      REAL              A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
 *> \verbatim
 *>
-*>     DGEMLQ overwrites the general real M-by-N matrix C with
+*>      SGEMLQ overwrites the general real M-by-N matrix C with
 *>
 *>
 *>                    SIDE = 'L'     SIDE = 'R'
@@ -26,7 +25,7 @@
 *>    TRANS = 'T':      Q**T * C       C * Q**T
 *>    where Q is a real orthogonal matrix defined as the product
 *>    of blocked elementary reflectors computed by short wide LQ
-*>    factorization (DGELQ)
+*>    factorization (SGELQ)
 *> \endverbatim
 *
 *  Arguments:
@@ -62,12 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is REAL array, dimension (LDA,K)
-*>          The i-th row must contain the vector which defines the blocked
-*>          elementary reflector H(i), for i = 1,2,...,k, as returned by
-*>          DLASWLQ in the first k rows of its array argument A.
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,41 +75,42 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is REAL array, dimension (MAX(1,LWORK1)) is
-*>          returned by GEQR.
+*>          T is REAL array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *>          C is REAL array, dimension (LDC,N)
 *>          On entry, the M-by-N matrix C.
 *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
 *> \param[in] LDC
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) REAL array, dimension (MAX(1,LWORK2))
+*>         (workspace) REAL array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -128,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LASWLQ or GELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*>  the LQ factorization.
+*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LAMSWLQ or GEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $      C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -157,10 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      REAL              A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+      REAL              A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -174,7 +180,7 @@
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     .. External Subroutines ..
-      EXTERNAL           STPMLQT, SGEMLQT, XERBLA
+      EXTERNAL           SLAMSWLQ, SGEMLQT, XERBLA
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -182,26 +188,27 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'T' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
-      IF (LEFT) THEN
+      MB = INT(T(2))
+      NB = INT(T(3))
+      IF ( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
         LW = M * MB
         MN = N
       END IF
-      IF ((NB.GT.K).AND.(MN.GT.K)) THEN
-        IF(MOD(MN-K, NB-K).EQ.0) THEN
-          NBLCKS = (MN-K)/(NB-K)
+*
+      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+          NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
-          NBLCKS = (MN-K)/(NB-K) + 1
+          NBLCKS = ( MN - K ) / ( NB - K ) + 1
         END IF
       ELSE
         NBLCKS = 1
@@ -220,40 +227,43 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+     $         .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
-         INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+        INFO = -11
+      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = REAL(LW)
       END IF
+*
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGEMLQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
-     $   (NB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
         CALL SGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
-     $        WORK1(6), MB, C, LDC, WORK2, INFO)
+     $        T(4), MB, C, LDC, WORK, INFO)
       ELSE
-        CALL SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $    MB, C, LDC, WORK2, LWORK2, INFO )
+        CALL SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $    MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = REAL(LW)
+*
       RETURN
 *
 *     End of SGEMLQ
diff --git a/SRC/sgemqr.f b/SRC/sgemqr.f
index cda7990c..78eb75c8 100644
--- a/SRC/sgemqr.f
+++ b/SRC/sgemqr.f
@@ -2,17 +2,16 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                     LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      REAL              A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      REAL              A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
@@ -26,7 +25,7 @@
 *>      TRANS = 'T':      Q**T * C       C * Q**T
 *>      where Q is a real orthogonal matrix defined as the product
 *>      of blocked elementary reflectors computed by tall skinny
-*>      QR factorization (DGEQR)
+*>      QR factorization (SGEQR)
 *> \endverbatim
 *
 *  Arguments:
@@ -62,13 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is REAL array, dimension (LDA,K)
-*>          The i-th column must contain the vector which defines the
-*>          blockedelementary reflector H(i), for i = 1,2,...,k, as
-*>          returned by DGETSQR in the first k columns of
-*>          its array argument A.
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -79,16 +75,16 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is REAL array, dimension (MAX(1,LWORK1)) as
-*>          it is returned by GEQR.
+*>          T is REAL array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
@@ -100,21 +96,21 @@
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) REAL array, dimension (MAX(1,LWORK2))
+*>         (workspace) REAL array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -130,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LATSQR or GEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
+*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LATMSQR or GEMQRT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $        C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -159,11 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      REAL              A( LDA, * ), WORK1( * ), C(LDC, * ),
-     $               WORK2( * )
+      REAL              A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -177,7 +180,7 @@
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     .. External Subroutines ..
-      EXTERNAL           SGEMQRT, STPMQRT, XERBLA
+      EXTERNAL           SGEMQRT, SLAMTSQR, XERBLA
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -185,14 +188,14 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'T' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
+      MB = INT(T(2))
+      NB = INT(T(3))
       IF(LEFT) THEN
         LW = N * NB
         MN = M
@@ -201,11 +204,11 @@
         MN = N
       END IF
 *
-      IF ((MB.GT.K).AND.(MN.GT.K)) THEN
-          IF(MOD(MN-K, MB-K).EQ.0) THEN
-             NBLCKS = (MN-K)/(MB-K)
+      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+          IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+             NBLCKS = ( MN - K ) / ( MB - K )
           ELSE
-             NBLCKS = (MN-K)/(MB-K) + 1
+             NBLCKS = ( MN - K ) / ( MB - K ) + 1
           END IF
       ELSE
         NBLCKS = 1
@@ -224,43 +227,44 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+     $        .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ).AND.MIN(M,N,K).NE.0 ) THEN
+      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
          INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
 *     Determine the block size if it is tall skinny or short and wide
 *
-      IF( INFO.EQ.0)  THEN
-        WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGEMQR', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
        RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
-     $   (MB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
         CALL SGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        WORK1(6), NB, C, LDC, WORK2, INFO)
+     $        T(4), NB, C, LDC, WORK, INFO )
       ELSE
-        CALL SLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $      NB, C, LDC, WORK2, LWORK2, INFO )
+        CALL SLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $      NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = LW
 *
       RETURN
 *
diff --git a/SRC/sgeqr.f b/SRC/sgeqr.f
index 41e04622..ad18e8a2 100644
--- a/SRC/sgeqr.f
+++ b/SRC/sgeqr.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE SGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+*       SUBROUTINE SGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
 *                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       REAL              A( LDA, * ), WORK1( * ), WORK2( * )
+*       REAL              A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> SGEQR computes a QR factorization of an M-by-N matrix A,
-*> using SLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise SGEQRT:
-*> A = Q * R .
+*> SGEQR computes a QR factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,7 +42,8 @@
 *>          On exit, the elements on and above the diagonal of the array
 *>          contain the min(M,N)-by-N upper trapezoidal matrix R
 *>          (R is upper triangular if M >= N);
-*>          the elements below the diagonal represent Q (see Further Details).
+*>          the elements below the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -55,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is REAL array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    DLATSQR or DGEQRT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): row block size
-*>          WORK1(5): column block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           SLATSQR or SGEQRT
+*>          T is REAL array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2), and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) REAL array, dimension (MAX(1,LWORK2))
+*>          (workspace) REAL array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -113,26 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLATSQR or DGEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE SGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE SGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      REAL              A( LDA, * ), WORK1( * ), WORK2( * )
+      REAL              A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -172,7 +196,18 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
@@ -183,34 +218,31 @@
         MB = M
         NB = 1
       END IF
-      IF( MB.GT.M.OR.MB.LE.N) MB = M
-      IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
-      MINLW1 = N + 5
-      IF ((MB.GT.N).AND.(M.GT.N)) THEN
-        IF(MOD(M-N, MB-N).EQ.0) THEN
-          NBLCKS = (M-N)/(MB-N)
+      IF( MB.GT.M .OR. MB.LE.N ) MB = M
+      IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+      MINTSZ = N + 5
+      IF ( MB.GT.N .AND. M.GT.N ) THEN
+        IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+          NBLCKS = ( M - N ) / ( MB - N )
         ELSE
-          NBLCKS = (M-N)/(MB-N) + 1
+          NBLCKS = ( M - N ) / ( MB - N ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
-     $    .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             NB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             MB = M
         END IF
-        IF (LWORK2.LT.NB*N) THEN
+        IF ( LWORK.LT.NB*N ) THEN
             LMINWS = .TRUE.
             NB = 1
         END IF
@@ -222,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
-     $   .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
-
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = NB * N * NBLCKS + 5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = NB * N
-        WORK2(2) = N
+*
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = NB*N*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, NB*N )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGEQR', -INFO )
         RETURN
       ELSE IF (LQUERY) THEN
-       RETURN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The QR Decomposition
 *
-      IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
-         CALL SGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
+      IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+        CALL SGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
       ELSE
-         CALL SLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
-     $                    LWORK2, INFO)
+        CALL SLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+     $                    LWORK, INFO )
       END IF
+      WORK(1) = MAX( 1, NB*N )
       RETURN
 *
 *     End of SGEQR
diff --git a/SRC/sgetsls.f b/SRC/sgetsls.f
index b7bcd0f0..dcbaf36a 100644
--- a/SRC/sgetsls.f
+++ b/SRC/sgetsls.f
@@ -1,16 +1,15 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-*     $                   , WORK, LWORK, INFO )
-
+*       SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+*     $                     WORK, LWORK, INFO )
 *
 *       .. Scalar Arguments ..
 *       CHARACTER          TRANS
 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
 *       ..
 *       .. Array Arguments ..
-*       REAL   A( LDA, * ), B( LDB, * ), WORK( * )
+*       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
 *       ..
 *
 *
@@ -77,7 +76,7 @@
 *>          On entry, the M-by-N matrix A.
 *>          On exit,
 *>          A is overwritten by details of its QR or LQ
-*>          factorization as returned by DGETSQR.
+*>          factorization as returned by SGEQR or SGEQL.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -112,18 +111,21 @@
 *>
 *> \param[out] WORK
 *> \verbatim
-*>          WORK is REAL array, dimension (MAX(1,LWORK))
-*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK,
-*>          and WORK(2) returns the minimum LWORK.
+*>          (workspace) REAL array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
 *>          The dimension of the array WORK.
-*>          IF LWORK=-1, workspace query is assumed, and
-*>          WORK(1) returns the optimal LWORK,
-*>          and WORK(2) returns the minimum LWORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -150,8 +152,8 @@
 *> \ingroup doubleGEsolve
 *
 *  =====================================================================
-      SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-     $                   , WORK, LWORK, INFO )
+      SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+     $                     WORK, LWORK, INFO )
 *
 *  -- LAPACK driver routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -160,7 +162,7 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER          TRANS
-      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, MB
+      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
 *     ..
 *     .. Array Arguments ..
       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
@@ -175,10 +177,11 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY, TRAN
-      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
-     $                   SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2, INFO2,
-     $                   NB
+      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+     $                   SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
+     $                   INFO2
       REAL               ANRM, BIGNUM, BNRM, SMLNUM
+      COMPLEX            WRK1(5), WRK2
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -191,7 +194,7 @@
      $                   STRTRS, XERBLA, SGELQ, SGEMLQ
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          REAL, MAX, MIN
+      INTRINSIC          REAL, MAX, MIN, INT
 *     ..
 *     .. Executable Statements ..
 *
@@ -200,10 +203,10 @@
       INFO=0
       MINMN = MIN( M, N )
       MAXMN = MAX( M, N )
-      MNK   = MAX(MINMN,NRHS)
+      MNK   = MAX( MINMN, NRHS )
       TRAN  = LSAME( TRANS, 'T' )
 *
-      LQUERY = ( LWORK.EQ.-1 )
+      LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
       IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
      $    LSAME( TRANS, 'T' ) ) ) THEN
          INFO = -1
@@ -219,53 +222,58 @@
          INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
+      IF( INFO.EQ.0 )  THEN
 *
 *     Determine the block size and minimum LWORK
 *
        IF ( M.GE.N ) THEN
-        CALL SGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        MB = INT(WORK(4))
-        NB = INT(WORK(5))
-        LW = INT(WORK(6))
-        CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+         CALL SGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1 , INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL SGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZM, B, LDB, WRK2, -2, INFO2 )
+         LWM = MAX( LWM, INT(WRK2) )
        ELSE
-        CALL SGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        MB = INT(WORK(4))
-        NB = INT(WORK(5))
-        LW = INT(WORK(6))
-        CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+         CALL SGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1, INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL SGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -2, INFO2 )
+         LWM  = MAX( LWM, INT(WRK2) )
        END IF
 *
-       IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+       IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
           INFO=-10
        END IF
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGETSLS', -INFO )
-        WORK( 1 ) = REAL( WSIZEO )
-        WORK( 2 ) = REAL( WSIZEM )
+        WORK( 1 ) = REAL( TSZ + LW )
         RETURN
-      ELSE IF (LQUERY) THEN
-        WORK( 1 ) = REAL( WSIZEO )
-        WORK( 2 ) = REAL( WSIZEM )
+      END IF
+      IF ( LQUERY ) THEN
+        WORK( 1 ) = REAL( TSZ + LW )
         RETURN
       END IF
-      IF(LWORK.LT.WSIZEO) THEN
-        LW1=INT(WORK(3))
-        LW2=MAX(LW,INT(WORK(6)))
+      IF( LWORK.LT.( TSZ + LW ) ) THEN
+        LW1 = TSZM
+        LW2 = LWM
       ELSE
-        LW1=INT(WORK(2))
-        LW2=MAX(LW,INT(WORK(6)))
+        LW1 = TSZ
+        LW2 = LW
       END IF
 *
 *     Quick return if possible
@@ -332,9 +340,9 @@
 *
 *        compute QR factorization of A
 *
-        CALL SGEQR( M, N, A, LDA, WORK(LW2+1), LW1
-     $    , WORK(1), LW2, INFO )
-        IF (.NOT.TRAN) THEN
+        CALL SGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
+     $              WORK(1), LW2, INFO )
+        IF ( .NOT.TRAN ) THEN
 *
 *           Least-Squares Problem min || A * X - B ||
 *
@@ -386,8 +394,8 @@
 *
 *        Compute LQ factorization of A
 *
-        CALL SGELQ( M, N, A, LDA, WORK(LW2+1), LW1
-     $    , WORK(1), LW2, INFO )
+         CALL SGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+     $              WORK(1), LW2, INFO )
 *
 *        workspace at least M, optimally M*NB.
 *
@@ -467,8 +475,7 @@
       END IF
 *
    50 CONTINUE
-       WORK( 1 ) = REAL( WSIZEO )
-       WORK( 2 ) = REAL( WSIZEM )
+      WORK( 1 ) = REAL( TSZ + LW )
       RETURN
 *
 *     End of SGETSLS
diff --git a/SRC/zgelq.f b/SRC/zgelq.f
index 33125b3d..5c51cf52 100644
--- a/SRC/zgelq.f
+++ b/SRC/zgelq.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE CGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-*                          INFO)
+*       SUBROUTINE ZGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+*                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX*16        A( LDA, * ), WORK1( * ), WORK2( * )
+*       COMPLEX*16  A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> ZGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using ZLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise ZGELQT:
-*> A = L * Q .
+*> ZGELQ computes a LQ factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,8 +42,8 @@
 *>          On exit, the elements on and below the diagonal of the array
 *>          contain the M-by-min(M,N) lower trapezoidal matrix L
 *>          (L is lower triangular if M <= N);
-*>          the elements above the diagonal are the rows of
-*>          blocked V representing Q (see Further Details).
+*>          the elements above the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -56,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    ZLASWLQ or ZGELQT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): horizontal block size
-*>          WORK1(5): vertical block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           ZLASWLQ or ZGELQT
+*>          T is COMPLEX*16 array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2),  and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
-*>
+*>          (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -114,26 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLASWLQ or DGELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*>  the LQ factorization.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE ZGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE ZGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -142,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16  A( LDA, * ), WORK1( * ), WORK2( * )
+      COMPLEX*16  A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -173,7 +196,18 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
@@ -184,34 +218,31 @@
         MB = 1
         NB = N
       END IF
-      IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
-      IF( NB.GT.N.OR.NB.LE.M) NB = N
-      MINLW1 = M + 5
-      IF ((NB.GT.M).AND.(N.GT.M)) THEN
-        IF(MOD(N-M, NB-M).EQ.0) THEN
-          NBLCKS = (N-M)/(NB-M)
+      IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+      IF( NB.GT.N .OR. NB.LE.M ) NB = N
+      MINTSZ = M + 5
+      IF ( NB.GT.M .AND. N.GT.M ) THEN
+        IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+          NBLCKS = ( N - M ) / ( NB - M )
         ELSE
-          NBLCKS = (N-M)/(NB-M) + 1
+          NBLCKS = ( N - M ) / ( NB - M ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
-     $    .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+     $    .AND. ( .NOT.LQUERY) ) THEN
+        IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             MB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             NB = N
         END IF
-        IF (LWORK2.LT.MB*M) THEN
+        IF ( LWORK.LT.MB*M ) THEN
             LMINWS = .TRUE.
             MB = 1
         END IF
@@ -223,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
-     $   .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = MB*M*NBLCKS+5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = MB * M
-        WORK2(2) = M
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = MB*M*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, MB*M )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGELQ', -INFO )
         RETURN
       ELSE IF (LQUERY) THEN
-       RETURN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The LQ Decomposition
 *
-      IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
-        CALL ZGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+      IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+        CALL ZGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
       ELSE
-         CALL ZLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
-     $                    LWORK2, INFO)
+        CALL ZLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+     $                    LWORK, INFO)
       END IF
+      WORK(1) = MAX( 1, MB*M )
       RETURN
 *
 *     End of ZGELQ
diff --git a/SRC/zgemlq.f b/SRC/zgemlq.f
index 10d3a5e4..f02d7b1a 100644
--- a/SRC/zgemlq.f
+++ b/SRC/zgemlq.f
@@ -2,31 +2,30 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      COMPLEX*16        A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      COMPLEX*16        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
 *> \verbatim
 *>
-*>     ZGEMLQ overwrites the general real M-by-N matrix C with
+*>      ZGEMLQ overwrites the general real M-by-N matrix C with
 *>
 *>
-*>                    SIDE = 'L'     SIDE = 'R'
-*>    TRANS = 'N':      Q * C          C * Q
-*>    TRANS = 'T':      Q**T * C       C * Q**T
-*>    where Q is a complex orthogonal matrix defined as the product
-*>    of blocked elementary reflectors computed by short wide LQ
-*>    factorization (DGELQ)
+*>                      SIDE = 'L'     SIDE = 'R'
+*>      TRANS = 'N':      Q * C          C * Q
+*>      TRANS = 'C':      Q**H * C       C * Q**H
+*>      where Q is a complex unitary matrix defined as the product
+*>      of blocked elementary reflectors computed by short wide
+*>      LQ factorization (ZGELQ)
 *> \endverbatim
 *
 *  Arguments:
@@ -62,12 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is COMPLEX*16 array, dimension (LDA,K)
-*>          The i-th row must contain the vector which defines the blocked
-*>          elementary reflector H(i), for i = 1,2,...,k, as returned by
-*>          DLASWLQ in the first k rows of its array argument A.
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,41 +75,42 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1)) is
-*>          returned by GEQR.
+*>          T is COMPLEX*16 array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *>          C is COMPLEX*16 array, dimension (LDC,N)
 *>          On entry, the M-by-N matrix C.
 *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
 *> \param[in] LDC
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
+*>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -128,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LASWLQ or GELQT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GELQ will use either
-*>  LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*>  the LQ decomposition.
-*>  The output of LASWLQ or GELQT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LASWLQ or GELQT was used is the same as used below in
-*>  GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LASWLQ or GELQT.
+*>  LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*>  the LQ factorization.
+*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LAMSWLQ or GEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $      C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -157,10 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16     A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+      COMPLEX*16        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -182,26 +188,27 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'C' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
-      IF (LEFT) THEN
+      MB = INT(T(2))
+      NB = INT(T(3))
+      IF ( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
         LW = M * MB
         MN = N
       END IF
-      IF ((NB.GT.K).AND.(MN.GT.K)) THEN
-        IF(MOD(MN-K, NB-K).EQ.0) THEN
-          NBLCKS = (MN-K)/(NB-K)
+*
+      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+          NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
-          NBLCKS = (MN-K)/(NB-K) + 1
+          NBLCKS = ( MN - K ) / ( NB - K ) + 1
         END IF
       ELSE
         NBLCKS = 1
@@ -220,40 +227,43 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+     $         .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
-         INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+        INFO = -11
+      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0)  THEN
-        WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = LW
       END IF
+*
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGEMLQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
-     $   (NB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
         CALL ZGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
-     $        WORK1(6), MB, C, LDC, WORK2, INFO)
+     $        T(4), MB, C, LDC, WORK, INFO)
       ELSE
-        CALL ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $    MB, C, LDC, WORK2, LWORK2, INFO )
+        CALL ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $    MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = LW
+*
       RETURN
 *
 *     End of ZGEMLQ
diff --git a/SRC/zgemqr.f b/SRC/zgemqr.f
index 3141067f..71614d5e 100644
--- a/SRC/zgemqr.f
+++ b/SRC/zgemqr.f
@@ -2,17 +2,16 @@
 *  Definition:
 *  ===========
 *
-*      SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-*     $                     LWORK1, C, LDC, WORK2, LWORK2, INFO )
+*      SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
 *      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      COMPLEX*16        A( LDA, * ), WORK1( * ), C(LDC, * ),
-*     $                  WORK2( * )
+*      COMPLEX*16        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *> \par Purpose:
 *  =============
 *>
@@ -23,8 +22,8 @@
 *>
 *>                      SIDE = 'L'     SIDE = 'R'
 *>      TRANS = 'N':      Q * C          C * Q
-*>      TRANS = 'T':      Q**T * C       C * Q**T
-*>      where Q is a complex orthogonal matrix defined as the product
+*>      TRANS = 'T':      Q**H * C       C * Q**H
+*>      where Q is a complex unitary matrix defined as the product
 *>      of blocked elementary reflectors computed by tall skinny
 *>      QR factorization (ZGEQR)
 *> \endverbatim
@@ -62,13 +61,10 @@
 *>
 *> \endverbatim
 *>
-*> \param[in,out] A
+*> \param[in] A
 *> \verbatim
 *>          A is COMPLEX*16 array, dimension (LDA,K)
-*>          The i-th column must contain the vector which defines the
-*>          blockedelementary reflector H(i), for i = 1,2,...,k, as
-*>          returned by DGETSQR in the first k columns of
-*>          its array argument A.
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -79,16 +75,16 @@
 *>          if SIDE = 'R', LDA >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK1
+*> \param[in] T
 *> \verbatim
-*>          WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1)) as
-*>          it is returned by GEQR.
+*>          T is COMPLEX*16 array, dimension (MAX(5,TSIZE)).
+*>          Part of the data structure to represent Q as returned by SGEQR.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
+*>          TSIZE is INTEGER
+*>          The dimension of the array T. TSIZE >= 5.
 *> \endverbatim
 *>
 *> \param[in,out] C
@@ -100,21 +96,21 @@
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
+*>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
 *>
 *> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK2 array, returns
-*>          this value as the third entry of the WORK2 array (WORK2(1)),
-*>          and no error message related to LWORK2 is issued by XERBLA.
-*>
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1, then a workspace query is assumed. The routine
+*>          only calculates the size of the WORK array, returns this
+*>          value as WORK(1), and no error message related to WORK 
+*>          is issued by XERBLA.
 *> \endverbatim
+*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -130,27 +126,35 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           LATSQR or GEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
+*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  multiply matrix Q by another matrix.
+*>  Further Details in LATMSQR or GEMQRT.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
-     $        C, LDC, WORK2, LWORK2, INFO )
+      SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+     $        C, LDC, WORK, LWORK, INFO )
 *
 *  -- LAPACK computational routine (version 3.5.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -159,11 +163,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16        A( LDA, * ), WORK1( * ), C(LDC, * ),
-     $               WORK2( * )
+      COMPLEX*16        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
@@ -185,14 +188,14 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK2.LT.0
+      LQUERY  = LWORK.LT.0
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'C' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(WORK1(4))
-      NB = INT(WORK1(5))
+      MB = INT(T(2))
+      NB = INT(T(3))
       IF(LEFT) THEN
         LW = N * NB
         MN = M
@@ -201,11 +204,11 @@
         MN = N
       END IF
 *
-      IF ((MB.GT.K).AND.(MN.GT.K)) THEN
-          IF(MOD(MN-K, MB-K).EQ.0) THEN
-             NBLCKS = (MN-K)/(MB-K)
+      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+          IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+             NBLCKS = ( MN - K ) / ( MB - K )
           ELSE
-             NBLCKS = (MN-K)/(MB-K) + 1
+             NBLCKS = ( MN - K ) / ( MB - K ) + 1
           END IF
       ELSE
         NBLCKS = 1
@@ -224,45 +227,47 @@
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+     $        .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
          INFO = -11
-      ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
 *     Determine the block size if it is tall skinny or short and wide
 *
-      IF( INFO.EQ.0)  THEN
-         WORK2(1) = LW
+      IF( INFO.EQ.0 )  THEN
+        WORK(1) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGEMQR', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF ( LQUERY ) THEN
        RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
-      IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
-     $   (MB.GE.MAX(M,N,K))) THEN
+      IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+     $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
         CALL ZGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        WORK1(6), NB, C, LDC, WORK2, INFO)
+     $        T(4), NB, C, LDC, WORK, INFO )
       ELSE
-        CALL ZLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
-     $      NB, C, LDC, WORK2, LWORK2, INFO )
+        CALL ZLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+     $      NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK2(1) = LW
+      WORK(1) = LW
+*
       RETURN
 *
-*     End of DGEMQR
+*     End of ZGEMQR
 *
       END
diff --git a/SRC/zgeqr.f b/SRC/zgeqr.f
index 10fab97f..226db33e 100644
--- a/SRC/zgeqr.f
+++ b/SRC/zgeqr.f
@@ -2,14 +2,14 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE ZGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+*       SUBROUTINE ZGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
 *                        INFO)
 *
 *       .. Scalar Arguments ..
-*       INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+*       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX*16        A( LDA, * ), WORK1( * ), WORK2( * )
+*       COMPLEX*16        A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -17,11 +17,7 @@
 *  =============
 *>
 *> \verbatim
-*>
-*> ZGEQR computes a QR factorization of an M-by-N matrix A,
-*> using ZLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise ZGEQRT:
-*> A = Q * R .
+*> ZGEQR computes a QR factorization of an M-by-N matrix A.
 *> \endverbatim
 *
 *  Arguments:
@@ -46,7 +42,8 @@
 *>          On exit, the elements on and above the diagonal of the array
 *>          contain the min(M,N)-by-N upper trapezoidal matrix R
 *>          (R is upper triangular if M >= N);
-*>          the elements below the diagonal represent Q (see Further Details).
+*>          the elements below the diagonal are used to store part of the 
+*>          data structure to represent Q.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -55,47 +52,50 @@
 *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
-*> \param[out] WORK1
+*> \param[out] T
 *> \verbatim
-*>          WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1))
-*>          WORK1 contains part of the data structure used to store Q.
-*>          WORK1(1): algorithm type = 1, to indicate output from
-*>                    ZLATSQR or ZGEQRT
-*>          WORK1(2): optimum size of WORK1
-*>          WORK1(3): minimum size of WORK1
-*>          WORK1(4): row block size
-*>          WORK1(5): column block size
-*>          WORK1(6:LWORK1): data structure needed for Q, computed by
-*>                           CLATSQR or CGEQRT
+*>          T is COMPLEX*16 array, dimension (MAX(5,TSIZE))
+*>          On exit, if INFO = 0, T(1) returns optimal (or either minimal 
+*>          or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*>          Remaining T contains part of the data structure used to represent Q.
+*>          If one wants to apply or construct Q, then one needs to keep T 
+*>          (in addition to A) and pass it to further subroutines.
 *> \endverbatim
 *>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
 *> \verbatim
-*>          LWORK1 is INTEGER
-*>          The dimension of the array WORK1.
-*>          If LWORK1 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK1 and
-*>          returns this value in WORK1(2), and calculates the minimum
-*>          size of WORK1 and returns this value in WORK1(3).
-*>          No error message related to LWORK1 is issued by XERBLA when
-*>          LWORK1 = -1.
+*>          TSIZE is INTEGER
+*>          If TSIZE >= 5, the dimension of the array T.
+*>          If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If TSIZE = -1, the routine calculates optimal size of T for the 
+*>          optimum performance and returns this value in T(1).
+*>          If TSIZE = -2, the routine calculates minimal size of T and 
+*>          returns this value in T(1).
 *> \endverbatim
 *>
-*> \param[out] WORK2
+*> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
+*>          (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
-*> \param[in] LWORK2
+*> \param[in] LWORK
 *> \verbatim
-*>          LWORK2 is INTEGER
-*>          The dimension of the array WORK2.
-*>          If LWORK2 = -1, then a query is assumed. In this case the
-*>          routine calculates the optimal size of WORK2 and
-*>          returns this value in WORK2(1), and calculates the minimum
-*>          size of WORK2 and returns this value in WORK2(2).
-*>          No error message related to LWORK2 is issued by XERBLA when
-*>          LWORK2 = -1.
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.
+*>          If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*>          only calculates the sizes of the T and WORK arrays, returns these
+*>          values as the first entries of the T and WORK arrays, and no error
+*>          message related to T or WORK is issued by XERBLA.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -113,26 +113,50 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \par Further Details:
-*  =====================
+*> \par Further Details
+*  ====================
 *>
 *> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular 
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one 
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However, 
+*> for a given execution, the size of T and WORK are fixed and will not change 
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+*  ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not 
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*>          T(2): row block size (MB)
+*>          T(3): column block size (NB)
+*>          T(4:TSIZE): data structure needed for Q, computed by
+*>                           DLATSQR or DGEQRT
+*>
 *>  Depending on the matrix dimensions M and N, and row and column
 *>  block sizes MB and NB returned by ILAENV, GEQR will use either
 *>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*>  the QR decomposition.
-*>  The output of LATSQR or GEQRT representing Q is stored in A and in
-*>  array WORK1(6:LWORK1) for later use.
-*>  WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*>  which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*>  WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*>  decide whether LATSQR or GEQRT was used is the same as used below in
-*>  GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*>  Further Details in LATSQR or GEQRT.
+*>  the QR factorization.
 *> \endverbatim
 *>
 *  =====================================================================
-      SUBROUTINE ZGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+      SUBROUTINE ZGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
      $   INFO)
 *
 *  -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
 *     November 2013
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, LWORK1, LWORK2
+      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16      A( LDA, * ), WORK1( * ), WORK2( * )
+      COMPLEX*16        A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS
-      INTEGER    MB, NB, I, II, KK, MINLW1, NBLCKS
+      LOGICAL    LQUERY, LMINWS, MINT, MINW
+      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
 *     ..
 *     .. EXTERNAL FUNCTIONS ..
       LOGICAL            LSAME
@@ -172,7 +196,18 @@
 *
       INFO = 0
 *
-      LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+     $           LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+      MINT = .FALSE.
+      IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINT = .TRUE.
+      ENDIF
+*
+      MINW = .FALSE.
+      IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+        MINW = .TRUE.
+      ENDIF
 *
 *     Determine the block size
 *
@@ -183,34 +218,31 @@
         MB = M
         NB = 1
       END IF
-      IF( MB.GT.M.OR.MB.LE.N) MB = M
-      IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
-      MINLW1 = N + 5
-      IF ((MB.GT.N).AND.(M.GT.N)) THEN
-        IF(MOD(M-N, MB-N).EQ.0) THEN
-          NBLCKS = (M-N)/(MB-N)
+      IF( MB.GT.M .OR. MB.LE.N ) MB = M
+      IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+      MINTSZ = N + 5
+      IF ( MB.GT.N .AND. M.GT.N ) THEN
+        IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+          NBLCKS = ( M - N ) / ( MB - N )
         ELSE
-          NBLCKS = (M-N)/(MB-N) + 1
+          NBLCKS = ( M - N ) / ( MB - N ) + 1
         END IF
       ELSE
         NBLCKS = 1
       END IF
 *
-*     Determine if the workspace size satisfies minimum size
+*     Determine if the workspace size satisfies minimal size
 *
       LMINWS = .FALSE.
-      IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
-     $    .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
-     $    .AND.(.NOT.LQUERY)) THEN
-        IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+      IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             NB = 1
-        END IF
-        IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
-            LMINWS = .TRUE.
             MB = M
         END IF
-        IF (LWORK2.LT.NB*N) THEN
+        IF ( LWORK.LT.NB*N ) THEN
             LMINWS = .TRUE.
             NB = 1
         END IF
@@ -222,44 +254,50 @@
         INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
-      ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
-     $   .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+      ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+     $   .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -6
-      ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
-     $   .AND.(.NOT.LMINWS)) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+     $   .AND. ( .NOT.LMINWS ) ) THEN
         INFO = -8
       END IF
-
-      IF( INFO.EQ.0)  THEN
-        WORK1(1) = 1
-        WORK1(2) = NB * N * NBLCKS + 5
-        WORK1(3) = MINLW1
-        WORK1(4) = MB
-        WORK1(5) = NB
-        WORK2(1) = NB * N
-        WORK2(2) = N
+*
+      IF( INFO.EQ.0 )  THEN
+        IF ( MINT ) THEN
+          T(1) = MINTSZ
+        ELSE
+          T(1) = NB*N*NBLCKS + 5
+        ENDIF
+        T(2) = MB
+        T(3) = NB
+        IF ( MINW ) THEN
+          WORK(1) = MAX( 1, N )
+        ELSE
+          WORK(1) = MAX( 1, NB*N )
+        ENDIF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGEQR', -INFO )
         RETURN
       ELSE IF (LQUERY) THEN
-       RETURN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN(M,N).EQ.0 ) THEN
-          RETURN
+        RETURN
       END IF
 *
 *     The QR Decomposition
 *
-      IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
-         CALL ZGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
+      IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+        CALL ZGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
       ELSE
-         CALL ZLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
-     $                    LWORK2, INFO)
+        CALL ZLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+     $                    LWORK, INFO )
       END IF
+      WORK(1) = MAX( 1, NB*N )
       RETURN
 *
 *     End of ZGEQR
diff --git a/SRC/zgetsls.f b/SRC/zgetsls.f
index d61b88c3..2a46fc1b 100644
--- a/SRC/zgetsls.f
+++ b/SRC/zgetsls.f
@@ -1,16 +1,15 @@
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-*     $                   , WORK, LWORK, INFO )
-
+*       SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+*     $                     WORK, LWORK, INFO )
 *
 *       .. Scalar Arguments ..
 *       CHARACTER          TRANS
 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX*16   A( LDA, * ), B( LDB, * ), WORK( * )
+*       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 *       ..
 *
 *
@@ -19,10 +18,11 @@
 *>
 *> \verbatim
 *>
-*> ZGETSLS solves overdetermined or underdetermined real linear systems
-*> involving an M-by-N matrix A, or its transpose, using a tall skinny
-*> QR or short wide LQfactorization of A.  It is assumed that A has
-*> full rank.
+*> ZGETSLS solves overdetermined or underdetermined complex linear systems
+*> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
+*> factorization of A.  It is assumed that A has full rank.
+*>
+*>
 *>
 *> The following options are provided:
 *>
@@ -80,10 +80,8 @@
 *>          A is COMPLEX*16 array, dimension (LDA,N)
 *>          On entry, the M-by-N matrix A.
 *>          On exit,
-*>            if M >= N, A is overwritten by details of its QR
-*>                       factorization as returned by DGEQRF;
-*>            if M <  N, A is overwritten by details of its LQ
-*>                       factorization as returned by DGELQF.
+*>          A is overwritten by details of its QR or LQ
+*>          factorization as returned by ZGEQR or ZGEQL.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -97,21 +95,17 @@
 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 *>          On entry, the matrix B of right hand side vectors, stored
 *>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
-*>          if TRANS = 'T'.
+*>          if TRANS = 'C'.
 *>          On exit, if INFO = 0, B is overwritten by the solution
 *>          vectors, stored columnwise:
 *>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
-*>          squares solution vectors; the residual sum of squares for the
-*>          solution in each column is given by the sum of squares of
-*>          elements N+1 to M in that column;
+*>          squares solution vectors.
 *>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
 *>          minimum norm solution vectors;
-*>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
+*>          if TRANS = 'C' and m >= n, rows 1 to M of B contain the
 *>          minimum norm solution vectors;
-*>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
-*>          least squares solution vectors; the residual sum of squares
-*>          for the solution in each column is given by the sum of
-*>          squares of elements M+1 to N in that column.
+*>          if TRANS = 'C' and m < n, rows 1 to M of B contain the
+*>          least squares solution vectors.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -122,23 +116,21 @@
 *>
 *> \param[out] WORK
 *> \verbatim
-*>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
-*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*>          (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*>          or optimal, if query was assumed) LWORK.
+*>          See LWORK for details.
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
 *>          The dimension of the array WORK.
-*>          LWORK >= max( 1, MN + max( MN, NRHS ) ).
-*>          For optimal performance,
-*>          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
-*>          where MN = min(M,N) and NB is the optimum block size.
-*>
-*>          If LWORK = -1, then a workspace query is assumed; the routine
-*>          only calculates the optimal size of the WORK array, returns
-*>          this value as the first entry of the WORK array, and no error
-*>          message related to LWORK is issued by XERBLA.
+*>          If LWORK = -1 or -2, then a workspace query is assumed.
+*>          If LWORK = -1, the routine calculates optimal size of WORK for the
+*>          optimal performance and returns this value in WORK(1).
+*>          If LWORK = -2, the routine calculates minimal size of WORK and 
+*>          returns this value in WORK(1).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -162,9 +154,11 @@
 *
 *> \date November 2011
 *
+*> \ingroup complex16GEsolve
+*
 *  =====================================================================
-      SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-     $                   , WORK, LWORK, INFO )
+      SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+     $                     WORK, LWORK, INFO )
 *
 *  -- LAPACK driver routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -173,10 +167,10 @@
 *
 *     .. Scalar Arguments ..
       CHARACTER          TRANS
-      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, MB, NB
+      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16   A( LDA, * ), B( LDB, * ), WORK( * )
+      COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 *
 *     ..
 *
@@ -190,9 +184,11 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LQUERY, TRAN
-      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
-     $                   SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2, INFO2
+      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+     $                   SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
+     $                   INFO2
       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMLNUM
+      COMPLEX*16         WRK1(5), WRK2
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -205,7 +201,7 @@
      $                   ZTRTRS, XERBLA, ZGELQ, ZGEMLQ
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          DBLE, MAX, MIN
+      INTRINSIC          DBLE, MAX, MIN, INT
 *     ..
 *     .. Executable Statements ..
 *
@@ -217,7 +213,7 @@
       MNK   = MAX(MINMN,NRHS)
       TRAN  = LSAME( TRANS, 'C' )
 *
-      LQUERY = ( LWORK.EQ.-1 )
+      LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
       IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
      $    LSAME( TRANS, 'C' ) ) ) THEN
          INFO = -1
@@ -233,53 +229,58 @@
          INFO = -8
       END IF
 *
-      IF( INFO.EQ.0)  THEN
+      IF( INFO.EQ.0 )  THEN
 *
 *     Determine the block size and minimum LWORK
 *
-      IF ( M.GE.N ) THEN
-        CALL ZGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        MB = INT(WORK(4))
-        NB = INT(WORK(5))
-        LW = INT(WORK(6))
-        CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
-      ELSE
-        CALL ZGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
-     $   INFO2)
-        MB = INT(WORK(4))
-        NB = INT(WORK(5))
-        LW = INT(WORK(6))
-        CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
-     $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
-        WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
-        WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
-      END IF
+       IF ( M.GE.N ) THEN
+         CALL ZGEQR( M, N, A, LDA, WRK1, -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL ZGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1 , INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+     $        TSZM, B, LDB, WRK2, -2, INFO2 )
+         LWM = MAX( LWM, INT(WRK2) )
+       ELSE
+         CALL ZGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+         TSZ = INT( WRK1(1) )
+         LW  = INT( WRK2 )
+         CALL ZGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+         TSZM = INT( WRK1(1) )
+         LWM  = INT( WRK2 )
+         CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -1, INFO2 )
+         LW  = MAX( LW, INT(WRK2) )
+         CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+     $        TSZ, B, LDB, WRK2, -2, INFO2 )
+         LWM  = MAX( LWM, INT(WRK2) )
+       END IF
 *
-       IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+       IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
           INFO=-10
        END IF
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGETSLS', -INFO )
-        WORK( 1 ) = DBLE( WSIZEO )
-        WORK( 2 ) = DBLE( WSIZEM )
+        WORK( 1 ) = DBLE( TSZ + LW )
         RETURN
-      ELSE IF (LQUERY) THEN
-        WORK( 1 ) = DBLE( WSIZEO )
-        WORK( 2 ) = DBLE( WSIZEM )
+      END IF
+      IF ( LQUERY ) THEN
+        WORK( 1 ) = DBLE( TSZ + LW )
         RETURN
       END IF
-      IF(LWORK.LT.WSIZEO) THEN
-        LW1=INT(WORK(3))
-        LW2=MAX(LW,INT(WORK(6)))
+      IF( LWORK.LT.( TSZ + LW ) ) THEN
+        LW1 = TSZM
+        LW2 = LWM
       ELSE
-        LW1=INT(WORK(2))
-        LW2=MAX(LW,INT(WORK(6)))
+        LW1 = TSZ
+        LW2 = LW
       END IF
 *
 *     Quick return if possible
@@ -346,9 +347,9 @@
 *
 *        compute QR factorization of A
 *
-        CALL ZGEQR( M, N, A, LDA, WORK(LW2+1), LW1
-     $    , WORK(1), LW2, INFO )
-        IF (.NOT.TRAN) THEN
+        CALL ZGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
+     $              WORK(1), LW2, INFO )
+        IF ( .NOT.TRAN ) THEN
 *
 *           Least-Squares Problem min || A * X - B ||
 *
@@ -400,8 +401,8 @@
 *
 *        Compute LQ factorization of A
 *
-        CALL ZGELQ( M, N, A, LDA, WORK(LW2+1), LW1
-     $    , WORK(1), LW2, INFO )
+         CALL ZGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+     $              WORK(1), LW2, INFO )
 *
 *        workspace at least M, optimally M*NB.
 *
@@ -481,8 +482,7 @@
       END IF
 *
    50 CONTINUE
-       WORK( 1 ) = DBLE( WSIZEO )
-       WORK( 2 ) = DBLE( WSIZEM )
+      WORK( 1 ) = DBLE( TSZ + LW )
       RETURN
 *
 *     End of ZGETSLS