Fix ?GEQR and ?GEMQR

8 files changed, 507 insertions, 451 deletions

diff --git a/SRC/cgemqr.f b/SRC/cgemqr.f
index a5d6263f..a5420553 100644
--- a/SRC/cgemqr.f
+++ b/SRC/cgemqr.f
@@ -3,43 +3,51 @@
 *  ===========
 *
 *      SUBROUTINE CGEMQR( 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:
 *  =============
 *>
 *> \verbatim
 *>
-*>      CGEMQR overwrites the general real M-by-N matrix C with
-*>
+*> CGEMQR overwrites the general real M-by-N matrix C with
 *>
 *>                      SIDE = 'L'     SIDE = 'R'
 *>      TRANS = 'N':      Q * C          C * Q
 *>      TRANS = 'T':      Q**H * C       C * Q**H
-*>      where Q is a complex unitary matrix defined as the product
-*>      of blocked elementary reflectors computed by tall skinny
-*>      QR factorization (CGEQR)
+*>
+*> where Q is a complex unitary matrix defined as the product
+*> of blocked elementary reflectors computed by tall skinny
+*> QR factorization (CGEQR)
+*>
 *> \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 +57,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. M >= N >= 0.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,14 +65,14 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          N >= 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 SGEQR.
+*>          Part of the data structure to represent Q as returned by CGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,7 +86,7 @@
 *> \param[in] T
 *> \verbatim
 *>          T is COMPLEX array, dimension (MAX(5,TSIZE)).
-*>          Part of the data structure to represent Q as returned by SGEQR.
+*>          Part of the data structure to represent Q as returned by CGEQR.
 *> \endverbatim
 *>
 *> \param[in] TSIZE
@@ -88,19 +96,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 +152,50 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
-*>                           LATSQR or GEQRT
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           CLATSQR or CGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, CGEQR will use either
+*>  CLATSQR (if the matrix is tall-and-skinny) or CGEQRT to compute
 *>  the QR factorization.
-*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  This version of CGEMQR will use either CLAMTSQR or CGEMQRT to
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LATMSQR or GEMQRT.
+*>  Further Details in CLAMTSQR or CGEMQRT.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE CGEMQR( 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    MB, NB, I, II, KK, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           CGEMQRT, CLAMTSQR, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,27 +203,27 @@
 *
 *     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 * NB
         MN = M
-      ELSE IF(RIGHT) THEN
+      ELSE
         LW = MB * NB
         MN = N
       END IF
 *
-      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
-        IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
-           NBLCKS = ( MN - K ) / ( MB - K )
+      IF( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, MB - K ).EQ.0 ) THEN
+          NBLCKS = ( MN - K ) / ( MB - K )
         ELSE
-           NBLCKS = ( MN - K ) / ( MB - K ) + 1
+          NBLCKS = ( MN - K ) / ( MB - K ) + 1
         END IF
       ELSE
         NBLCKS = 1
@@ -216,55 +231,52 @@
 *
       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
         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
+      ELSE IF( LDA.LT.MAX( 1, MN ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
-     $        .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
-         INFO = -11
+      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+        INFO = -11
       ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-*     Determine the block size if it is tall skinny or short and wide
-*
-      IF( INFO.EQ.0 )  THEN
-        WORK(1) = LW
+      IF( INFO.EQ.0 ) THEN
+        WORK( 1 ) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGEMQR', -INFO )
         RETURN
-      ELSE IF ( LQUERY ) THEN
-       RETURN
+      ELSE IF( LQUERY ) THEN
+        RETURN
       END IF
 *
 *     Quick return if possible
 *
-      IF( MIN(M,N,K).EQ.0 ) THEN
+      IF( MIN( M, N, K ).EQ.0 ) THEN
         RETURN
       END IF
 *
       IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
      $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
-        CALL CGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        T(4), NB, C, LDC, WORK, INFO )
+        CALL CGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA, T( 6 ),
+     $                NB, C, LDC, WORK, INFO )
       ELSE
-        CALL CLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $      NB, C, LDC, WORK, LWORK, INFO )
+        CALL CLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = LW
+      WORK( 1 ) = LW
 *
       RETURN
 *
