Fix ?GELQ and ?GEMLQ

8 files changed, 457 insertions, 406 deletions

diff --git a/SRC/cgelq.f b/SRC/cgelq.f
index 0abd2d72..497851f5 100644
--- a/SRC/cgelq.f
+++ b/SRC/cgelq.f
@@ -3,13 +3,13 @@
 *  ===========
 *
 *       SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX  A( LDA, * ), T( * ), WORK( * )
+*       COMPLEX           A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -120,7 +120,7 @@
 *>
 *> 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
+*> (trapezoidal) L has to be stored in the lower 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.
 *>
@@ -146,74 +146,73 @@
 *>
 *>          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
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           CLASWLQ or CGELQT
 *>
 *>  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
+*>  block sizes MB and NB returned by ILAENV, CGELQ will use either
+*>  CLASWLQ (if the matrix is short-and-wide) or CGELQT to compute
 *>  the LQ factorization.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-     $   INFO)
+     $                  INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
+      INTEGER            INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      COMPLEX  A( LDA, * ), T( * ), WORK( * )
+      COMPLEX            A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS, MINT, MINW
-      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
+      LOGICAL            LQUERY, LMINWS, MINT, MINW
+      INTEGER            MB, NB, MINTSZ, NBLCKS
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
-*     .. EXTERNAL SUBROUTINES ..
+*     ..
+*     .. External Subroutines ..
       EXTERNAL           CGELQT, CLASWLQ, XERBLA
-*     .. INTRINSIC FUNCTIONS ..
+*     ..
+*     .. Intrinsic Functions ..
       INTRINSIC          MAX, MIN, MOD
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       INTEGER            ILAENV
       EXTERNAL           ILAENV
 *     ..
-*     .. EXECUTABLE STATEMENTS ..
+*     .. Executable Statements ..
 *
-*     TEST THE INPUT ARGUMENTS
+*     Test the input arguments
 *
       INFO = 0
 *
-      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+      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
+      IF( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) THEN
+        IF( TSIZE.NE.-1 ) MINT = .TRUE.
+        IF( LWORK.NE.-1 ) MINW = .TRUE.
+      END IF
 *
 *     Determine the block size
 *
-      IF ( MIN(M,N).GT.0 ) THEN
-        MB = ILAENV( 1, 'CGELQ ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'CGELQ ', ' ', M, N, 2, -1)
+      IF( MIN( M, N ).GT.0 ) THEN
+        MB = ILAENV( 1, 'CGELQ ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'CGELQ ', ' ', M, N, 2, -1 )
       ELSE
         MB = 1
         NB = N
@@ -221,7 +220,7 @@
       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( NB.GT.M .AND. N.GT.M ) THEN
         IF( MOD( N - M, NB - M ).EQ.0 ) THEN
           NBLCKS = ( N - M ) / ( NB - M )
         ELSE
@@ -235,16 +234,16 @@
 *
       LMINWS = .FALSE.
       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
-            NB = N
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.MINTSZ )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
+          LMINWS = .TRUE.
+          MB = 1
+          NB = N
         END IF
-        IF ( LWORK.LT.MB*M ) THEN
-            LMINWS = .TRUE.
-            MB = 1
+        IF( LWORK.LT.MB*M ) THEN
+          LMINWS = .TRUE.
+          MB = 1
         END IF
       END IF
 *
@@ -262,42 +261,44 @@
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        IF ( MINT ) THEN
-          T(1) = MINTSZ
+      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 )
+          T( 1 ) = MB*M*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, MB*M )
-        ENDIF
+          WORK( 1 ) = MAX( 1, MB*M )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGELQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N).EQ.0 ) THEN
+      IF( MIN( M, N ).EQ.0 ) THEN
         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, T(4), MB, WORK, INFO)
+        CALL CGELQT( M, N, MB, A, LDA, T( 6 ), MB, WORK, INFO )
       ELSE
-        CALL CLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
-     $                    LWORK, INFO)
+        CALL CLASWLQ( M, N, MB, NB, A, LDA, T( 6 ), MB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, MB*M )
+*
+      WORK( 1 ) = MAX( 1, MB*M )
+*
       RETURN
 *
 *     End of CGELQ
diff --git a/SRC/cgemlq.f b/SRC/cgemlq.f
index 03dae76d..59df3ddf 100644
--- a/SRC/cgemlq.f
+++ b/SRC/cgemlq.f
@@ -3,15 +3,17 @@
 *  ===========
 *
 *      SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
-*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
+*     $                   TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
-*      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
+*     CHARACTER         SIDE, TRANS
+*     INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      COMPLEX           A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*     COMPLEX           A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*     ..
+*
 *> \par Purpose:
 *  =============
 *>
@@ -19,27 +21,32 @@
 *>
 *>     CGEMLQ overwrites the general real M-by-N matrix C with
 *>
