Merge pull request #101 from karturov/master

TS QR: changed API, added LAPACKE interfaces and fixes.

16 files changed, 602 insertions, 262 deletions

diff --git a/TESTING/LIN/cchkaa.f b/TESTING/LIN/cchkaa.f
index f2ef59f1..5881043f 100644
--- a/TESTING/LIN/cchkaa.f
+++ b/TESTING/LIN/cchkaa.f
@@ -1047,8 +1047,7 @@
             CALL CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A( 1, 1 ),
      $                   A( 1, 2 ), A( 1, 3 ), A( 1, 4 ), A( 1, 5 ),
-     $                   S( 1 ), S( NMAX+1 ), WORK, RWORK, IWORK,
-     $                   NOUT )
+     $                   S( 1 ), S( NMAX+1 ), NOUT )
          ELSE
             WRITE( NOUT, FMT = 9989 )PATH
          END IF
diff --git a/TESTING/LIN/cdrvls.f b/TESTING/LIN/cdrvls.f
index ededde5a..49c99354 100644
--- a/TESTING/LIN/cdrvls.f
+++ b/TESTING/LIN/cdrvls.f
@@ -10,8 +10,7 @@
 *
 *       SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
 *                          NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-*                          COPYB, C, S, COPYS, WORK, RWORK, IWORK,
-*                          NOUT )
+*                          COPYB, C, S, COPYS, NOUT )
 *
 *       .. Scalar Arguments ..
 *       LOGICAL            TSTERR
@@ -20,11 +19,10 @@
 *       ..
 *       .. Array Arguments ..
 *       LOGICAL            DOTYPE( * )
-*       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+*       INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
 *      $                   NVAL( * ), NXVAL( * )
-*       REAL               COPYS( * ), RWORK( * ), S( * )
-*       COMPLEX            A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-*      $                   WORK( * )
+*       REAL               COPYS( * ), S( * )
+*       COMPLEX            A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
 *       ..
 *
 *
@@ -33,7 +31,7 @@
 *>
 *> \verbatim
 *>
-*> CDRVLS tests the least squares driver routines CGELS, CGELSS, CGELSY
+*> CDRVLS tests the least squares driver routines CGELS, CGETSLS, CGELSS, CGELSY
 *> and CGELSD.
 *> \endverbatim
 *
@@ -171,22 +169,6 @@
 *>                      (min(MMAX,NMAX))
 *> \endverbatim
 *>
-*> \param[out] WORK
-*> \verbatim
-*>          WORK is COMPLEX array, dimension
-*>                      (MMAX*NMAX + 4*NMAX + MMAX).
-*> \endverbatim
-*>
-*> \param[out] RWORK
-*> \verbatim
-*>          RWORK is REAL array, dimension (5*NMAX-1)
-*> \endverbatim
-*>
-*> \param[out] IWORK
-*> \verbatim
-*>          IWORK is INTEGER array, dimension (15*NMAX)
-*> \endverbatim
-*>
 *> \param[in] NOUT
 *> \verbatim
 *>          NOUT is INTEGER
@@ -208,8 +190,7 @@
 *  =====================================================================
       SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-     $                   COPYB, C, S, COPYS, WORK, RWORK, IWORK,
-     $                   NOUT )
+     $                   COPYB, C, S, COPYS, NOUT )
 *
 *  -- LAPACK test routine (version 3.6.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -223,11 +204,10 @@
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
-      INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+      INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
      $                   NVAL( * ), NXVAL( * )
-      REAL               COPYS( * ), RWORK( * ), S( * )
-      COMPLEX            A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-     $                   WORK( * )
+      REAL               COPYS( * ), S( * )
+      COMPLEX            A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
 *     ..
 *
 *  =====================================================================
@@ -249,12 +229,22 @@
       INTEGER            CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
      $                   ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
      $                   LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
-     $                   NFAIL, NRHS, NROWS, NRUN, RANK, MB, LWTS
+     $                   NFAIL, NRHS, NROWS, NRUN, RANK, MB,
+     $                   MMAX, NMAX, NSMAX, LIWORK, LRWORK,
+     $                   LWORK_CGELS, LWORK_CGETSLS, LWORK_CGELSS,
+     $                   LWORK_CGELSY, LWORK_CGELSD,
+     $                   LRWORK_CGELSY, LRWORK_CGELSS, LRWORK_CGELSD
       REAL               EPS, NORMA, NORMB, RCOND
 *     ..
 *     .. Local Arrays ..
-      INTEGER            ISEED( 4 ), ISEEDY( 4 )
-      REAL               RESULT( NTESTS )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 ), IWORKQUERY
+      REAL               RESULT( NTESTS ), RWORKQUERY
+      COMPLEX            WORKQUERY
+*     ..
+*     .. Allocatable Arrays ..
+      COMPLEX, ALLOCATABLE :: WORK (:)
+      REAL, ALLOCATABLE :: RWORK (:)
+      INTEGER, ALLOCATABLE :: IWORK (:)
 *     ..
 *     .. External Functions ..
       REAL               CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
@@ -267,7 +257,7 @@
      $                   SAXPY, XLAENV
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          MAX, MIN, REAL, SQRT
+      INTRINSIC          MAX, MIN, INT, REAL, SQRT
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -311,6 +301,77 @@
      $   CALL ALAHD( NOUT, PATH )
       INFOT = 0
 *
+*     Compute maximal workspace needed for all routines
+*
+      NMAX = 0
+      MMAX = 0
+      NSMAX = 0
+      DO I = 1, NM
+         IF ( MVAL( I ).GT.MMAX ) THEN
+            MMAX = MVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NN
+         IF ( NVAL( I ).GT.NMAX ) THEN
+            NMAX = NVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NNS
+         IF ( NSVAL( I ).GT.NSMAX ) THEN
+            NSMAX = NSVAL( I )
+         END IF
+      ENDDO
+      M = MMAX
+      N = NMAX
+      NRHS = NSMAX
+      LDA = MAX( 1, M )
+      LDB = MAX( 1, M, N )
+      MNMIN = MAX( MIN( M, N ), 1 )
+*
+*     Compute workspace needed for routines
+*     CQRT14, CQRT17 (two side cases), CQRT15 and CQRT12
+*
+      LWORK = MAX( ( M+N )*NRHS,
+     $      ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ),
+     $      MAX( M+MNMIN, NRHS*MNMIN,2*N+M ),
+     $      MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) )
+*
+*     Compute workspace needed for CGELS
+      CALL CGELS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $            WORKQUERY, -1, INFO )
+      LWORK_CGELS = INT( WORKQUERY )
+*     Compute workspace needed for CGETSLS
+      CALL CGETSLS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $              WORKQUERY, -1, INFO )
+      LWORK_CGETSLS = INT( WORKQUERY )
+*     Compute workspace needed for CGELSY
+      CALL CGELSY( M, N, NRHS, A, LDA, B, LDB, IWORKQUERY,
+     $             RCOND, CRANK, WORKQUERY, -1, RWORK, INFO )
+      LWORK_CGELSY = INT( WORKQUERY )
+      LRWORK_CGELSY = 2*N
+*     Compute workspace needed for CGELSS
+      CALL CGELSS( M, N, NRHS, A, LDA, B, LDB, S,
+     $             RCOND, CRANK, WORKQUERY, -1, RWORK, INFO )
+      LWORK_CGELSS = INT( WORKQUERY )
+      LRWORK_CGELSS = 5*MNMIN
+*     Compute workspace needed for CGELSD
+      CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, CRANK,
+     $             WORKQUERY, -1, RWORKQUERY, IWORKQUERY, INFO )
+      LWORK_CGELSD = INT( WORKQUERY )
+      LRWORK_CGELSD = INT( RWORKQUERY )
+*     Compute LIWORK workspace needed for CGELSY and CGELSD
+      LIWORK = MAX( 1, N, IWORKQUERY )
+*     Compute LRWORK workspace needed for CGELSY, CGELSS and CGELSD
+      LRWORK = MAX( 1, LRWORK_CGELSY, LRWORK_CGELSS, LRWORK_CGELSD )
+*     Compute LWORK workspace needed for all functions
+      LWORK = MAX( 1, LWORK, LWORK_CGELS, LWORK_CGETSLS, LWORK_CGELSY,
+     $             LWORK_CGELSS, LWORK_CGELSD )
+      LWLSY = LWORK
+*
+      ALLOCATE( WORK( LWORK ) )
+      ALLOCATE( IWORK( LIWORK ) )
+      ALLOCATE( RWORK( LRWORK ) )
+*
       DO 140 IM = 1, NM
          M = MVAL( IM )
          LDA = MAX( 1, M )