diff --git a/SRC/cgeqr.f b/SRC/cgeqr.f
index 1743960f..e336a916 100644
--- a/SRC/cgeqr.f
+++ b/SRC/cgeqr.f
@@ -3,7 +3,7 @@
 *  ===========
 *
 *       SUBROUTINE CGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
@@ -146,74 +146,74 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
-*>                           DLATSQR or DGEQRT
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           CLATSQR or CGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, CGEQR will use either
+*>  CLATSQR (if the matrix is tall-and-skinny) or CGEQRT to compute
 *>  the QR factorization.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE CGEQR( 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           CLATSQR, CGEQRT, 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, 'CGEQR ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'CGEQR ', ' ', M, N, 2, -1)
+      IF( MIN( M, N ).GT.0 ) THEN
+        MB = ILAENV( 1, 'CGEQR ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'CGEQR ', ' ', M, N, 2, -1 )
       ELSE
         MB = M
         NB = 1
@@ -221,7 +221,7 @@
       IF( MB.GT.M .OR. MB.LE.N ) MB = M
       IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
       MINTSZ = N + 5
-      IF ( MB.GT.N .AND. M.GT.N ) THEN
+      IF( MB.GT.N .AND. M.GT.N ) THEN
         IF( MOD( M - N, MB - N ).EQ.0 ) THEN
           NBLCKS = ( M - N ) / ( MB - N )
         ELSE
@@ -235,16 +235,16 @@
 *
       LMINWS = .FALSE.
       IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
-     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.MINTSZ )
      $    .AND. ( .NOT.LQUERY ) ) THEN
-        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
-            LMINWS = .TRUE.
-            NB = 1
-            MB = M
+        IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
+          LMINWS = .TRUE.
+          NB = 1
+          MB = M
         END IF
-        IF ( LWORK.LT.NB*N ) THEN
-            LMINWS = .TRUE.
-            NB = 1
+        IF( LWORK.LT.NB*N ) THEN
+          LMINWS = .TRUE.
+          NB = 1
         END IF
       END IF
 *
@@ -262,42 +262,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) = NB*N*NBLCKS + 5
-        ENDIF
-        T(2) = MB
-        T(3) = NB
-        IF ( MINW ) THEN
-          WORK(1) = MAX( 1, N )
+          T( 1 ) = NB*N*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, NB*N )
-        ENDIF
+          WORK( 1 ) = MAX( 1, NB*N )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGEQR', -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 QR Decomposition
 *
       IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
-        CALL CGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
+        CALL CGEQRT( M, N, NB, A, LDA, T( 6 ), NB, WORK, INFO )
       ELSE
-        CALL CLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
-     $                    LWORK, INFO )
+        CALL CLATSQR( M, N, MB, NB, A, LDA, T( 6 ), NB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, NB*N )
+*
+      WORK( 1 ) = MAX( 1, NB*N )
+*
       RETURN
 *
 *     End of CGEQR
diff --git a/SRC/dgemqr.f b/SRC/dgemqr.f
index ee02d98e..6032be6d 100644
--- a/SRC/dgemqr.f
+++ b/SRC/dgemqr.f
@@ -3,43 +3,51 @@
 *  ===========
 *
 *      SUBROUTINE DGEMQR( 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:
 *  =============
 *>
 *> \verbatim
 *>
-*>      DGEMQR overwrites the general real M-by-N matrix C with
-*>
+*> DGEMQR 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 tall skinny
-*>      QR factorization (DGEQR)
+*>
+*> where Q is a real orthogonal matrix defined as the product
+*> of blocked elementary reflectors computed by tall skinny
+*> QR factorization (DGEQR)
+*>
 *> \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 +57,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. M >= N >= 0.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,14 +65,14 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          N >= 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 SGEQR.
+*>          Part of the data structure to represent Q as returned by DGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,7 +86,7 @@
 *> \param[in] T
 *> \verbatim
 *>          T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE)).