-*>
 *>                      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:
 *  ==========
 *
 *> \param[in] SIDE
+*> \verbatim
 *>          SIDE is CHARACTER*1
 *>          = 'L': apply Q or Q**T from the Left;
 *>          = 'R': apply Q or Q**T from the Right.
+*> \endverbatim
 *>
 *> \param[in] TRANS
+*> \verbatim
 *>          TRANS is CHARACTER*1
 *>          = 'N':  No transpose, apply Q;
 *>          = 'T':  Transpose, apply Q**T.
+*> \endverbatim
+*>
 *> \param[in] M
 *> \verbatim
 *>          M is INTEGER
@@ -49,7 +56,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. N >= M.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,28 +64,28 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          M >= K >= 0;
-*>
+*>          If SIDE = 'L', M >= K >= 0;
+*>          if SIDE = 'R', N >= K >= 0.
 *> \endverbatim
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is COMPLEX array, dimension (LDA,K)
-*>          Part of the data structure to represent Q as returned by ZGELQ.
+*>          A is COMPLEX array, dimension
+*>                               (LDA,M) if SIDE = 'L',
+*>                               (LDA,N) if SIDE = 'R'
+*>          Part of the data structure to represent Q as returned by CGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.
-*>          If SIDE = 'L', LDA >= max(1,M);
-*>          if SIDE = 'R', LDA >= max(1,N).
+*>          The leading dimension of the array A. LDA >= max(1,K).
 *> \endverbatim
 *>
 *> \param[in] T
 *> \verbatim
 *>          T is COMPLEX array, dimension (MAX(5,TSIZE)).
-*>          Part of the data structure to represent Q as returned by ZGELQ.
+*>          Part of the data structure to represent Q as returned by CGELQ.
 *> \endverbatim
 *>
 *> \param[in] TSIZE
@@ -88,19 +95,23 @@
 *> \endverbatim
 *>
 *> \param[in,out] C
+*> \verbatim
 *>          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.
+*> \endverbatim
 *>
 *> \param[in] LDC
+*> \verbatim
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
 *>
 *> \param[out] WORK
 *> \verbatim
 *>         (workspace) COMPLEX array, dimension (MAX(1,LWORK))
-*>
 *> \endverbatim
+*>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
@@ -140,47 +151,49 @@
 *>
 *>          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
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           CLASWQR or CGELQT
 *>
 *>  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 wide-and-short) or GELQT to compute
+*>  block sizes MB and NB returned by ILAENV, CGELQ will use either
+*>  CLASWLQ (if the matrix is wide-and-short) or CGELQT to compute
 *>  the LQ factorization.
-*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  This version of CGEMLQ will use either CLAMSWLQ or CGEMLQT to 
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LAMSWLQ or GEMLQT.
+*>  Further Details in CLAMSWLQ or CGEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
-     $        C, LDC, WORK, LWORK, INFO )
+     $                   C, LDC, WORK, LWORK, INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
+      CHARACTER          SIDE, TRANS
+      INTEGER            INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      COMPLEX           A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+      COMPLEX            A( LDA, * ), T( * ), C( LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LEFT, RIGHT, TRAN, NOTRAN, LQUERY
-      INTEGER    I, II, KK, MB, NB, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
-      EXTERNAL           ZLAMSWLQ, ZGEMLQT, XERBLA
+      EXTERNAL           CLAMSWLQ, CGEMLQT, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,15 +201,15 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK.LT.0
+      LQUERY  = LWORK.EQ.-1
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'C' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(T(2))
-      NB = INT(T(3))
-      IF ( LEFT ) THEN
+      MB = INT( T( 2 ) )
+      NB = INT( T( 3 ) )
+      IF( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
@@ -204,7 +217,7 @@
         MN = N
       END IF
 *
-      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+      IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
         IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
           NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
@@ -216,34 +229,33 @@
 *
       INFO = 0
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
-         INFO = -1
+        INFO = -1
       ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
-         INFO = -2
+        INFO = -2
       ELSE IF( M.LT.0 ) THEN
         INFO = -3
-      ELSE IF( N.LT.0) THEN
+      ELSE IF( N.LT.0 ) THEN
         INFO = -4
-      ELSE IF( K.LT.0 ) THEN
+      ELSE IF( K.LT.0 .OR. K.GT.MN ) THEN
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
-     $         .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
         INFO = -11
-      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        WORK(1) = REAL( 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
 *
@@ -256,13 +268,13 @@
       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,
-     $        T(4), MB, C, LDC, WORK, INFO)
+     $                T( 6 ), MB, C, LDC, WORK, INFO )
       ELSE
-        CALL CLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $    MB, C, LDC, WORK, LWORK, INFO )
+        CALL CLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = REAL ( LW )
+      WORK( 1 ) = REAL( LW )
 *
       RETURN
 *