@@ -320,16 +381,9 @@
             MNMIN = MAX(MIN( M, N ),1)
             LDB = MAX( 1, M, N )
             MB = (MNMIN+1)
-            IF(MNMIN.NE.MB) THEN
-              LWTS = (((LDB-MB)/(MB-MNMIN))*MNMIN+MNMIN*2)*MB+5
-            ELSE
-              LWTS = 2*MNMIN+5
-            END IF
 *
             DO 120 INS = 1, NNS
                NRHS = NSVAL( INS )
-               LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
-     $                 M*N+4*MNMIN+MAX( M, N ), 2*N+M, LWTS )
 *
                DO 110 IRANK = 1, 2
                   DO 100 ISCALE = 1, 3
@@ -580,12 +634,6 @@
                            IWORK( J ) = 0
    70                   CONTINUE
 *
-*                       Set LWLSY to the adequate value.
-*
-                        LWLSY = MNMIN + MAX( 2*MNMIN, NB*( N+1 ),
-     $                          MNMIN+NB*NRHS )
-                        LWLSY = MAX( 1, LWLSY )
-*
                         SRNAMT = 'CGELSY'
                         CALL CGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
      $                               RCOND, CRANK, WORK, LWLSY, RWORK,
@@ -776,6 +824,10 @@
  9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
      $      ', MB=', I4,', NB=', I4,', type', I2,
      $      ', test(', I2, ')=', G12.5 )
+*
+      DEALLOCATE( WORK )
+      DEALLOCATE( RWORK )
+      DEALLOCATE( IWORK )
       RETURN
 *
 *     End of CDRVLS