-*>          Part of the data structure to represent Q as returned by SGEQR.
+*>          Part of the data structure to represent Q as returned by DGEQR.
 *> \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,50 @@
 *>
 *>          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
 *>                           DLATSQR or DGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, DGEQR will use either
+*>  DLATSQR (if the matrix is tall-and-skinny) or DGEQRT to compute
 *>  the QR factorization.
-*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  This version of DGEMQR will use either DLAMTSQR or DGEMQRT to 
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LATMSQR or GEMQRT.
+*>  Further Details in DLATMSQR or DGEMQRT.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE DGEMQR( 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    MB, NB, I, II, KK, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           DGEMQRT, DLAMTSQR, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,65 +203,62 @@
 *
 *     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 * NB
         MN = M
-      ELSE IF(RIGHT) THEN
+      ELSE
         LW = MB * NB
         MN = N
       END IF
 *
-      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
-          IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
-             NBLCKS = ( MN - K ) / ( MB - K )
-          ELSE
-             NBLCKS = ( MN - K ) / ( MB - K ) + 1
-          END IF
+      IF( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, MB - K ).EQ.0 ) THEN
+          NBLCKS = ( MN - K ) / ( MB - K )
+        ELSE
+          NBLCKS = ( MN - K ) / ( MB - K ) + 1
+        END IF
       ELSE
         NBLCKS = 1
       END IF
 *
       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
         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
+      ELSE IF( LDA.LT.MAX( 1, MN ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
-     $        .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
-         INFO = -11
+      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+        INFO = -11
       ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-*     Determine the block size if it is tall skinny or short and wide
-*
-      IF( INFO.EQ.0 )  THEN
-        WORK(1) = LW
+      IF( INFO.EQ.0 ) THEN
+        WORK( 1 ) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGEMQR', -INFO )
         RETURN
-      ELSE IF ( LQUERY ) THEN
-       RETURN
+      ELSE IF( LQUERY ) THEN
+        RETURN
       END IF
 *
 *     Quick return if possible
@@ -257,14 +269,14 @@
 *
       IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
      $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
-        CALL DGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        T(4), NB, C, LDC, WORK, INFO )
+        CALL DGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA, T( 6 ),
+     $                NB, C, LDC, WORK, INFO )
       ELSE
-        CALL DLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $      NB, C, LDC, WORK, LWORK, INFO )
+        CALL DLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = LW
+      WORK( 1 ) = LW
 *
       RETURN
 *
diff --git a/SRC/dgeqr.f b/SRC/dgeqr.f
index 93f02bd8..5212c69e 100644
--- a/SRC/dgeqr.f
+++ b/SRC/dgeqr.f
@@ -3,7 +3,7 @@
 *  ===========
 *
 *       SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
@@ -146,72 +146,72 @@
 *>
 *>          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
 *>                           DLATSQR or DGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, DGEQR will use either
+*>  DLATSQR (if the matrix is tall-and-skinny) or DGEQRT to compute
 *>  the QR factorization.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE DGEQR( 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           DLATSQR, DGEQRT, 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, 'DGEQR ', ' ', M, N, 1, -1 )
         NB = ILAENV( 1, 'DGEQR ', ' ', M, N, 2, -1 )
       ELSE
@@ -221,7 +221,7 @@
       IF( MB.GT.M .OR. MB.LE.N ) MB = M
       IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
       MINTSZ = N + 5
-      IF ( MB.GT.N .AND. M.GT.N ) THEN
+      IF( MB.GT.N .AND. M.GT.N ) THEN
         IF( MOD( M - N, MB - N ).EQ.0 ) THEN
           NBLCKS = ( M - N ) / ( MB - N )
         ELSE
@@ -235,16 +235,16 @@
 *
       LMINWS = .FALSE.
       IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
-     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.MINTSZ )
      $    .AND. ( .NOT.LQUERY ) ) THEN
-        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
-            LMINWS = .TRUE.
-            NB = 1
-            MB = M
+        IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
+          LMINWS = .TRUE.
+          NB = 1
+          MB = M
         END IF