diff --git a/SRC/dgelq.f b/SRC/dgelq.f
index 59d9fa91..a9af9006 100644
--- a/SRC/dgelq.f
+++ b/SRC/dgelq.f
@@ -3,7 +3,7 @@
 *  ===========
 *
 *       SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
@@ -120,7 +120,7 @@
 *>
 *> 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
+*> (trapezoidal) L has to be stored in the lower 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.
 *>
@@ -146,72 +146,71 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
+*>          T(6: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
+*>  block sizes MB and NB returned by ILAENV, DGELQ will use either
+*>  DLASWLQ (if the matrix is short-and-wide) or DGELQT to compute
 *>  the LQ factorization.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-     $   INFO)
+     $                  INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
+      INTEGER            INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION  A( LDA, * ), T( * ), WORK( * )
+      DOUBLE PRECISION   A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS, MINT, MINW
-      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
+      LOGICAL            LQUERY, LMINWS, MINT, MINW
+      INTEGER            MB, NB, MINTSZ, NBLCKS
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
-*     .. EXTERNAL SUBROUTINES ..
+*     ..
+*     .. External Subroutines ..
       EXTERNAL           DGELQT, DLASWLQ, XERBLA
-*     .. INTRINSIC FUNCTIONS ..
+*     ..
+*     .. Intrinsic Functions ..
       INTRINSIC          MAX, MIN, MOD
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       INTEGER            ILAENV
       EXTERNAL           ILAENV
 *     ..
-*     .. EXECUTABLE STATEMENTS ..
+*     .. Executable Statements ..
 *
-*     TEST THE INPUT ARGUMENTS
+*     Test the input arguments
 *
       INFO = 0
 *
-      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+      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
+      IF( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) THEN
+        IF( TSIZE.NE.-1 ) MINT = .TRUE.
+        IF( LWORK.NE.-1 ) MINW = .TRUE.
+      END IF
 *
 *     Determine the block size
 *
-      IF ( MIN(M,N).GT.0 ) THEN
+      IF( MIN( M, N ).GT.0 ) THEN
         MB = ILAENV( 1, 'DGELQ ', ' ', M, N, 1, -1 )
         NB = ILAENV( 1, 'DGELQ ', ' ', M, N, 2, -1 )
       ELSE
@@ -235,14 +234,14 @@
 *
       LMINWS = .FALSE.
       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
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.MINTSZ )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             MB = 1
             NB = N
         END IF
-        IF ( LWORK.LT.MB*M ) THEN
+        IF( LWORK.LT.MB*M ) THEN
             LMINWS = .TRUE.
             MB = 1
         END IF
@@ -262,42 +261,44 @@
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        IF ( MINT ) THEN
-          T(1) = MINTSZ
+      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 )
+          T( 1 ) = MB*M*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, MB*M )
-        ENDIF
+          WORK( 1 ) = MAX( 1, MB*M )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGELQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N).EQ.0 ) THEN
+      IF( MIN( M, N ).EQ.0 ) THEN
         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, T(4), MB, WORK, INFO)
+        CALL DGELQT( M, N, MB, A, LDA, T( 6 ), MB, WORK, INFO )
       ELSE
-        CALL DLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
-     $                    LWORK, INFO)
+        CALL DLASWLQ( M, N, MB, NB, A, LDA, T( 6 ), MB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, MB*M )
+*
+      WORK( 1 ) = MAX( 1, MB*M )
+*
       RETURN
 *
 *     End of DGELQ
diff --git a/SRC/dgemlq.f b/SRC/dgemlq.f
index 17c4de5c..203ca7ec 100644
--- a/SRC/dgemlq.f
+++ b/SRC/dgemlq.f
@@ -3,15 +3,17 @@
 *  ===========
 *
 *      SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
-*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
+*     $                   TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
-*      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
+*      CHARACTER          SIDE, TRANS
+*      INTEGER            INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      DOUBLE PRECISION  A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*      DOUBLE PRECISION   A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*     ..
+*
 *> \par Purpose:
 *  =============
 *>
@@ -19,27 +21,32 @@
 *>
 *>     DGEMLQ 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 real orthogonal matrix defined as the product
 *>    of blocked elementary reflectors computed by short wide LQ
 *>    factorization (DGELQ)
+*>
 *> \endverbatim
 *
 *  Arguments:
 *  ==========
 *
 *> \param[in] SIDE
+*> \verbatim
 *>          SIDE is CHARACTER*1
 *>          = 'L': apply Q or Q**T from the Left;
 *>          = 'R': apply Q or Q**T from the Right.
+*> \endverbatim
 *>
 *> \param[in] TRANS
+*> \verbatim
 *>          TRANS is CHARACTER*1
 *>          = 'N':  No transpose, apply Q;
 *>          = 'T':  Transpose, apply Q**T.
+*> \endverbatim
+*>
 *> \param[in] M
 *> \verbatim
 *>          M is INTEGER