diff --git a/TESTING/LIN/cerrtsqr.f b/TESTING/LIN/cerrtsqr.f
index b8b42dcc..399fcc41 100644
--- a/TESTING/LIN/cerrtsqr.f
+++ b/TESTING/LIN/cerrtsqr.f
@@ -161,13 +161,13 @@
       CALL CGEMQR( 'L', 'N', 2, 1, 0, A, 0, TAU, 1, C, 1, W, 1,INFO)
       CALL CHKXER( 'CGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL CGEMQR( 'R', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL CGEMQR( 'R', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'CGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL CGEMQR( 'L', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL CGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'CGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 11
-      CALL CGEMQR( 'L', 'N', 2, 1, 1, A, 1, TAU, 6, C, 0, W, 1,INFO)
+      CALL CGEMQR( 'L', 'N', 2, 1, 1, A, 2, TAU, 6, C, 0, W, 1,INFO)
       CALL CHKXER( 'CGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 13
       CALL CGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 6, C, 2, W, 0,INFO)
diff --git a/TESTING/LIN/ctsqr01.f b/TESTING/LIN/ctsqr01.f
index a437386b..6fbeb837 100644
--- a/TESTING/LIN/ctsqr01.f
+++ b/TESTING/LIN/ctsqr01.f
@@ -109,11 +109,12 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL TESTZEROS, TS
-      INTEGER INFO, J, K, L, LWORK, LT ,MNB
+      INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
       REAL    ANORM, EPS, RESID, CNORM, DNORM
 *     ..
 *     .. Local Arrays ..
       INTEGER            ISEED( 4 )
+      COMPLEX            TQUERY( 5 ), WORKQUERY
 *     ..
 *     .. External Functions ..
       REAL     SLAMCH, CLANGE, CLANSY
@@ -145,17 +146,11 @@
       L = MAX(M,N,1)
       MNB = MAX ( MB, NB)
       LWORK = MAX(3,L)*MNB
-      IF((K.GE.MNB).OR.(MNB.GE.L))THEN
-         LT=MAX(1,L)*MNB+5
-      ELSE
-         LT=MAX(1,(L-K)/(MNB-K)+1)*L*MNB+5
-      END IF
-
 *
 *     Dynamically allocate local arrays
 *
       ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
-     $           WORK(LWORK), T(LT), C(M,N), CF(M,N),
+     $           C(M,N), CF(M,N),
      $           D(N,M), DF(N,M), LQ(L,N) )
 *
 *     Put random numbers into A and copy to AF
@@ -176,14 +171,34 @@
 *
 *     Factor the matrix A in the array AF.
 *
+      CALL CGEQR( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL CGEMQR( 'L', 'N', M, M, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMQR( 'L', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMQR( 'L', 'C', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMQR( 'R', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMQR( 'R', 'C', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'CGEQR'
-      CALL CGEQR( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL CGEQR( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *     Generate the m-by-m matrix Q
 *
       CALL CLASET( 'Full', M, M, CZERO, ONE, Q, M )
       srnamt = 'CGEMQR'
-      CALL CGEMQR( 'L', 'N', M, M, K, AF, M, T, LT, Q, M,
+      CALL CGEMQR( 'L', 'N', M, M, K, AF, M, T, TSIZE, Q, M,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -220,7 +235,7 @@
 *     Apply Q to C as Q*C
 *
       srnamt = 'CGEMQR'
-      CALL CGEMQR( 'L', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL CGEMQR( 'L', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*C - Q*C| / |C|
@@ -240,7 +255,7 @@
 *     Apply Q to C as QT*C
 *
       srnamt = 'CGEMQR'
-      CALL CGEMQR( 'L', 'C', M, N, K, AF, M, T, LT, CF, M,
+      CALL CGEMQR( 'L', 'C', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*C - QT*C| / |C|
@@ -264,7 +279,7 @@
 *     Apply Q to D as D*Q
 *
       srnamt = 'CGEMQR'
-      CALL CGEMQR( 'R', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL CGEMQR( 'R', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*Q - D*Q| / |D|
@@ -283,7 +298,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL CGEMQR( 'R', 'C', N, M, K, AF, M, T, LT, DF, N,
+      CALL CGEMQR( 'R', 'C', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*QT - D*QT| / |D|
@@ -299,15 +314,35 @@
 *     Short and wide
 *
       ELSE
+      CALL CGELQ( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL CGEMLQ( 'R', 'N', N, N, K, AF, M, TQUERY, TSIZE, Q, N,
+     $              WORKQUERY, -1, INFO )
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMLQ( 'L', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMLQ( 'L', 'C', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMLQ( 'R', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL CGEMLQ( 'R', 'C', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'CGELQ'
-      CALL CGELQ( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL CGELQ( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *
 *     Generate the n-by-n matrix Q
 *
       CALL CLASET( 'Full', N, N, CZERO, ONE, Q, N )
       srnamt = 'CGEMLQ'
-      CALL CGEMLQ( 'R', 'N', N, N, K, AF, M, T, LT, Q, N,
+      CALL CGEMLQ( 'R', 'N', N, N, K, AF, M, T, TSIZE, Q, N,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -343,7 +378,7 @@
 *
 *     Apply Q to C as Q*C
 *
-      CALL CGEMLQ( 'L', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL CGEMLQ( 'L', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*D - Q*D| / |D|
@@ -362,7 +397,7 @@
 *
 *     Apply Q to D as QT*D
 *
-      CALL CGEMLQ( 'L', 'C', N, M, K, AF, M, T, LT, DF, N,
+      CALL CGEMLQ( 'L', 'C', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*D - QT*D| / |D|
@@ -385,7 +420,7 @@
 *
 *     Apply Q to C as C*Q
 *
-      CALL CGEMLQ( 'R', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL CGEMLQ( 'R', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*Q - C*Q| / |C|
@@ -404,7 +439,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL CGEMLQ( 'R', 'C', M, N, K, AF, M, T, LT, CF, M,
+      CALL CGEMLQ( 'R', 'C', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*QT - C*QT| / |C|
diff --git a/TESTING/LIN/dchkaa.f b/TESTING/LIN/dchkaa.f
index 5d122d38..707517d2 100644
--- a/TESTING/LIN/dchkaa.f
+++ b/TESTING/LIN/dchkaa.f
@@ -907,7 +907,7 @@
             CALL DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A( 1, 1 ),
      $                   A( 1, 2 ), B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
-     $                   RWORK, RWORK( NMAX+1 ), WORK, IWORK, NOUT )
+     $                   RWORK, RWORK( NMAX+1 ), NOUT )
          ELSE
             WRITE( NOUT, FMT = 9988 )PATH
          END IF
diff --git a/TESTING/LIN/ddrvls.f b/TESTING/LIN/ddrvls.f
index d11f910e..74b0c336 100644
--- a/TESTING/LIN/ddrvls.f
+++ b/TESTING/LIN/ddrvls.f
@@ -10,7 +10,7 @@
 *
 *       SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
 *                          NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-*                          COPYB, C, S, COPYS, WORK, IWORK, NOUT )
+*                          COPYB, C, S, COPYS, NOUT )
 *
 *       .. Scalar Arguments ..
 *       LOGICAL            TSTERR
@@ -19,10 +19,10 @@
 *       ..
 *       .. Array Arguments ..
 *       LOGICAL            DOTYPE( * )
-*       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+*       INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
 *      $                   NVAL( * ), NXVAL( * )
 *       DOUBLE PRECISION   A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-*      $                   COPYS( * ), S( * ), WORK( * )
+*      $                   COPYS( * ), S( * )
 *       ..
 *
 *
@@ -169,17 +169,6 @@
 *>                      (min(MMAX,NMAX))
 *> \endverbatim
 *>
-*> \param[out] WORK
-*> \verbatim
-*>          WORK is DOUBLE PRECISION array,
-*>                      dimension (MMAX*NMAX + 4*NMAX + MMAX).
-*> \endverbatim
-*>
-*> \param[out] IWORK
-*> \verbatim
-*>          IWORK is INTEGER array, dimension (15*NMAX)
-*> \endverbatim
-*>
 *> \param[in] NOUT
 *> \verbatim
 *>          NOUT is INTEGER
@@ -201,7 +190,7 @@
 *  =====================================================================
       SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-     $                   COPYB, C, S, COPYS, WORK, IWORK, NOUT )
+     $                   COPYB, C, S, COPYS, NOUT )
 *
 *  -- LAPACK test routine (version 3.6.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -215,10 +204,10 @@
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
-      INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+      INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
      $                   NVAL( * ), NXVAL( * )
       DOUBLE PRECISION   A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-     $                   COPYS( * ), S( * ), WORK( * )
+     $                   COPYS( * ), S( * )
 *     ..
 *
 *  =====================================================================
@@ -237,12 +226,19 @@
       INTEGER            CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
      $                   ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
      $                   LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
-     $                   NFAIL, NLVL, NRHS, NROWS, NRUN, RANK, MB, LWTS
+     $                   NFAIL, NRHS, NROWS, NRUN, RANK, MB,
+     $                   MMAX, NMAX, NSMAX, LIWORK,
+     $                   LWORK_DGELS, LWORK_DGETSLS, LWORK_DGELSS,
+     $                   LWORK_DGELSY, LWORK_DGELSD
       DOUBLE PRECISION   EPS, NORMA, NORMB, RCOND
 *     ..
 *     .. Local Arrays ..
-      INTEGER            ISEED( 4 ), ISEEDY( 4 )
-      DOUBLE PRECISION   RESULT( NTESTS )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 ), IWORKQUERY
+      DOUBLE PRECISION   RESULT( NTESTS ), WORKQUERY
+*     ..
+*     .. Allocatable Arrays ..
+      DOUBLE PRECISION, ALLOCATABLE :: WORK (:)
+      INTEGER, ALLOCATABLE :: IWORK (:)
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
@@ -302,6 +298,71 @@
       CALL XLAENV( 2, 2 )
       CALL XLAENV( 9, SMLSIZ )
 *
+*     Compute maximal workspace needed for all routines
+*
+      NMAX = 0
+      MMAX = 0
+      NSMAX = 0
+      DO I = 1, NM
+         IF ( MVAL( I ).GT.MMAX ) THEN
+            MMAX = MVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NN
+         IF ( NVAL( I ).GT.NMAX ) THEN
+            NMAX = NVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NNS
+         IF ( NSVAL( I ).GT.NSMAX ) THEN
+            NSMAX = NSVAL( I )
+         END IF
+      ENDDO
+      M = MMAX
+      N = NMAX
+      NRHS = NSMAX
+      LDA = MAX( 1, M )
+      LDB = MAX( 1, M, N )
+      MNMIN = MAX( MIN( M, N ), 1 )
+*
+*     Compute workspace needed for routines
+*     DQRT14, DQRT17 (two side cases), DQRT15 and DQRT12
+*
+      LWORK = MAX( ( M+N )*NRHS,
+     $      ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ),
+     $      MAX( M+MNMIN, NRHS*MNMIN,2*N+M ),
+     $      MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) )
+*
+*     Compute workspace needed for DGELS
+      CALL DGELS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $            WORKQUERY, -1, INFO )
+      LWORK_DGELS = INT ( WORKQUERY )
+*     Compute workspace needed for DGETSLS
+      CALL DGETSLS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $              WORKQUERY, -1, INFO )
+      LWORK_DGETSLS = INT( WORKQUERY )
+*     Compute workspace needed for DGELSY
+      CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORKQUERY,
+     $             RCOND, CRANK, WORKQUERY, -1, INFO )
+      LWORK_DGELSY = INT( WORKQUERY )
+*     Compute workspace needed for DGELSS
+      CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S,
+     $             RCOND, CRANK, WORKQUERY, -1 , INFO )
+      LWORK_DGELSS = INT( WORKQUERY )
+*     Compute workspace needed for DGELSD
+      CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S,
+     $             RCOND, CRANK, WORKQUERY, -1, IWORKQUERY, INFO )
+      LWORK_DGELSD = INT( WORKQUERY )
+*     Compute LIWORK workspace needed for DGELSY and DGELSD
+      LIWORK = MAX( 1, N, IWORKQUERY )
+*     Compute LWORK workspace needed for all functions
+      LWORK = MAX( 1, LWORK, LWORK_DGELS, LWORK_DGETSLS, LWORK_DGELSY,
+     $             LWORK_DGELSS, LWORK_DGELSD )
+      LWLSY = LWORK
+*
+      ALLOCATE( WORK( LWORK ) )
+      ALLOCATE( IWORK( LIWORK ) )
+*
       DO 150 IM = 1, NM
          M = MVAL( IM )
          LDA = MAX( 1, M )
@@ -311,20 +372,9 @@
             MNMIN = MAX(MIN( M, N ),1)
             LDB = MAX( 1, M, N )
             MB = (MNMIN+1)
-            IF(MNMIN.NE.MB) THEN
-              LWTS = (((LDB-MB)/(MB-MNMIN))*MNMIN+MNMIN*2)*MB+5
-            ELSE
-              LWTS = 2*MNMIN+5
-            END IF
 *
             DO 130 INS = 1, NNS
                NRHS = NSVAL( INS )
-               NLVL = MAX( INT( LOG( MAX( ONE, DBLE( MNMIN ) ) /
-     $                DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
-               LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
-     $                 M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
-     $                 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2,LWTS)
-     $
 *
                DO 120 IRANK = 1, 2
                   DO 110 ISCALE = 1, 3
@@ -570,11 +620,6 @@
                            IWORK( J ) = 0
    70                   CONTINUE
 *
-*                       Set LWLSY to the adequate value.
-*
-                        LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
-     $                          2*MNMIN+NB*NRHS )
-*
                         CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
      $                               LDB )
@@ -768,6 +813,9 @@
  9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
      $      ', MB=', I4,', NB=', I4,', type', I2,
      $      ', test(', I2, ')=', G12.5 )
+*
+      DEALLOCATE( WORK )
+      DEALLOCATE( IWORK )
       RETURN
 *
 *     End of DDRVLS
diff --git a/TESTING/LIN/derrtsqr.f b/TESTING/LIN/derrtsqr.f
index 4a5ad5e6..a83ed1fb 100644
--- a/TESTING/LIN/derrtsqr.f
+++ b/TESTING/LIN/derrtsqr.f
@@ -161,13 +161,13 @@
       CALL DGEMQR( 'L', 'N', 2, 1, 0, A, 0, TAU, 1, C, 1, W, 1,INFO)
       CALL CHKXER( 'DGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL DGEMQR( 'R', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL DGEMQR( 'R', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'DGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL DGEMQR( 'L', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL DGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'DGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 11
-      CALL DGEMQR( 'L', 'N', 2, 1, 1, A, 1, TAU, 6, C, 0, W, 1,INFO)
+      CALL DGEMQR( 'L', 'N', 2, 1, 1, A, 2, TAU, 6, C, 0, W, 1,INFO)
       CALL CHKXER( 'DGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 13
       CALL DGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 6, C, 2, W, 0,INFO)
diff --git a/TESTING/LIN/dtsqr01.f b/TESTING/LIN/dtsqr01.f
index a9ac1635..d8f34cba 100644
--- a/TESTING/LIN/dtsqr01.f
+++ b/TESTING/LIN/dtsqr01.f
@@ -110,11 +110,12 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL TESTZEROS, TS
-      INTEGER INFO, J, K, L, LWORK, LT ,MNB
+      INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
       DOUBLE PRECISION   ANORM, EPS, RESID, CNORM, DNORM
 *     ..
 *     .. Local Arrays ..
       INTEGER            ISEED( 4 )
+      DOUBLE PRECISION   TQUERY( 5 ), WORKQUERY
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
@@ -146,17 +147,11 @@
       L = MAX(M,N,1)
       MNB = MAX ( MB, NB)
       LWORK = MAX(3,L)*MNB
-      IF((K.GE.MNB).OR.(MNB.GE.L))THEN
-         LT=MAX(1,L)*MNB+5
-      ELSE
-         LT=MAX(1,(L-K)/(MNB-K)+1)*L*MNB+5
-      END IF
-
 *
 *     Dynamically allocate local arrays
 *
       ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
-     $           WORK(LWORK), T(LT), C(M,N), CF(M,N),
+     $           C(M,N), CF(M,N),
      $           D(N,M), DF(N,M), LQ(L,N) )
 *
 *     Put random numbers into A and copy to AF
@@ -177,14 +172,34 @@
 *
 *     Factor the matrix A in the array AF.
 *
+      CALL DGEQR( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL DGEMQR( 'L', 'N', M, M, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMQR( 'L', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMQR( 'L', 'T', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMQR( 'R', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMQR( 'R', 'T', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'DGEQR'
-      CALL DGEQR( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL DGEQR( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *     Generate the m-by-m matrix Q
 *
       CALL DLASET( 'Full', M, M, ZERO, ONE, Q, M )
       srnamt = 'DGEMQR'
-      CALL DGEMQR( 'L', 'N', M, M, K, AF, M, T, LT, Q, M,
+      CALL DGEMQR( 'L', 'N', M, M, K, AF, M, T, TSIZE, Q, M,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -221,7 +236,7 @@
 *     Apply Q to C as Q*C
 *
       srnamt = 'DGEMQR'
-      CALL DGEMQR( 'L', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL DGEMQR( 'L', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*C - Q*C| / |C|
@@ -241,7 +256,7 @@
 *     Apply Q to C as QT*C
 *
       srnamt = 'DGEMQR'
-      CALL DGEMQR( 'L', 'T', M, N, K, AF, M, T, LT, CF, M,
+      CALL DGEMQR( 'L', 'T', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*C - QT*C| / |C|
@@ -265,7 +280,7 @@
 *     Apply Q to D as D*Q
 *
       srnamt = 'DGEMQR'
-      CALL DGEMQR( 'R', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL DGEMQR( 'R', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*Q - D*Q| / |D|
@@ -284,7 +299,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL DGEMQR( 'R', 'T', N, M, K, AF, M, T, LT, DF, N,
+      CALL DGEMQR( 'R', 'T', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*QT - D*QT| / |D|
@@ -300,15 +315,35 @@
 *     Short and wide
 *
       ELSE
+      CALL DGELQ( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL DGEMLQ( 'R', 'N', N, N, K, AF, M, TQUERY, TSIZE, Q, N,
+     $              WORKQUERY, -1, INFO )
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMLQ( 'L', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMLQ( 'L', 'T', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMLQ( 'R', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL DGEMLQ( 'R', 'T', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'DGELQ'
-      CALL DGELQ( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL DGELQ( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *
 *     Generate the n-by-n matrix Q
 *
       CALL DLASET( 'Full', N, N, ZERO, ONE, Q, N )
       srnamt = 'DGEMLQ'
-      CALL DGEMLQ( 'R', 'N', N, N, K, AF, M, T, LT, Q, N,
+      CALL DGEMLQ( 'R', 'N', N, N, K, AF, M, T, TSIZE, Q, N,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -344,7 +379,7 @@
 *
 *     Apply Q to C as Q*C
 *
-      CALL DGEMLQ( 'L', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL DGEMLQ( 'L', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*D - Q*D| / |D|
@@ -363,7 +398,7 @@
 *
 *     Apply Q to D as QT*D
 *
-      CALL DGEMLQ( 'L', 'T', N, M, K, AF, M, T, LT, DF, N,
+      CALL DGEMLQ( 'L', 'T', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*D - QT*D| / |D|
@@ -386,7 +421,7 @@
 *
 *     Apply Q to C as C*Q
 *
-      CALL DGEMLQ( 'R', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL DGEMLQ( 'R', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*Q - C*Q| / |C|
@@ -405,7 +440,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL DGEMLQ( 'R', 'T', M, N, K, AF, M, T, LT, CF, M,
+      CALL DGEMLQ( 'R', 'T', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*QT - C*QT| / |C|
diff --git a/TESTING/LIN/schkaa.f b/TESTING/LIN/schkaa.f
index 675e32f1..410379b3 100644
--- a/TESTING/LIN/schkaa.f
+++ b/TESTING/LIN/schkaa.f
@@ -904,7 +904,7 @@
             CALL SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A( 1, 1 ),
      $                   A( 1, 2 ), B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
-     $                   RWORK, RWORK( NMAX+1 ), WORK, IWORK, NOUT )
+     $                   RWORK, RWORK( NMAX+1 ), NOUT )
          ELSE
             WRITE( NOUT, FMT = 9988 )PATH
          END IF
diff --git a/TESTING/LIN/sdrvls.f b/TESTING/LIN/sdrvls.f
index 03598937..d6a55708 100644
--- a/TESTING/LIN/sdrvls.f
+++ b/TESTING/LIN/sdrvls.f
@@ -10,7 +10,7 @@
 *
 *       SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
 *                          NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-*                          COPYB, C, S, COPYS, WORK, IWORK, NOUT )
+*                          COPYB, C, S, COPYS, NOUT )
 *
 *       .. Scalar Arguments ..
 *       LOGICAL            TSTERR
@@ -19,10 +19,10 @@
 *       ..
 *       .. Array Arguments ..
 *       LOGICAL            DOTYPE( * )
-*       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+*       INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
 *      $                   NVAL( * ), NXVAL( * )
 *       REAL               A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-*      $                   COPYS( * ), S( * ), WORK( * )
+*      $                   COPYS( * ), S( * )
 *       ..
 *
 *
@@ -169,17 +169,6 @@
 *>                      (min(MMAX,NMAX))
 *> \endverbatim
 *>
-*> \param[out] WORK
-*> \verbatim
-*>          WORK is REAL array,
-*>                      dimension (MMAX*NMAX + 4*NMAX + MMAX).
-*> \endverbatim
-*>
-*> \param[out] IWORK
-*> \verbatim
-*>          IWORK is INTEGER array, dimension (15*NMAX)
-*> \endverbatim
-*>
 *> \param[in] NOUT
 *> \verbatim
 *>          NOUT is INTEGER
@@ -201,7 +190,7 @@
 *  =====================================================================
       SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-     $                   COPYB, C, S, COPYS, WORK, IWORK, NOUT )
+     $                   COPYB, C, S, COPYS, NOUT )
 *
 *  -- LAPACK test routine (version 3.6.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -215,10 +204,10 @@
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
-      INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+      INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
      $                   NVAL( * ), NXVAL( * )
       REAL               A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-     $                   COPYS( * ), S( * ), WORK( * )
+     $                   COPYS( * ), S( * )
 *     ..
 *
 *  =====================================================================
@@ -237,12 +226,19 @@
       INTEGER            CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
      $                   ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
      $                   LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
-     $                   NFAIL, NLVL, NRHS, NROWS, NRUN, RANK, MB, LWTS
+     $                   NFAIL, NRHS, NROWS, NRUN, RANK, MB,
+     $                   MMAX, NMAX, NSMAX, LIWORK,
+     $                   LWORK_SGELS, LWORK_SGETSLS, LWORK_SGELSS,
+     $                   LWORK_SGELSY, LWORK_SGELSD
       REAL               EPS, NORMA, NORMB, RCOND
 *     ..
 *     .. Local Arrays ..
-      INTEGER            ISEED( 4 ), ISEEDY( 4 )
-      REAL               RESULT( NTESTS )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 ), IWORKQUERY
+      REAL               RESULT( NTESTS ), WORKQUERY
+*     ..
+*     .. Allocatable Arrays ..
+      REAL, ALLOCATABLE :: WORK (:)
+      INTEGER, ALLOCATABLE :: IWORK (:)
 *     ..
 *     .. External Functions ..
       REAL               SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
@@ -302,6 +298,71 @@
       CALL XLAENV( 2, 2 )
       CALL XLAENV( 9, SMLSIZ )
 *
+*     Compute maximal workspace needed for all routines
+*
+      NMAX = 0
+      MMAX = 0
+      NSMAX = 0
+      DO I = 1, NM
+         IF ( MVAL( I ).GT.MMAX ) THEN
+            MMAX = MVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NN
+         IF ( NVAL( I ).GT.NMAX ) THEN
+            NMAX = NVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NNS
+         IF ( NSVAL( I ).GT.NSMAX ) THEN
+            NSMAX = NSVAL( I )
+         END IF
+      ENDDO
+      M = MMAX
+      N = NMAX
+      NRHS = NSMAX
+      LDA = MAX( 1, M )
+      LDB = MAX( 1, M, N )
+      MNMIN = MAX( MIN( M, N ), 1 )
+*
+*     Compute workspace needed for routines
+*     SQRT14, SQRT17 (two side cases), SQRT15 and SQRT12
+*
+      LWORK = MAX( ( M+N )*NRHS,
+     $      ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ),
+     $      MAX( M+MNMIN, NRHS*MNMIN,2*N+M ),
+     $      MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) )
+*
+*     Compute workspace needed for SGELS
+      CALL SGELS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $            WORKQUERY, -1, INFO )
+      LWORK_SGELS = INT ( WORKQUERY )
+*     Compute workspace needed for SGETSLS
+      CALL SGETSLS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $              WORKQUERY, -1, INFO )
+      LWORK_SGETSLS = INT( WORKQUERY )
+*     Compute workspace needed for SGELSY
+      CALL SGELSY( M, N, NRHS, A, LDA, B, LDB, IWORKQUERY,
+     $             RCOND, CRANK, WORKQUERY, -1, INFO )
+      LWORK_SGELSY = INT( WORKQUERY )
+*     Compute workspace needed for SGELSS
+      CALL SGELSS( M, N, NRHS, A, LDA, B, LDB, S,
+     $             RCOND, CRANK, WORKQUERY, -1 , INFO )
+      LWORK_SGELSS = INT( WORKQUERY )
+*     Compute workspace needed for SGELSD
+      CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S,
+     $             RCOND, CRANK, WORKQUERY, -1, IWORKQUERY, INFO )
+      LWORK_SGELSD = INT( WORKQUERY )
+*     Compute LIWORK workspace needed for SGELSY and SGELSD
+      LIWORK = MAX( 1, N, IWORKQUERY )
+*     Compute LWORK workspace needed for all functions
+      LWORK = MAX( 1, LWORK, LWORK_SGELS, LWORK_SGETSLS, LWORK_SGELSY,
+     $             LWORK_SGELSS, LWORK_SGELSD )
+      LWLSY = LWORK
+*
+      ALLOCATE( WORK( LWORK ) )
+      ALLOCATE( IWORK( LIWORK ) )
+*
       DO 150 IM = 1, NM
          M = MVAL( IM )
          LDA = MAX( 1, M )
@@ -311,20 +372,9 @@
             MNMIN = MAX(MIN( M, N ),1)
             LDB = MAX( 1, M, N )
             MB = (MNMIN+1)
-            IF(MNMIN.NE.MB) THEN
-              LWTS = (((LDB-MB)/(MB-MNMIN))*MNMIN+MNMIN*2)*MB+5
-            ELSE
-              LWTS = 2*MNMIN+5
-            END IF
 *
             DO 130 INS = 1, NNS
                NRHS = NSVAL( INS )
-               NLVL = MAX( INT( LOG( MAX( ONE, REAL( MNMIN ) ) /
-     $                REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
-               LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
-     $                 M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
-     $                 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2,LWTS)
-     $
 *
                DO 120 IRANK = 1, 2
                   DO 110 ISCALE = 1, 3
@@ -570,11 +620,6 @@
                            IWORK( J ) = 0
    70                   CONTINUE
 *
-*                       Set LWLSY to the adequate value.
-*
-                        LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
-     $                          2*MNMIN+NB*NRHS )
-*
                         CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
      $                               LDB )
@@ -768,6 +813,9 @@
  9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
      $      ', MB=', I4,', NB=', I4,', type', I2,
      $      ', test(', I2, ')=', G12.5 )
+*
+      DEALLOCATE( WORK )
+      DEALLOCATE( IWORK )
       RETURN
 *
 *     End of SDRVLS
diff --git a/TESTING/LIN/serrtsqr.f b/TESTING/LIN/serrtsqr.f
index eddadbee..cb7c3266 100644
--- a/TESTING/LIN/serrtsqr.f
+++ b/TESTING/LIN/serrtsqr.f
@@ -161,13 +161,13 @@
       CALL SGEMQR( 'L', 'N', 2, 1, 0, A, 0, TAU, 1, C, 1, W, 1,INFO)
       CALL CHKXER( 'SGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL SGEMQR( 'R', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL SGEMQR( 'R', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'SGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL SGEMQR( 'L', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL SGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'SGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 11
-      CALL SGEMQR( 'L', 'N', 2, 1, 1, A, 1, TAU, 6, C, 0, W, 1,INFO)
+      CALL SGEMQR( 'L', 'N', 2, 1, 1, A, 2, TAU, 6, C, 0, W, 1,INFO)
       CALL CHKXER( 'SGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 13
       CALL SGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 6, C, 2, W, 0,INFO)
diff --git a/TESTING/LIN/stsqr01.f b/TESTING/LIN/stsqr01.f
index 4cebfc88..3e4e3d09 100644
--- a/TESTING/LIN/stsqr01.f
+++ b/TESTING/LIN/stsqr01.f
@@ -110,11 +110,12 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL TESTZEROS, TS
-      INTEGER INFO, J, K, L, LWORK, LT ,MNB
+      INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
       REAL   ANORM, EPS, RESID, CNORM, DNORM
 *     ..
 *     .. Local Arrays ..
       INTEGER            ISEED( 4 )
+      REAL               TQUERY( 5 ), WORKQUERY
 *     ..
 *     .. External Functions ..
       REAL     SLAMCH, SLANGE, SLANSY
@@ -146,17 +147,11 @@
       L = MAX(M,N,1)
       MNB = MAX ( MB, NB)
       LWORK = MAX(3,L)*MNB
-      IF((K.GE.MNB).OR.(MNB.GE.L))THEN
-         LT=MAX(1,L)*MNB+5
-      ELSE
-         LT=MAX(1,(L-K)/(MNB-K)+1)*L*MNB+5
-      END IF
-
 *
 *     Dynamically allocate local arrays
 *
       ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
-     $           WORK(LWORK), T(LT), C(M,N), CF(M,N),
+     $           C(M,N), CF(M,N),
      $           D(N,M), DF(N,M), LQ(L,N) )
 *
 *     Put random numbers into A and copy to AF
@@ -177,14 +172,34 @@
 *
 *     Factor the matrix A in the array AF.
 *
+      CALL SGEQR( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL SGEMQR( 'L', 'N', M, M, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMQR( 'L', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMQR( 'L', 'T', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMQR( 'R', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMQR( 'R', 'T', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'SGEQR'
-      CALL SGEQR( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL SGEQR( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *     Generate the m-by-m matrix Q
 *
       CALL SLASET( 'Full', M, M, ZERO, ONE, Q, M )
       srnamt = 'SGEMQR'
-      CALL SGEMQR( 'L', 'N', M, M, K, AF, M, T, LT, Q, M,
+      CALL SGEMQR( 'L', 'N', M, M, K, AF, M, T, TSIZE, Q, M,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -221,7 +236,7 @@
 *     Apply Q to C as Q*C
 *
       srnamt = 'DGEQR'
-      CALL SGEMQR( 'L', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL SGEMQR( 'L', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*C - Q*C| / |C|
@@ -241,7 +256,7 @@
 *     Apply Q to C as QT*C
 *
       srnamt = 'DGEQR'
-      CALL SGEMQR( 'L', 'T', M, N, K, AF, M, T, LT, CF, M,
+      CALL SGEMQR( 'L', 'T', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*C - QT*C| / |C|
@@ -265,7 +280,7 @@
 *     Apply Q to D as D*Q
 *
       srnamt = 'DGEQR'
-      CALL SGEMQR( 'R', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL SGEMQR( 'R', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*Q - D*Q| / |D|
@@ -284,7 +299,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL SGEMQR( 'R', 'T', N, M, K, AF, M, T, LT, DF, N,
+      CALL SGEMQR( 'R', 'T', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*QT - D*QT| / |D|
@@ -300,15 +315,35 @@
 *     Short and wide
 *
       ELSE
+      CALL SGELQ( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL SGEMLQ( 'R', 'N', N, N, K, AF, M, TQUERY, TSIZE, Q, N,
+     $              WORKQUERY, -1, INFO )
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMLQ( 'L', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMLQ( 'L', 'T', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMLQ( 'R', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL SGEMLQ( 'R', 'T', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'SGELQ'
-      CALL SGELQ( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL SGELQ( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *
 *     Generate the n-by-n matrix Q
 *
       CALL SLASET( 'Full', N, N, ZERO, ONE, Q, N )
-      srnamt = 'SGEMQR'
-      CALL SGEMLQ( 'R', 'N', N, N, K, AF, M, T, LT, Q, N,
+      srnamt = 'SGEMLQ'
+      CALL SGEMLQ( 'R', 'N', N, N, K, AF, M, T, TSIZE, Q, N,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -344,7 +379,7 @@
 *
 *     Apply Q to C as Q*C
 *
-      CALL SGEMLQ( 'L', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL SGEMLQ( 'L', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*D - Q*D| / |D|
@@ -363,7 +398,7 @@
 *
 *     Apply Q to D as QT*D
 *
-      CALL SGEMLQ( 'L', 'T', N, M, K, AF, M, T, LT, DF, N,
+      CALL SGEMLQ( 'L', 'T', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*D - QT*D| / |D|
@@ -386,7 +421,7 @@
 *
 *     Apply Q to C as C*Q
 *
-      CALL SGEMLQ( 'R', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL SGEMLQ( 'R', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*Q - C*Q| / |C|
@@ -405,7 +440,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL SGEMLQ( 'R', 'T', M, N, K, AF, M, T, LT, CF, M,
+      CALL SGEMLQ( 'R', 'T', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*QT - C*QT| / |C|
diff --git a/TESTING/LIN/zchkaa.f b/TESTING/LIN/zchkaa.f
index 818f1e63..5a41ab32 100644
--- a/TESTING/LIN/zchkaa.f
+++ b/TESTING/LIN/zchkaa.f
@@ -1049,8 +1049,7 @@
             CALL ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A( 1, 1 ),
      $                   A( 1, 2 ), A( 1, 3 ), A( 1, 4 ), A( 1, 5 ),
-     $                   S( 1 ), S( NMAX+1 ), WORK, RWORK, IWORK,
-     $                   NOUT )
+     $                   S( 1 ), S( NMAX+1 ), NOUT )
          ELSE
             WRITE( NOUT, FMT = 9989 )PATH
          END IF
diff --git a/TESTING/LIN/zdrvls.f b/TESTING/LIN/zdrvls.f
index fe63b540..13a9263e 100644
--- a/TESTING/LIN/zdrvls.f
+++ b/TESTING/LIN/zdrvls.f
@@ -10,7 +10,7 @@
 *
 *       SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
 *                          NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-*                          COPYB, C, S, COPYS, WORK, RWORK, IWORK, NOUT )
+*                          COPYB, C, S, COPYS, NOUT )
 *
 *       .. Scalar Arguments ..
 *       LOGICAL            TSTERR
@@ -19,11 +19,10 @@
 *       ..
 *       .. Array Arguments ..
 *       LOGICAL            DOTYPE( * )
-*       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+*       INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
 *      $                   NVAL( * ), NXVAL( * )
-*       DOUBLE PRECISION   COPYS( * ), RWORK( * ), S( * )
-*       COMPLEX*16         A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-*      $                   WORK( * )
+*       DOUBLE PRECISION   COPYS( * ), S( * )
+*       COMPLEX*16         A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
 *       ..
 *
 *
@@ -32,8 +31,8 @@
 *>
 *> \verbatim
 *>
-*> ZDRVLS tests the least squares driver routines ZGELS, CGELSS, ZGELSY
-*> and CGELSD.
+*> ZDRVLS tests the least squares driver routines ZGELS, ZGETSLS, ZGELSS, ZGELSY
+*> and ZGELSD.
 *> \endverbatim
 *
 *  Arguments:
@@ -170,22 +169,6 @@
 *>                      (min(MMAX,NMAX))
 *> \endverbatim
 *>
-*> \param[out] WORK
-*> \verbatim
-*>          WORK is COMPLEX*16 array, dimension
-*>                      (MMAX*NMAX + 4*NMAX + MMAX).
-*> \endverbatim
-*>
-*> \param[out] RWORK
-*> \verbatim
-*>          RWORK is DOUBLE PRECISION array, dimension (5*NMAX-1)
-*> \endverbatim
-*>
-*> \param[out] IWORK
-*> \verbatim
-*>          IWORK is INTEGER array, dimension (15*NMAX)
-*> \endverbatim
-*>
 *> \param[in] NOUT
 *> \verbatim
 *>          NOUT is INTEGER
@@ -207,7 +190,7 @@
 *  =====================================================================
       SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
      $                   NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
-     $                   COPYB, C, S, COPYS, WORK, RWORK, IWORK, NOUT )
+     $                   COPYB, C, S, COPYS, NOUT )
 *
 *  -- LAPACK test routine (version 3.6.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
@@ -221,11 +204,10 @@
 *     ..
 *     .. Array Arguments ..
       LOGICAL            DOTYPE( * )
-      INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
+      INTEGER            MVAL( * ), NBVAL( * ), NSVAL( * ),
      $                   NVAL( * ), NXVAL( * )
-      DOUBLE PRECISION   COPYS( * ), RWORK( * ), S( * )
-      COMPLEX*16         A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
-     $                   WORK( * )
+      DOUBLE PRECISION   COPYS( * ), S( * )
+      COMPLEX*16         A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
 *     ..
 *
 *  =====================================================================
@@ -247,12 +229,22 @@
       INTEGER            CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
      $                   ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
      $                   LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
-     $                   NFAIL, NRHS, NROWS, NRUN, RANK, MB, LWTS
+     $                   NFAIL, NRHS, NROWS, NRUN, RANK, MB,
+     $                   MMAX, NMAX, NSMAX, LIWORK, LRWORK,
+     $                   LWORK_ZGELS, LWORK_ZGETSLS, LWORK_ZGELSS,
+     $                   LWORK_ZGELSY, LWORK_ZGELSD,
+     $                   LRWORK_ZGELSY, LRWORK_ZGELSS, LRWORK_ZGELSD
       DOUBLE PRECISION   EPS, NORMA, NORMB, RCOND
 *     ..
 *     .. Local Arrays ..
-      INTEGER            ISEED( 4 ), ISEEDY( 4 )
-      DOUBLE PRECISION   RESULT( NTESTS )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 ), IWORKQUERY
+      DOUBLE PRECISION   RESULT( NTESTS ), RWORKQUERY
+      COMPLEX*16         WORKQUERY
+*     ..
+*     .. Allocatable Arrays ..
+      COMPLEX*16, ALLOCATABLE :: WORK (:)
+      DOUBLE PRECISION, ALLOCATABLE :: RWORK (:)
+      INTEGER, ALLOCATABLE :: IWORK (:)
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
@@ -265,7 +257,7 @@
      $                   ZQRT16, ZGETSLS
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          DBLE, MAX, MIN, SQRT
+      INTRINSIC          DBLE, MAX, MIN, INT, SQRT
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -309,6 +301,77 @@
      $   CALL ALAHD( NOUT, PATH )
       INFOT = 0
 *
+*     Compute maximal workspace needed for all routines
+*
+      NMAX = 0
+      MMAX = 0
+      NSMAX = 0
+      DO I = 1, NM
+         IF ( MVAL( I ).GT.MMAX ) THEN
+            MMAX = MVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NN
+         IF ( NVAL( I ).GT.NMAX ) THEN
+            NMAX = NVAL( I )
+         END IF
+      ENDDO
+      DO I = 1, NNS
+         IF ( NSVAL( I ).GT.NSMAX ) THEN
+            NSMAX = NSVAL( I )
+         END IF
+      ENDDO
+      M = MMAX
+      N = NMAX
+      NRHS = NSMAX
+      LDA = MAX( 1, M )
+      LDB = MAX( 1, M, N )
+      MNMIN = MAX( MIN( M, N ), 1 )
+*
+*     Compute workspace needed for routines
+*     ZQRT14, ZQRT17 (two side cases), ZQRT15 and ZQRT12
+*
+      LWORK = MAX( ( M+N )*NRHS,
+     $      ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ),
+     $      MAX( M+MNMIN, NRHS*MNMIN,2*N+M ),
+     $      MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) )
+*
+*     Compute workspace needed for ZGELS
+      CALL ZGELS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $            WORKQUERY, -1, INFO )
+      LWORK_ZGELS = INT ( WORKQUERY )
+*     Compute workspace needed for ZGETSLS
+      CALL ZGETSLS( 'N', M, N, NRHS, A, LDA, B, LDB,
+     $              WORKQUERY, -1, INFO )
+      LWORK_ZGETSLS = INT( WORKQUERY )
+*     Compute workspace needed for ZGELSY
+      CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWORKQUERY,
+     $             RCOND, CRANK, WORKQUERY, -1, RWORK, INFO )
+      LWORK_ZGELSY = INT( WORKQUERY )
+      LRWORK_ZGELSY = 2*N
+*     Compute workspace needed for ZGELSS
+      CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
+     $             RCOND, CRANK, WORKQUERY, -1 , RWORK, INFO )
+      LWORK_ZGELSS = INT( WORKQUERY )
+      LRWORK_ZGELSS = 5*MNMIN
+*     Compute workspace needed for ZGELSD
+      CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, CRANK,
+     $             WORKQUERY, -1, RWORKQUERY, IWORKQUERY, INFO )
+      LWORK_ZGELSD = INT( WORKQUERY )
+      LRWORK_ZGELSD = INT( RWORKQUERY )
+*     Compute LIWORK workspace needed for ZGELSY and ZGELSD
+      LIWORK = MAX( 1, N, IWORKQUERY )
+*     Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
+      LRWORK = MAX( 1, LRWORK_ZGELSY, LRWORK_ZGELSS, LRWORK_ZGELSD )
+*     Compute LWORK workspace needed for all functions
+      LWORK = MAX( 1, LWORK, LWORK_ZGELS, LWORK_ZGETSLS, LWORK_ZGELSY,
+     $             LWORK_ZGELSS, LWORK_ZGELSD )
+      LWLSY = LWORK
+*
+      ALLOCATE( WORK( LWORK ) )
+      ALLOCATE( IWORK( LIWORK ) )
+      ALLOCATE( RWORK( LRWORK ) )
+*
       DO 140 IM = 1, NM
          M = MVAL( IM )
          LDA = MAX( 1, M )
@@ -318,16 +381,9 @@
             MNMIN = MAX(MIN( M, N ),1)
             LDB = MAX( 1, M, N )
             MB = (MNMIN+1)
-            IF(MNMIN.NE.MB) THEN
-              LWTS = (((LDB-MB)/(MB-MNMIN))*MNMIN+LDB*2)*MB+5
-            ELSE
-              LWTS = 2*MNMIN+5
-            END IF
 *
             DO 120 INS = 1, NNS
                NRHS = NSVAL( INS )
-               LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
-     $                 M*N+4*MNMIN+MAX( M, N ), 2*N+M, LWTS )
 *
                DO 110 IRANK = 1, 2
                   DO 100 ISCALE = 1, 3
@@ -578,12 +634,6 @@
                            IWORK( J ) = 0
    70                   CONTINUE
 *
-*                       Set LWLSY to the adequate value.
-*
-                        LWLSY = MNMIN + MAX( 2*MNMIN, NB*( N+1 ),
-     $                          MNMIN+NB*NRHS )
-                        LWLSY = MAX( 1, LWLSY )
-*
                         SRNAMT = 'ZGELSY'
                         CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
      $                               RCOND, CRANK, WORK, LWLSY, RWORK,
@@ -774,6 +824,10 @@
  9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
      $      ', MB=', I4,', NB=', I4,', type', I2,
      $      ', test(', I2, ')=', G12.5 )
+*
+      DEALLOCATE( WORK )
+      DEALLOCATE( IWORK )
+      DEALLOCATE( RWORK )
       RETURN
 *
 *     End of ZDRVLS
diff --git a/TESTING/LIN/zerrtsqr.f b/TESTING/LIN/zerrtsqr.f
index 3aa3e4a5..5550e219 100644
--- a/TESTING/LIN/zerrtsqr.f
+++ b/TESTING/LIN/zerrtsqr.f
@@ -161,13 +161,13 @@
       CALL ZGEMQR( 'L', 'N', 2, 1, 0, A, 0, TAU, 1, C, 1, W, 1,INFO)
       CALL CHKXER( 'ZGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL ZGEMQR( 'R', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL ZGEMQR( 'R', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'ZGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 9
-      CALL ZGEMQR( 'L', 'N', 2, 2, 1, A, 1, TAU, 0, C, 1, W, 1,INFO)
+      CALL ZGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 0, C, 1, W, 1,INFO)
       CALL CHKXER( 'ZGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 11
-      CALL ZGEMQR( 'L', 'N', 2, 1, 1, A, 1, TAU, 6, C, 0, W, 1,INFO)
+      CALL ZGEMQR( 'L', 'N', 2, 1, 1, A, 2, TAU, 6, C, 0, W, 1,INFO)
       CALL CHKXER( 'ZGEMQR', INFOT, NOUT, LERR, OK )
       INFOT = 13
       CALL ZGEMQR( 'L', 'N', 2, 2, 1, A, 2, TAU, 6, C, 2, W, 0,INFO)
diff --git a/TESTING/LIN/ztsqr01.f b/TESTING/LIN/ztsqr01.f
index 38ace9c8..acc027e0 100644
--- a/TESTING/LIN/ztsqr01.f
+++ b/TESTING/LIN/ztsqr01.f
@@ -109,11 +109,12 @@
 *     ..
 *     .. Local Scalars ..
       LOGICAL TESTZEROS, TS
-      INTEGER INFO, J, K, L, LWORK, LT ,MNB
+      INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
       DOUBLE PRECISION   ANORM, EPS, RESID, CNORM, DNORM
 *     ..
 *     .. Local Arrays ..
       INTEGER            ISEED( 4 )
+      COMPLEX*16         TQUERY( 5 ), WORKQUERY
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
@@ -145,17 +146,11 @@
       L = MAX(M,N,1)
       MNB = MAX ( MB, NB)
       LWORK = MAX(3,L)*MNB
-      IF((K.GE.MNB).OR.(MNB.GE.L))THEN
-         LT=MAX(1,L)*MNB+5
-      ELSE
-         LT=MAX(1,(L-K)/(MNB-K)+1)*L*MNB+5
-      END IF
-
 *
 *     Dynamically allocate local arrays
 *
       ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L),
-     $           WORK(LWORK), T(LT), C(M,N), CF(M,N),
+     $           C(M,N), CF(M,N),
      $           D(N,M), DF(N,M), LQ(L,N) )
 *
 *     Put random numbers into A and copy to AF
@@ -176,14 +171,34 @@
 *
 *     Factor the matrix A in the array AF.
 *
+      CALL ZGEQR( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL ZGEMQR( 'L', 'N', M, M, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMQR( 'L', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMQR( 'L', 'C', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMQR( 'R', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMQR( 'R', 'C', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'ZGEQR'
-      CALL ZGEQR( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL ZGEQR( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *     Generate the m-by-m matrix Q
 *
       CALL ZLASET( 'Full', M, M, CZERO, ONE, Q, M )
       srnamt = 'ZGEMQR'
-      CALL ZGEMQR( 'L', 'N', M, M, K, AF, M, T, LT, Q, M,
+      CALL ZGEMQR( 'L', 'N', M, M, K, AF, M, T, TSIZE, Q, M,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -220,7 +235,7 @@
 *     Apply Q to C as Q*C
 *
       srnamt = 'ZGEMQR'
-      CALL ZGEMQR( 'L', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL ZGEMQR( 'L', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*C - Q*C| / |C|
@@ -240,7 +255,7 @@
 *     Apply Q to C as QT*C
 *
       srnamt = 'ZGEMQR'
-      CALL ZGEMQR( 'L', 'C', M, N, K, AF, M, T, LT, CF, M,
+      CALL ZGEMQR( 'L', 'C', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*C - QT*C| / |C|
@@ -264,7 +279,7 @@
 *     Apply Q to D as D*Q
 *
       srnamt = 'ZGEMQR'
-      CALL ZGEMQR( 'R', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL ZGEMQR( 'R', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*Q - D*Q| / |D|
@@ -283,7 +298,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL ZGEMQR( 'R', 'C', N, M, K, AF, M, T, LT, DF, N,
+      CALL ZGEMQR( 'R', 'C', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |D*QT - D*QT| / |D|
@@ -299,15 +314,35 @@
 *     Short and wide
 *
       ELSE
+      CALL ZGELQ( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
+      TSIZE = INT( TQUERY( 1 ) )
+      LWORK = INT( WORKQUERY )
+      CALL ZGEMLQ( 'R', 'N', N, N, K, AF, M, TQUERY, TSIZE, Q, N,
+     $              WORKQUERY, -1, INFO )
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMLQ( 'L', 'N', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMLQ( 'L', 'C', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMLQ( 'R', 'N', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      CALL ZGEMLQ( 'R', 'C', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
+     $             WORKQUERY, -1, INFO)
+      LWORK = MAX( LWORK, INT( WORKQUERY ) )
+      ALLOCATE ( T( TSIZE ) )
+      ALLOCATE ( WORK( LWORK ) )
       srnamt = 'ZGELQ'
-      CALL ZGELQ( M, N, AF, M, T, LT, WORK, LWORK, INFO )
+      CALL ZGELQ( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
 *
 *
 *     Generate the n-by-n matrix Q
 *
       CALL ZLASET( 'Full', N, N, CZERO, ONE, Q, N )
       srnamt = 'ZGEMLQ'
-      CALL ZGEMLQ( 'R', 'N', N, N, K, AF, M, T, LT, Q, N,
+      CALL ZGEMLQ( 'R', 'N', N, N, K, AF, M, T, TSIZE, Q, N,
      $              WORK, LWORK, INFO )
 *
 *     Copy R
@@ -343,7 +378,7 @@
 *
 *     Apply Q to C as Q*C
 *
-      CALL ZGEMLQ( 'L', 'N', N, M, K, AF, M, T, LT, DF, N,
+      CALL ZGEMLQ( 'L', 'N', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |Q*D - Q*D| / |D|
@@ -362,7 +397,7 @@
 *
 *     Apply Q to D as QT*D
 *
-      CALL ZGEMLQ( 'L', 'C', N, M, K, AF, M, T, LT, DF, N,
+      CALL ZGEMLQ( 'L', 'C', N, M, K, AF, M, T, TSIZE, DF, N,
      $             WORK, LWORK, INFO)
 *
 *     Compute |QT*D - QT*D| / |D|
@@ -385,7 +420,7 @@
 *
 *     Apply Q to C as C*Q
 *
-      CALL ZGEMLQ( 'R', 'N', M, N, K, AF, M, T, LT, CF, M,
+      CALL ZGEMLQ( 'R', 'N', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*Q - C*Q| / |C|
@@ -404,7 +439,7 @@
 *
 *     Apply Q to D as D*QT
 *
-      CALL ZGEMLQ( 'R', 'C', M, N, K, AF, M, T, LT, CF, M,
+      CALL ZGEMLQ( 'R', 'C', M, N, K, AF, M, T, TSIZE, CF, M,
      $             WORK, LWORK, INFO)
 *
 *     Compute |C*QT - C*QT| / |C|