-        IF ( LWORK.LT.NB*N ) THEN
-            LMINWS = .TRUE.
-            NB = 1
+        IF( LWORK.LT.NB*N ) THEN
+          LMINWS = .TRUE.
+          NB = 1
         END IF
       END IF
 *
@@ -262,42 +262,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) = NB*N*NBLCKS + 5
-        ENDIF
-        T(2) = MB
-        T(3) = NB
-        IF ( MINW ) THEN
-          WORK(1) = MAX( 1, N )
+          T( 1 ) = NB*N*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, NB*N )
-        ENDIF
+          WORK( 1 ) = MAX( 1, NB*N )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGEQR', -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 QR Decomposition
 *
       IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
-        CALL DGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
+        CALL DGEQRT( M, N, NB, A, LDA, T( 6 ), NB, WORK, INFO )
       ELSE
-        CALL DLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
-     $                    LWORK, INFO )
+        CALL DLATSQR( M, N, MB, NB, A, LDA, T( 6 ), NB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, NB*N )
+*
+      WORK( 1 ) = MAX( 1, NB*N )
+*
       RETURN
 *
 *     End of DGEQR
diff --git a/SRC/sgemqr.f b/SRC/sgemqr.f
index 78eb75c8..6bf3a1bd 100644
--- a/SRC/sgemqr.f
+++ b/SRC/sgemqr.f
@@ -3,43 +3,51 @@
 *  ===========
 *
 *      SUBROUTINE SGEMQR( 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
 *>
-*>      SGEMQR overwrites the general real M-by-N matrix C with
-*>
+*> SGEMQR 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 tall skinny
-*>      QR factorization (SGEQR)
+*>     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 tall skinny
+*> QR factorization (SGEQR)
+*>
 *> \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 +57,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. M >= N >= 0.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,8 +65,8 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          N >= K >= 0;
-*>
+*>          If SIDE = 'L', M >= K >= 0;
+*>          if SIDE = 'R', N >= K >= 0.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -88,19 +96,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 +152,50 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
-*>                           LATSQR or GEQRT
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           SLATSQR or SGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, SGEQR will use either
+*>  SLATSQR (if the matrix is tall-and-skinny) or SGEQRT to compute
 *>  the QR factorization.
-*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  This version of SGEMQR will use either SLAMTSQR or SGEMQRT to 
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LATMSQR or GEMQRT.
+*>  Further Details in SLAMTSQR or SGEMQRT.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE SGEMQR( 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    MB, NB, I, II, KK, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           SGEMQRT, SLAMTSQR, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,65 +203,62 @@
 *
 *     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 * NB
         MN = M
-      ELSE IF(RIGHT) THEN
+      ELSE
         LW = MB * NB
         MN = N
       END IF
 *
-      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
-          IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
-             NBLCKS = ( MN - K ) / ( MB - K )
-          ELSE
-             NBLCKS = ( MN - K ) / ( MB - K ) + 1
-          END IF
+      IF( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, MB - K ).EQ.0 ) THEN
+          NBLCKS = ( MN - K ) / ( MB - K )
+        ELSE
+          NBLCKS = ( MN - K ) / ( MB - K ) + 1
+        END IF
       ELSE
         NBLCKS = 1
       END IF
 *
       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
         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
+      ELSE IF( LDA.LT.MAX( 1, MN ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
-     $        .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
-         INFO = -11
+      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+        INFO = -11
       ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-*     Determine the block size if it is tall skinny or short and wide
-*
-      IF( INFO.EQ.0 )  THEN
-        WORK(1) = LW
+      IF( INFO.EQ.0 ) THEN
+        WORK( 1 ) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGEMQR', -INFO )
         RETURN
-      ELSE IF ( LQUERY ) THEN
-       RETURN
+      ELSE IF( LQUERY ) THEN
+        RETURN
       END IF
 *
 *     Quick return if possible
@@ -257,14 +269,14 @@
 *
       IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
      $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
-        CALL SGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        T(4), NB, C, LDC, WORK, INFO )
+        CALL SGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA, T( 6 ),
+     $                NB, C, LDC, WORK, INFO )
       ELSE
-        CALL SLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $      NB, C, LDC, WORK, LWORK, INFO )
+        CALL SLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = LW
+      WORK( 1 ) = LW
 *
       RETURN
 *