@@ -49,7 +56,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. N >= M.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,28 +64,29 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          M >= K >= 0;
+*>          If SIDE = 'L', M >= K >= 0;
+*>          if SIDE = 'R', N >= K >= 0.
 *>
 *> \endverbatim
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array, dimension (LDA,K)
-*>          Part of the data structure to represent Q as returned by ZGELQ.
+*>          A is DOUBLE PRECISION array, dimension
+*>                               (LDA,M) if SIDE = 'L',
+*>                               (LDA,N) if SIDE = 'R'
+*>          Part of the data structure to represent Q as returned by DGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.
-*>          If SIDE = 'L', LDA >= max(1,M);
-*>          if SIDE = 'R', LDA >= max(1,N).
+*>          The leading dimension of the array A. LDA >= max(1,K).
 *> \endverbatim
 *>
 *> \param[in] T
 *> \verbatim
 *>          T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE)).
-*>          Part of the data structure to represent Q as returned by ZGELQ.
+*>          Part of the data structure to represent Q as returned by DGELQ.
 *> \endverbatim
 *>
 *> \param[in] TSIZE
@@ -88,19 +96,23 @@
 *> \endverbatim
 *>
 *> \param[in,out] C
+*> \verbatim
 *>          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.
+*> \endverbatim
 *>
 *> \param[in] LDC
+*> \verbatim
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
 *>
 *> \param[out] WORK
 *> \verbatim
 *>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
-*>
 *> \endverbatim
+*>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
@@ -140,47 +152,49 @@
 *>
 *>          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
+*>          T(6: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 wide-and-short) or GELQT to compute
+*>  block sizes MB and NB returned by ILAENV, DGELQ will use either
+*>  DLASWLQ (if the matrix is wide-and-short) or DGELQT to compute
 *>  the LQ factorization.
-*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  This version of DGEMLQ will use either DLAMSWLQ or DGEMLQT to 
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LAMSWLQ or GEMLQT.
+*>  Further Details in DLAMSWLQ or DGEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
-     $        C, LDC, WORK, LWORK, INFO )
+     $                   C, LDC, WORK, LWORK, INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
+      CHARACTER          SIDE, TRANS
+      INTEGER            INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION  A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+      DOUBLE PRECISION   A( LDA, * ), T( * ), C( LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LEFT, RIGHT, TRAN, NOTRAN, LQUERY
-      INTEGER    I, II, KK, MB, NB, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           DLAMSWLQ, DGEMLQT, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,15 +202,15 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK.LT.0
+      LQUERY  = LWORK.EQ.-1
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'T' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(T(2))
-      NB = INT(T(3))
-      IF ( LEFT ) THEN
+      MB = INT( T( 2 ) )
+      NB = INT( T( 3 ) )
+      IF( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
@@ -204,7 +218,7 @@
         MN = N
       END IF
 *
-      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+      IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
         IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
           NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
@@ -216,34 +230,33 @@
 *
       INFO = 0
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
-         INFO = -1
+        INFO = -1
       ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
-         INFO = -2
+        INFO = -2
       ELSE IF( M.LT.0 ) THEN
         INFO = -3
-      ELSE IF( N.LT.0) THEN
+      ELSE IF( N.LT.0 ) THEN
         INFO = -4
-      ELSE IF( K.LT.0 ) THEN
+      ELSE IF( K.LT.0 .OR. K.GT.MN ) THEN
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
-     $         .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
         INFO = -11
-      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        WORK(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
 *
@@ -256,13 +269,13 @@
       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,
-     $        T(4), MB, C, LDC, WORK, INFO)
+     $                T( 6 ), MB, C, LDC, WORK, INFO )
       ELSE
-        CALL DLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $    MB, C, LDC, WORK, LWORK, INFO )
+        CALL DLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = LW
+      WORK( 1 ) = LW
 *
       RETURN
 *
diff --git a/SRC/sgelq.f b/SRC/sgelq.f
index adc606d9..1ae47d15 100644
--- a/SRC/sgelq.f
+++ b/SRC/sgelq.f
@@ -3,13 +3,13 @@
 *  ===========
 *
 *       SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       REAL  A( LDA, * ), T( * ), WORK( * )
+*       REAL              A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -120,7 +120,7 @@
 *>
 *> 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
+*> (trapezoidal) L has to be stored in the lower 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.
 *>
@@ -146,74 +146,73 @@
 *>
 *>          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
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           SLASWLQ or SGELQT
 *>
 *>  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
+*>  block sizes MB and NB returned by ILAENV, SGELQ will use either
+*>  SLASWLQ (if the matrix is short-and-wide) or SGELQT to compute
 *>  the LQ factorization.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-     $   INFO)
+     $                  INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
+      INTEGER            INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      REAL              A( LDA, * ), T( * ), WORK( * )
+      REAL               A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS, MINT, MINW
-      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
+      LOGICAL            LQUERY, LMINWS, MINT, MINW
+      INTEGER            MB, NB, MINTSZ, NBLCKS
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
-*     .. EXTERNAL SUBROUTINES ..
+*     ..
+*     .. External Subroutines ..
       EXTERNAL           SGELQT, SLASWLQ, XERBLA
-*     .. INTRINSIC FUNCTIONS ..
+*     ..
+*     .. Intrinsic Functions ..
       INTRINSIC          MAX, MIN, MOD
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       INTEGER            ILAENV
       EXTERNAL           ILAENV
 *     ..
-*     .. EXECUTABLE STATEMENTS ..
+*     .. Executable statements ..
 *
-*     TEST THE INPUT ARGUMENTS
+*     Test the input arguments
 *
       INFO = 0
 *
-      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+      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
+      IF( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) THEN
+        IF( TSIZE.NE.-1 ) MINT = .TRUE.
+        IF( LWORK.NE.-1 ) MINW = .TRUE.
+      END IF
 *
 *     Determine the block size
 *
-      IF ( MIN(M,N).GT.0 ) THEN
-        MB = ILAENV( 1, 'SGELQ ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'SGELQ ', ' ', M, N, 2, -1)
+      IF( MIN( M, N ).GT.0 ) THEN
+        MB = ILAENV( 1, 'SGELQ ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'SGELQ ', ' ', M, N, 2, -1 )
       ELSE
         MB = 1
         NB = N
@@ -235,16 +234,16 @@
 *
       LMINWS = .FALSE.
       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
-            NB = N
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.MINTSZ )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
+          LMINWS = .TRUE.
+          MB = 1
+          NB = N
         END IF
-        IF ( LWORK.LT.MB*M ) THEN
-            LMINWS = .TRUE.
-            MB = 1
+        IF( LWORK.LT.MB*M ) THEN
+          LMINWS = .TRUE.
+          MB = 1
         END IF
       END IF
 *
@@ -262,42 +261,43 @@
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        IF ( MINT ) THEN
-          T(1) = MINTSZ
+      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 )
+          T( 1 ) = MB*M*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, MB*M )
-        ENDIF
+          WORK( 1 ) = MAX( 1, MB*M )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGELQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N).EQ.0 ) THEN
+      IF( MIN( M, N ).EQ.0 ) THEN
         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, T(4), MB, WORK, INFO)
+        CALL SGELQT( M, N, MB, A, LDA, T( 6 ), MB, WORK, INFO )
       ELSE
-        CALL SLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
-     $                    LWORK, INFO)
+        CALL SLASWLQ( M, N, MB, NB, A, LDA, T( 6 ), MB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, MB*M )
+*
+      WORK( 1 ) = MAX( 1, MB*M )
       RETURN
 *
 *     End of SGELQ
diff --git a/SRC/sgemlq.f b/SRC/sgemlq.f
index a9cd54bd..42306ae4 100644
--- a/SRC/sgemlq.f
+++ b/SRC/sgemlq.f
@@ -3,22 +3,23 @@
 *  ===========
 *
 *      SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
-*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
+*     $                   TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
-*      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
+*      CHARACTER          SIDE, TRANS
+*      INTEGER            INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      REAL              A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*      REAL               A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*     ..
+*
 *> \par Purpose:
 *  =============
 *>
 *> \verbatim
 *>
-*>      SGEMLQ 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'
 *>    TRANS = 'N':      Q * C          C * Q
@@ -26,20 +27,26 @@
 *>    where Q is a real orthogonal matrix defined as the product
 *>    of blocked elementary reflectors computed by short wide LQ
 *>    factorization (SGELQ)
+*>
 *> \endverbatim
 *
 *  Arguments:
 *  ==========
 *
 *> \param[in] SIDE
+*> \verbatim
 *>          SIDE is CHARACTER*1
 *>          = 'L': apply Q or Q**T from the Left;
 *>          = 'R': apply Q or Q**T from the Right.
+*> \endverbatim
 *>
 *> \param[in] TRANS
+*> \verbatim
 *>          TRANS is CHARACTER*1
 *>          = 'N':  No transpose, apply Q;
 *>          = 'T':  Transpose, apply Q**T.
+*> \endverbatim
+*>
 *> \param[in] M
 *> \verbatim
 *>          M is INTEGER
@@ -49,7 +56,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. N >= M.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,28 +64,28 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          M >= K >= 0;
-*>
+*>          If SIDE = 'L', M >= K >= 0;
+*>          if SIDE = 'R', N >= K >= 0.
 *> \endverbatim
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is REAL array, dimension (LDA,K)
-*>          Part of the data structure to represent Q as returned by ZGELQ.
+*>          A is REAL array, dimension
+*>                               (LDA,M) if SIDE = 'L',
+*>                               (LDA,N) if SIDE = 'R'
+*>          Part of the data structure to represent Q as returned by DGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.
-*>          If SIDE = 'L', LDA >= max(1,M);
-*>          if SIDE = 'R', LDA >= max(1,N).
+*>          The leading dimension of the array A. LDA >= max(1,K).
 *> \endverbatim
 *>
 *> \param[in] T
 *> \verbatim
 *>          T is REAL array, dimension (MAX(5,TSIZE)).