diff --git a/SRC/sgeqr.f b/SRC/sgeqr.f
index ad18e8a2..3f561a45 100644
--- a/SRC/sgeqr.f
+++ b/SRC/sgeqr.f
@@ -3,7 +3,7 @@
 *  ===========
 *
 *       SUBROUTINE SGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
@@ -146,74 +146,74 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
-*>                           DLATSQR or DGEQRT
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           SLATSQR or SGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, SGEQR will use either
+*>  SLATSQR (if the matrix is tall-and-skinny) or SGEQRT to compute
 *>  the QR factorization.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE SGEQR( 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           SLATSQR, SGEQRT, 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, 'SGEQR ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'SGEQR ', ' ', M, N, 2, -1)
+      IF( MIN( M, N ).GT.0 ) THEN
+        MB = ILAENV( 1, 'SGEQR ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'SGEQR ', ' ', M, N, 2, -1 )
       ELSE
         MB = M
         NB = 1
@@ -235,16 +235,16 @@
 *
       LMINWS = .FALSE.
       IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
-     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.MINTSZ )
      $    .AND. ( .NOT.LQUERY ) ) THEN
-        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
-            LMINWS = .TRUE.
-            NB = 1
-            MB = M
+        IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
+          LMINWS = .TRUE.
+          NB = 1
+          MB = M
         END IF
-        IF ( LWORK.LT.NB*N ) THEN
-            LMINWS = .TRUE.
-            NB = 1
+        IF( LWORK.LT.NB*N ) THEN
+          LMINWS = .TRUE.
+          NB = 1
         END IF
       END IF
 *
@@ -262,42 +262,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) = NB*N*NBLCKS + 5
-        ENDIF
-        T(2) = MB
-        T(3) = NB
-        IF ( MINW ) THEN
-          WORK(1) = MAX( 1, N )
+          T( 1 ) = NB*N*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, NB*N )
-        ENDIF
+          WORK( 1 ) = MAX( 1, NB*N )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGEQR', -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 QR Decomposition
 *
       IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
-        CALL SGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
+        CALL SGEQRT( M, N, NB, A, LDA, T( 6 ), NB, WORK, INFO )
       ELSE
-        CALL SLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
-     $                    LWORK, INFO )
+        CALL SLATSQR( M, N, MB, NB, A, LDA, T( 6 ), NB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, NB*N )
+*
+      WORK( 1 ) = MAX( 1, NB*N )
+*
       RETURN
 *
 *     End of SGEQR
diff --git a/SRC/zgemqr.f b/SRC/zgemqr.f
index 71614d5e..ed67b45f 100644
--- a/SRC/zgemqr.f
+++ b/SRC/zgemqr.f
@@ -3,43 +3,51 @@
 *  ===========
 *
 *      SUBROUTINE ZGEMQR( 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
 *>
-*>      ZGEMQR overwrites the general real M-by-N matrix C with
-*>
+*> ZGEMQR overwrites the general real M-by-N matrix C with
 *>
 *>                      SIDE = 'L'     SIDE = 'R'
 *>      TRANS = 'N':      Q * C          C * Q
 *>      TRANS = 'T':      Q**H * C       C * Q**H
-*>      where Q is a complex unitary matrix defined as the product
-*>      of blocked elementary reflectors computed by tall skinny
-*>      QR factorization (ZGEQR)
+*>
+*> where Q is a complex unitary matrix defined as the product
+*> of blocked elementary reflectors computed by tall skinny
+*> QR factorization (ZGEQR)
+*>
 *> \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 +57,7 @@
 *> \param[in] N
 *> \verbatim
 *>          N is INTEGER
-*>          The number of columns of the matrix C. M >= N >= 0.
+*>          The number of columns of the matrix C. N >= 0.
 *> \endverbatim
 *>
 *> \param[in] K
@@ -57,14 +65,14 @@
 *>          K is INTEGER
 *>          The number of elementary reflectors whose product defines
 *>          the matrix Q.
-*>          N >= 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)
-*>          Part of the data structure to represent Q as returned by SGEQR.
+*>          Part of the data structure to represent Q as returned by ZGEQR.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -78,7 +86,7 @@
 *> \param[in] T
 *> \verbatim
 *>          T is COMPLEX*16 array, dimension (MAX(5,TSIZE)).