-*>          Part of the data structure to represent Q as returned by ZGELQ.
+*>          Part of the data structure to represent Q as returned by SGELQ.
 *> \endverbatim
 *>
 *> \param[in] TSIZE
@@ -88,19 +95,23 @@
 *> \endverbatim
 *>
 *> \param[in,out] C
+*> \verbatim
 *>          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.
+*> \endverbatim
 *>
 *> \param[in] LDC
+*> \verbatim
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
 *>
 *> \param[out] WORK
 *> \verbatim
 *>         (workspace) REAL array, dimension (MAX(1,LWORK))
-*>
 *> \endverbatim
+*>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
@@ -140,47 +151,49 @@
 *>
 *>          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
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           SLASWLQ or SGELQT
 *>
 *>  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 wide-and-short) or GELQT to compute
+*>  block sizes MB and NB returned by ILAENV, SGELQ will use either
+*>  SLASWLQ (if the matrix is wide-and-short) or SGELQT to compute
 *>  the LQ factorization.
-*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  This version of SGEMLQ will use either SLAMSWLQ or SGEMLQT to 
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LAMSWLQ or GEMLQT.
+*>  Further Details in SLAMSWLQ or SGEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
-     $        C, LDC, WORK, LWORK, INFO )
+     $                   C, LDC, WORK, LWORK, INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
+      CHARACTER          SIDE, TRANS
+      INTEGER            INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      REAL              A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+      REAL               A( LDA, * ), T( * ), C( LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LEFT, RIGHT, TRAN, NOTRAN, LQUERY
-      INTEGER    I, II, KK, MB, NB, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           SLAMSWLQ, SGEMLQT, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,15 +201,15 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK.LT.0
+      LQUERY  = LWORK.EQ.-1
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'T' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(T(2))
-      NB = INT(T(3))
-      IF ( LEFT ) THEN
+      MB = INT( T( 2 ) )
+      NB = INT( T( 3 ) )
+      IF( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
@@ -204,7 +217,7 @@
         MN = N
       END IF
 *
-      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+      IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
         IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
           NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
@@ -216,34 +229,33 @@
 *
       INFO = 0
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
-         INFO = -1
+        INFO = -1
       ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
-         INFO = -2
+        INFO = -2
       ELSE IF( M.LT.0 ) THEN
         INFO = -3
-      ELSE IF( N.LT.0) THEN
+      ELSE IF( N.LT.0 ) THEN
         INFO = -4
-      ELSE IF( K.LT.0 ) THEN
+      ELSE IF( K.LT.0 .OR. K.GT.MN ) THEN
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
-     $         .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
         INFO = -11
-      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        WORK(1) = REAL(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
 *
@@ -256,13 +268,13 @@
       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,
-     $        T(4), MB, C, LDC, WORK, INFO)
+     $                T( 6 ), MB, C, LDC, WORK, INFO )
       ELSE
-        CALL SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $    MB, C, LDC, WORK, LWORK, INFO )
+        CALL SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = REAL(LW)
+      WORK( 1 ) = REAL( LW )
 *
       RETURN
 *
diff --git a/SRC/zgelq.f b/SRC/zgelq.f
index 5c51cf52..73d54771 100644
--- a/SRC/zgelq.f
+++ b/SRC/zgelq.f
@@ -3,13 +3,13 @@
 *  ===========
 *
 *       SUBROUTINE ZGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX*16  A( LDA, * ), T( * ), WORK( * )
+*       COMPLEX*16       A( LDA, * ), T( * ), WORK( * )
 *       ..
 *
 *
@@ -120,7 +120,7 @@
 *>
 *> 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
+*> (trapezoidal) L has to be stored in the lower 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.
 *>
@@ -146,74 +146,73 @@
 *>
 *>          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
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           ZLASWLQ or ZGELQT
 *>
 *>  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
+*>  block sizes MB and NB returned by ILAENV, ZGELQ will use either
+*>  ZLASWLQ (if the matrix is short-and-wide) or ZGELQT to compute
 *>  the LQ factorization.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE ZGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-     $   INFO)
+     $                  INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      INTEGER           INFO, LDA, M, N, TSIZE, LWORK
+      INTEGER            INFO, LDA, M, N, TSIZE, LWORK
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16  A( LDA, * ), T( * ), WORK( * )
+      COMPLEX*16         A( LDA, * ), T( * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LQUERY, LMINWS, MINT, MINW
-      INTEGER    MB, NB, I, II, KK, MINTSZ, NBLCKS
+      LOGICAL            LQUERY, LMINWS, MINT, MINW
+      INTEGER            MB, NB, MINTSZ, NBLCKS
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
-*     .. EXTERNAL SUBROUTINES ..
+*     ..
+*     .. External Subroutines ..
       EXTERNAL           ZGELQT, ZLASWLQ, XERBLA
-*     .. INTRINSIC FUNCTIONS ..
+*     ..
+*     .. Intrinsic Functions ..
       INTRINSIC          MAX, MIN, MOD
 *     ..
-*     .. EXTERNAL FUNCTIONS ..
+*     .. External Functions ..
       INTEGER            ILAENV
       EXTERNAL           ILAENV
 *     ..
-*     .. EXECUTABLE STATEMENTS ..
+*     .. Executable Statements ..
 *
-*     TEST THE INPUT ARGUMENTS
+*     Test the input arguments
 *
       INFO = 0
 *
-      LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR. 
+      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
+      IF( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) THEN
+        IF( TSIZE.NE.-1 ) MINT = .TRUE.
+        IF( LWORK.NE.-1 ) MINW = .TRUE.
+      END IF
 *
 *     Determine the block size
 *
-      IF ( MIN(M,N).GT.0 ) THEN
-        MB = ILAENV( 1, 'ZGELQ ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'ZGELQ ', ' ', M, N, 2, -1)
+      IF( MIN( M, N ).GT.0 ) THEN
+        MB = ILAENV( 1, 'ZGELQ ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'ZGELQ ', ' ', M, N, 2, -1 )
       ELSE
         MB = 1
         NB = N
@@ -235,14 +234,14 @@
 *
       LMINWS = .FALSE.
       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
+     $    .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.MINTSZ )
+     $    .AND. ( .NOT.LQUERY ) ) THEN
+        IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
             LMINWS = .TRUE.
             MB = 1
             NB = N
         END IF
-        IF ( LWORK.LT.MB*M ) THEN
+        IF( LWORK.LT.MB*M ) THEN
             LMINWS = .TRUE.
             MB = 1
         END IF
@@ -262,42 +261,44 @@
         INFO = -8
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        IF ( MINT ) THEN
-          T(1) = MINTSZ
+      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 )
+          T( 1 ) = MB*M*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, MB*M )
-        ENDIF
+          WORK( 1 ) = MAX( 1, MB*M )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGELQ', -INFO )
         RETURN
-      ELSE IF (LQUERY) THEN
+      ELSE IF( LQUERY ) THEN
         RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N).EQ.0 ) THEN
+      IF( MIN( M, N ).EQ.0 ) THEN
         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, T(4), MB, WORK, INFO)
+        CALL ZGELQT( M, N, MB, A, LDA, T( 6 ), MB, WORK, INFO )
       ELSE
-        CALL ZLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
-     $                    LWORK, INFO)
+        CALL ZLASWLQ( M, N, MB, NB, A, LDA, T( 6 ), MB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, MB*M )
+*
+      WORK( 1 ) = MAX( 1, MB*M )
+*
       RETURN
 *
 *     End of ZGELQ
diff --git a/SRC/zgemlq.f b/SRC/zgemlq.f
index f02d7b1a..5602d872 100644
--- a/SRC/zgemlq.f
+++ b/SRC/zgemlq.f
@@ -3,22 +3,21 @@
 *  ===========
 *
 *      SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
-*     $                     TSIZE, C, LDC, WORK, LWORK, INFO )
+*     $                   TSIZE, C, LDC, WORK, LWORK, INFO )
 *
 *
 *     .. Scalar Arguments ..
-*      CHARACTER         SIDE, TRANS
-*      INTEGER           INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
+*      CHARACTER          SIDE, TRANS
+*      INTEGER            INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-*      COMPLEX*16        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+*      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
@@ -26,20 +25,26 @@
 *>      where Q is a complex unitary matrix defined as the product
 *>      of blocked elementary reflectors computed by short wide
 *>      LQ factorization (ZGELQ)
+*>
 *> \endverbatim
 *
 *  Arguments:
 *  ==========
 *
 *> \param[in] SIDE
+*> \verbatim
 *>          SIDE is CHARACTER*1
 *>          = 'L': apply Q or Q**T from the Left;
 *>          = 'R': apply Q or Q**T from the Right.
+*> \endverbatim
 *>
 *> \param[in] TRANS
+*> \verbatim
 *>          TRANS is CHARACTER*1
 *>          = 'N':  No transpose, apply Q;
 *>          = 'T':  Transpose, apply Q**T.
+*> \endverbatim
+*>
 *> \param[in] M
 *> \verbatim
 *>          M is INTEGER