-*>          Part of the data structure to represent Q as returned by SGEQR.
+*>          Part of the data structure to represent Q as returned by ZGEQR.
 *> \endverbatim
 *>
 *> \param[in] TSIZE
@@ -88,19 +96,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 +152,50 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
-*>                           LATSQR or GEQRT
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           ZLATSQR or ZGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, ZGEQR will use either
+*>  ZLATSQR (if the matrix is tall-and-skinny) or ZGEQRT to compute
 *>  the QR factorization.
-*>  This version of GEMQR will use either LAMTSQR or GEMQRT to 
+*>  This version of ZGEMQR will use either ZLAMTSQR or ZGEMQRT to
 *>  multiply matrix Q by another matrix.
-*>  Further Details in LATMSQR or GEMQRT.
+*>  Further Details in ZLAMTSQR or ZGEMQRT.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE ZGEMQR( 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    MB, NB, I, II, KK, LW, NBLCKS, MN
+      LOGICAL            LEFT, RIGHT, TRAN, NOTRAN, LQUERY
+      INTEGER            MB, NB, LW, NBLCKS, MN
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           ZGEMQRT, ZLAMTSQR, XERBLA
+*     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, MOD
 *     ..
@@ -188,65 +203,62 @@
 *
 *     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 * NB
         MN = M
-      ELSE IF(RIGHT) THEN
+      ELSE
         LW = MB * NB
         MN = N
       END IF
 *
-      IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
-          IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
-             NBLCKS = ( MN - K ) / ( MB - K )
-          ELSE
-             NBLCKS = ( MN - K ) / ( MB - K ) + 1
-          END IF
+      IF( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+        IF( MOD( MN - K, MB - K ).EQ.0 ) THEN
+          NBLCKS = ( MN - K ) / ( MB - K )
+        ELSE
+          NBLCKS = ( MN - K ) / ( MB - K ) + 1
+        END IF
       ELSE
         NBLCKS = 1
       END IF
 *
       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
         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
+      ELSE IF( LDA.LT.MAX( 1, MN ) ) THEN
         INFO = -7
-      ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
-     $        .AND. ( .NOT.LQUERY ) ) THEN
+      ELSE IF( TSIZE.LT.5 ) THEN
         INFO = -9
-      ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
-         INFO = -11
+      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+        INFO = -11
       ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
         INFO = -13
       END IF
 *
-*     Determine the block size if it is tall skinny or short and wide
-*
-      IF( INFO.EQ.0 )  THEN
-        WORK(1) = LW
+      IF( INFO.EQ.0 ) THEN
+        WORK( 1 ) = LW
       END IF
 *
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGEMQR', -INFO )
         RETURN
-      ELSE IF ( LQUERY ) THEN
-       RETURN
+      ELSE IF( LQUERY ) THEN
+        RETURN
       END IF
 *
 *     Quick return if possible
@@ -257,14 +269,14 @@
 *
       IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
      $     .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
-        CALL ZGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
-     $        T(4), NB, C, LDC, WORK, INFO )
+        CALL ZGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA, T( 6 ),
+     $                NB, C, LDC, WORK, INFO )
       ELSE
-        CALL ZLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
-     $      NB, C, LDC, WORK, LWORK, INFO )
+        CALL ZLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T( 6 ),
+     $                 NB, C, LDC, WORK, LWORK, INFO )
       END IF
 *
-      WORK(1) = LW
+      WORK( 1 ) = LW
 *
       RETURN
 *
diff --git a/SRC/zgeqr.f b/SRC/zgeqr.f
index 226db33e..a38e47b5 100644
--- a/SRC/zgeqr.f
+++ b/SRC/zgeqr.f
@@ -3,7 +3,7 @@
 *  ===========
 *
 *       SUBROUTINE ZGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
-*                        INFO)
+*                         INFO )
 *
 *       .. Scalar Arguments ..
 *       INTEGER           INFO, LDA, M, N, TSIZE, LWORK
@@ -146,74 +146,74 @@
 *>
 *>          T(2): row block size (MB)
 *>          T(3): column block size (NB)
-*>          T(4:TSIZE): data structure needed for Q, computed by
-*>                           DLATSQR or DGEQRT
+*>          T(6:TSIZE): data structure needed for Q, computed by
+*>                           ZLATSQR or ZGEQRT
 *>
 *>  Depending on the matrix dimensions M and N, and row and column
-*>  block sizes MB and NB returned by ILAENV, GEQR will use either
-*>  LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
+*>  block sizes MB and NB returned by ILAENV, ZGEQR will use either
+*>  ZLATSQR (if the matrix is tall-and-skinny) or ZGEQRT to compute
 *>  the QR factorization.
+*>
 *> \endverbatim
 *>
 *  =====================================================================
       SUBROUTINE ZGEQR( 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    ZLATSQR, ZGEQRT, XERBLA
-*     .. INTRINSIC FUNCTIONS ..
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZLATSQR, ZGEQRT, XERBLA
+*     ..
+*     .. 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, 'ZGEQR ', ' ', M, N, 1, -1)
-        NB = ILAENV( 1, 'ZGEQR ', ' ', M, N, 2, -1)
+      IF( MIN ( M, N ).GT.0 ) THEN
+        MB = ILAENV( 1, 'ZGEQR ', ' ', M, N, 1, -1 )
+        NB = ILAENV( 1, 'ZGEQR ', ' ', M, N, 2, -1 )
       ELSE
         MB = M
         NB = 1
@@ -221,7 +221,7 @@
       IF( MB.GT.M .OR. MB.LE.N ) MB = M
       IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
       MINTSZ = N + 5
-      IF ( MB.GT.N .AND. M.GT.N ) THEN
+      IF( MB.GT.N .AND. M.GT.N ) THEN
         IF( MOD( M - N, MB - N ).EQ.0 ) THEN
           NBLCKS = ( M - N ) / ( MB - N )
         ELSE
@@ -235,16 +235,16 @@
 *
       LMINWS = .FALSE.
       IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
-     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+     $    .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.MINTSZ )
      $    .AND. ( .NOT.LQUERY ) ) THEN
-        IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
-            LMINWS = .TRUE.
-            NB = 1
-            MB = M
+        IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
+          LMINWS = .TRUE.
+          NB = 1
+          MB = M
         END IF
-        IF ( LWORK.LT.NB*N ) THEN
-            LMINWS = .TRUE.
-            NB = 1
+        IF( LWORK.LT.NB*N ) THEN
+          LMINWS = .TRUE.
+          NB = 1
         END IF
       END IF
 *
@@ -262,42 +262,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) = NB*N*NBLCKS + 5
-        ENDIF
-        T(2) = MB
-        T(3) = NB
-        IF ( MINW ) THEN
-          WORK(1) = MAX( 1, N )
+          T( 1 ) = NB*N*NBLCKS + 5
+        END IF
+        T( 2 ) = MB
+        T( 3 ) = NB
+        IF( MINW ) THEN
+          WORK( 1 ) = MAX( 1, N )
         ELSE
-          WORK(1) = MAX( 1, NB*N )
-        ENDIF
+          WORK( 1 ) = MAX( 1, NB*N )
+        END IF
       END IF
       IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZGEQR', -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 QR Decomposition
 *
       IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
-        CALL ZGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
+        CALL ZGEQRT( M, N, NB, A, LDA, T( 6 ), NB, WORK, INFO )
       ELSE
-        CALL ZLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
-     $                    LWORK, INFO )
+        CALL ZLATSQR( M, N, MB, NB, A, LDA, T( 6 ), NB, WORK,
+     $                LWORK, INFO )
       END IF
-      WORK(1) = MAX( 1, NB*N )
+*
+      WORK( 1 ) = MAX( 1, NB*N )
+*
       RETURN
 *
 *     End of ZGEQR