@@ -49,7 +54,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. N >= M.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,22 +62,23 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          M >= K >= 0;
+*>          If SIDE = 'L', M >= K >= 0;
+*>          if SIDE = 'R', N >= K >= 0.
 *>
 *> \endverbatim
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is COMPLEX*16 array, dimension (LDA,K)
+*>          A is COMPLEX*16 array, dimension
+*>                               (LDA,M) if SIDE = 'L',
+*>                               (LDA,N) if SIDE = 'R'
 *>          Part of the data structure to represent Q as returned by ZGELQ.
 *> \endverbatim
 *>
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.
-*>          If SIDE = 'L', LDA >= max(1,M);
-*>          if SIDE = 'R', LDA >= max(1,N).
+*>          The leading dimension of the array A. LDA >= max(1,K).
 *> \endverbatim
 *>
 *> \param[in] T
@@ -88,19 +94,23 @@
 *> \endverbatim
 *>
 *> \param[in,out] C
+*> \verbatim
 *>          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.
+*> \endverbatim
 *>
 *> \param[in] LDC
+*> \verbatim
 *>          LDC is INTEGER
 *>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
 *>
 *> \param[out] WORK
 *> \verbatim
 *>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
-*>
 *> \endverbatim
+*>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER
@@ -140,47 +150,49 @@
 *>
 *>          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
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           ZLASWLQ or ZGELQT
 *>
 *>  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 wide-and-short) or GELQT to compute
+*>  block sizes MB and NB returned by ILAENV, ZGELQ will use either
+*>  ZLASWLQ (if the matrix is wide-and-short) or ZGELQT to compute
 *>  the LQ factorization.
-*>  This version of GEMLQ will use either LAMSWLQ or GEMLQT to 
+*>  This version of ZGEMLQ will use either ZLAMSWLQ or ZGEMLQT to 
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LAMSWLQ or GEMLQT.
+*>  Further Details in ZLAMSWLQ or ZGEMLQT.
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
-     $        C, LDC, WORK, LWORK, INFO )
+     $                   C, LDC, WORK, LWORK, INFO )
 *
-*  -- LAPACK computational routine (version 3.5.0) --
+*  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2013
+*     November 2016
 *
 *     .. Scalar Arguments ..
-      CHARACTER         SIDE, TRANS
-      INTEGER           INFO, LDA, M, N, K, TSIZE, LWORK, LDC
+      CHARACTER          SIDE, TRANS
+      INTEGER            INFO, LDA, M, N, K, TSIZE, LWORK, LDC
 *     ..
 *     .. Array Arguments ..
-      COMPLEX*16        A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
+      COMPLEX*16         A( LDA, * ), T( * ), C( LDC, * ), WORK( * )
 *     ..
 *
 * =====================================================================
 *
 *     ..
 *     .. Local Scalars ..
-      LOGICAL    LEFT, RIGHT, TRAN, NOTRAN, LQUERY
-      INTEGER    I, II, KK, MB, NB, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           ZLAMSWLQ, ZGEMLQT, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,15 +200,15 @@
 *
 *     Test the input arguments
 *
-      LQUERY  = LWORK.LT.0
+      LQUERY  = LWORK.EQ.-1
       NOTRAN  = LSAME( TRANS, 'N' )
       TRAN    = LSAME( TRANS, 'C' )
       LEFT    = LSAME( SIDE, 'L' )
       RIGHT   = LSAME( SIDE, 'R' )
 *
-      MB = INT(T(2))
-      NB = INT(T(3))
-      IF ( LEFT ) THEN
+      MB = INT( T( 2 ) )
+      NB = INT( T( 3 ) )
+      IF( LEFT ) THEN
         LW = N * MB
         MN = M
       ELSE
@@ -204,7 +216,7 @@
         MN = N
       END IF
 *
-      IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+      IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
         IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
           NBLCKS = ( MN - K ) / ( NB - K )
         ELSE
@@ -216,34 +228,33 @@
 *
       INFO = 0
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
-         INFO = -1
+        INFO = -1
       ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
-         INFO = -2
+        INFO = -2
       ELSE IF( M.LT.0 ) THEN
         INFO = -3
-      ELSE IF( N.LT.0) THEN
+      ELSE IF( N.LT.0 ) THEN
         INFO = -4
-      ELSE IF( K.LT.0 ) THEN
+      ELSE IF( K.LT.0 .OR. K.GT.MN ) THEN
         INFO = -5
       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
-     $         .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
         INFO = -11
-      ELSE IF((  LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-      IF( INFO.EQ.0 )  THEN
-        WORK(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
 *
@@ -256,13 +267,13 @@
       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,
-     $        T(4), MB, C, LDC, WORK, INFO)
+     $                T( 6 ), MB, C, LDC, WORK, INFO )
       ELSE
-        CALL ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $    MB, C, LDC, WORK, LWORK, INFO )
+        CALL ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 MB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = LW
+      WORK( 1 ) = LW
 *
       RETURN
 *