Merge pull request #136 from echeresh/370_fix

A number of fixes for routines and testing sources integrated in LAPACK 3.7.0

8 files changed, 391 insertions, 220 deletions

diff --git a/SRC/chb2st_kernels.f b/SRC/chb2st_kernels.f
index 77ddaed5..fef90ad3 100644
--- a/SRC/chb2st_kernels.f
+++ b/SRC/chb2st_kernels.f
@@ -47,45 +47,87 @@
 *  Arguments:
 *  ==========
 *
-*> @param[in] n
-*>          The order of the matrix A.
-*>
-*> @param[in] nb
-*>          The size of the band.
-*>
-*> @param[in, out] A
-*>          A pointer to the matrix A.
-*>
-*> @param[in] lda
-*>          The leading dimension of the matrix A.
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*> \endverbatim
 *>
-*> @param[out] V
-*>          COMPLEX array, dimension 2*n if eigenvalues only are
-*>          requested or to be queried for vectors.
+*> \param[in] WANTZ
+*> \verbatim
+*>          WANTZ is LOGICAL which indicate if Eigenvalue are requested or both
+*>          Eigenvalue/Eigenvectors.
+*> \endverbatim
 *>
-*> @param[out] TAU
-*>          COMPLEX array, dimension (2*n).
-*>          The scalar factors of the Householder reflectors are stored
-*>          in this array.
+*> \param[in] TTYPE
+*> \verbatim
+*>          TTYPE is INTEGER
+*> \endverbatim
 *>
-*> @param[in] st
+*> \param[in] ST
+*> \verbatim
+*>          ST is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] ed
+*> \param[in] ED
+*> \verbatim
+*>          ED is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] sweep
+*> \param[in] SWEEP
+*> \verbatim
+*>          SWEEP is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] Vblksiz
-*>          internal parameter for indices.
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER. The order of the matrix A.
+*> \endverbatim
 *>
-*> @param[in] wantz
-*>          logical which indicate if Eigenvalue are requested or both
-*>          Eigenvalue/Eigenvectors.
+*> \param[in] NB
+*> \verbatim
+*>          NB is INTEGER. The size of the band.
+*> \endverbatim
+*>
+*> \param[in] IB
+*> \verbatim
+*>          IB is INTEGER.
+*> \endverbatim
+*>
+*> \param[in, out] A
+*> \verbatim
+*>          A is COMPLEX array. A pointer to the matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER. The leading dimension of the matrix A.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*>          V is COMPLEX array, dimension 2*n if eigenvalues only are
+*>          requested or to be queried for vectors.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is COMPLEX array, dimension (2*n).
+*>          The scalar factors of the Householder reflectors are stored
+*>          in this array.
+*> \endverbatim
+*>
+*> \param[in] LDVT
+*> \verbatim
+*>          LDVT is INTEGER.
+*> \endverbatim
 *>
-*> @param[in] work
-*>          Workspace of size nb.
+*> \param[in] WORK
+*> \verbatim
+*>          WORK is COMPLEX array. Workspace of size nb.
+*> \endverbatim
 *>
 *> \par Further Details:
 *  =====================
diff --git a/SRC/clasyf_aa.f b/SRC/clasyf_aa.f
index 2c8cdc46..860a0234 100644
--- a/SRC/clasyf_aa.f
+++ b/SRC/clasyf_aa.f
@@ -99,12 +99,12 @@
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.  LDA >= max(1,N).
+*>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
 *> \param[out] IPIV
 *> \verbatim
-*>          IPIV is INTEGER array, dimension (N)
+*>          IPIV is INTEGER array, dimension (M)
 *>          Details of the row and column interchanges,
 *>          the row and column k were interchanged with the row and
 *>          column IPIV(k).
@@ -217,8 +217,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
-*         where H(J:N, J) has been initialized to be A(J, J:N)
+*        H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
+*         where H(J:M, J) has been initialized to be A(J, J:M)
 *
          IF( K.GT.2 ) THEN
 *
@@ -234,14 +234,14 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(i:n, i) into WORK
+*        Copy H(i:M, i) into WORK
 *
          CALL CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
-*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
+*           Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
+*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
 *
             ALPHA = -A( K-1, J )
             CALL CAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
@@ -253,8 +253,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
-*            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
+*           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
+*            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( K, J )
@@ -262,7 +262,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -277,14 +277,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1, I1+1:N) with A(I1+1:N, I2)
+*              Swap A(I1, I1+1:M) with A(I1+1:M, I2)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL CSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
      $                              A( J1+I1, I2 ), 1 )
 *
-*              Swap A(I1, I2+1:N) with A(I2, I2+1:N)
+*              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
 *
                CALL CSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
      $                           A( J1+I2-1, I2+1 ), LDA )
@@ -324,14 +324,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J:N, J),
+*              Copy A(J+1:M, J+1) into H(J:M, J),
 *
                CALL CCOPY( M-J, A( K+1, J+1 ), LDA,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( K, J+1 ).NE.ZERO ) THEN
                ALPHA = ONE / A( K, J+1 )
@@ -367,8 +367,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
-*         where H(J:N, J) has been initialized to be A(J:N, J)
+*        H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
+*         where H(J:M, J) has been initialized to be A(J:M, J)
 *
          IF( K.GT.2 ) THEN
 *
@@ -384,13 +384,13 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(J:N, J) into WORK
+*        Copy H(J:M, J) into WORK
 *
          CALL CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
+*           Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 *
             ALPHA = -A( J, K-1 )
@@ -403,8 +403,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L((J+1):N, J)
-*            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
+*           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
+*            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( J, K )
@@ -412,7 +412,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -427,14 +427,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1+1:N, I1) with A(I2, I1+1:N)
+*              Swap A(I1+1:M, I1) with A(I2, I1+1:M)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL CSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
      $                              A( I2, J1+I1 ), LDA )
 *
-*              Swap A(I2+1:N, I1) with A(I2+1:N, I2)
+*              Swap A(I2+1:M, I1) with A(I2+1:M, I2)
 *
                CALL CSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
      $                           A( I2+1, J1+I2-1 ), 1 )
@@ -473,14 +473,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J+1:N, J),
+*              Copy A(J+1:M, J+1) into H(J+1:M, J),
 *
                CALL CCOPY( M-J, A( J+1, K+1 ), 1,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( J+1, K ).NE.ZERO ) THEN
                ALPHA = ONE / A( J+1, K )
diff --git a/SRC/dlasyf_aa.f b/SRC/dlasyf_aa.f
index 0bd2d6de..d9550e11 100644
--- a/SRC/dlasyf_aa.f
+++ b/SRC/dlasyf_aa.f
@@ -99,12 +99,12 @@
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.  LDA >= max(1,N).
+*>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
 *> \param[out] IPIV
 *> \verbatim
-*>          IPIV is INTEGER array, dimension (N)
+*>          IPIV is INTEGER array, dimension (M)
 *>          Details of the row and column interchanges,
 *>          the row and column k were interchanged with the row and
 *>          column IPIV(k).
@@ -217,8 +217,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
-*         where H(J:N, J) has been initialized to be A(J, J:N)
+*        H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
+*         where H(J:M, J) has been initialized to be A(J, J:M)
 *
          IF( K.GT.2 ) THEN
 *
@@ -234,14 +234,14 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(i:n, i) into WORK
+*        Copy H(i:M, i) into WORK
 *
          CALL DCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
-*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
+*           Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
+*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
 *
             ALPHA = -A( K-1, J )
             CALL DAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
@@ -253,8 +253,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
-*            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
+*           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
+*            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( K, J )
@@ -262,7 +262,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -277,14 +277,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1, I1+1:N) with A(I1+1:N, I2)
+*              Swap A(I1, I1+1:M) with A(I1+1:M, I2)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
      $                              A( J1+I1, I2 ), 1 )
 *
-*              Swap A(I1, I2+1:N) with A(I2, I2+1:N)
+*              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
 *
                CALL DSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
      $                           A( J1+I2-1, I2+1 ), LDA )
@@ -324,14 +324,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J:N, J),
+*              Copy A(J+1:M, J+1) into H(J:M, J),
 *
                CALL DCOPY( M-J, A( K+1, J+1 ), LDA,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( K, J+1 ).NE.ZERO ) THEN
                ALPHA = ONE / A( K, J+1 )
@@ -367,8 +367,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
-*         where H(J:N, J) has been initialized to be A(J:N, J)
+*        H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
+*         where H(J:M, J) has been initialized to be A(J:M, J)
 *
          IF( K.GT.2 ) THEN
 *
@@ -384,13 +384,13 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(J:N, J) into WORK
+*        Copy H(J:M, J) into WORK
 *
          CALL DCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
+*           Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 *
             ALPHA = -A( J, K-1 )
@@ -403,8 +403,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L((J+1):N, J)
-*            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
+*           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
+*            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( J, K )
@@ -412,7 +412,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -427,14 +427,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1+1:N, I1) with A(I2, I1+1:N)
+*              Swap A(I1+1:M, I1) with A(I2, I1+1:M)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
      $                              A( I2, J1+I1 ), LDA )
 *
-*              Swap A(I2+1:N, I1) with A(I2+1:N, I2)
+*              Swap A(I2+1:M, I1) with A(I2+1:M, I2)
 *
                CALL DSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
      $                           A( I2+1, J1+I2-1 ), 1 )
@@ -473,14 +473,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J+1:N, J),
+*              Copy A(J+1:M, J+1) into H(J+1:M, J),
 *
                CALL DCOPY( M-J, A( J+1, K+1 ), 1,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( J+1, K ).NE.ZERO ) THEN
                ALPHA = ONE / A( J+1, K )
diff --git a/SRC/dsb2st_kernels.f b/SRC/dsb2st_kernels.f
index afed5265..988ffc8d 100644
--- a/SRC/dsb2st_kernels.f
+++ b/SRC/dsb2st_kernels.f
@@ -47,45 +47,87 @@
 *  Arguments:
 *  ==========
 *
-*> @param[in] n
-*>          The order of the matrix A.
-*>
-*> @param[in] nb
-*>          The size of the band.
-*>
-*> @param[in, out] A
-*>          A pointer to the matrix A.
-*>
-*> @param[in] lda
-*>          The leading dimension of the matrix A.
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*> \endverbatim
 *>
-*> @param[out] V
-*>          DOUBLE PRECISION array, dimension 2*n if eigenvalues only are
-*>          requested or to be queried for vectors.
+*> \param[in] WANTZ
+*> \verbatim
+*>          WANTZ is LOGICAL which indicate if Eigenvalue are requested or both
+*>          Eigenvalue/Eigenvectors.
+*> \endverbatim
 *>
-*> @param[out] TAU
-*>          DOUBLE PRECISION array, dimension (2*n).
-*>          The scalar factors of the Householder reflectors are stored
-*>          in this array.
+*> \param[in] TTYPE
+*> \verbatim
+*>          TTYPE is INTEGER
+*> \endverbatim
 *>
-*> @param[in] st
+*> \param[in] ST
+*> \verbatim
+*>          ST is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] ed
+*> \param[in] ED
+*> \verbatim
+*>          ED is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] sweep
+*> \param[in] SWEEP
+*> \verbatim
+*>          SWEEP is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] Vblksiz
-*>          internal parameter for indices.
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER. The order of the matrix A.
+*> \endverbatim
 *>
-*> @param[in] wantz
-*>          logical which indicate if Eigenvalue are requested or both
-*>          Eigenvalue/Eigenvectors.
+*> \param[in] NB
+*> \verbatim
+*>          NB is INTEGER. The size of the band.
+*> \endverbatim
+*>
+*> \param[in] IB
+*> \verbatim
+*>          IB is INTEGER.
+*> \endverbatim
+*>
+*> \param[in, out] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array. A pointer to the matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER. The leading dimension of the matrix A.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*>          V is DOUBLE PRECISION array, dimension 2*n if eigenvalues only are
+*>          requested or to be queried for vectors.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is DOUBLE PRECISION array, dimension (2*n).
+*>          The scalar factors of the Householder reflectors are stored
+*>          in this array.
+*> \endverbatim
+*>
+*> \param[in] LDVT
+*> \verbatim
+*>          LDVT is INTEGER.
+*> \endverbatim
 *>
-*> @param[in] work
-*>          Workspace of size nb.
+*> \param[in] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION array. Workspace of size nb.
+*> \endverbatim
 *>
 *> \par Further Details:
 *  =====================
diff --git a/SRC/slasyf_aa.f b/SRC/slasyf_aa.f
index 5fb3cc9a..3b4bc4b1 100644
--- a/SRC/slasyf_aa.f
+++ b/SRC/slasyf_aa.f
@@ -99,12 +99,12 @@
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.  LDA >= max(1,N).
+*>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
 *> \param[out] IPIV
 *> \verbatim
-*>          IPIV is INTEGER array, dimension (N)
+*>          IPIV is INTEGER array, dimension (M)
 *>          Details of the row and column interchanges,
 *>          the row and column k were interchanged with the row and
 *>          column IPIV(k).
@@ -217,8 +217,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
-*         where H(J:N, J) has been initialized to be A(J, J:N)
+*        H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
+*         where H(J:M, J) has been initialized to be A(J, J:M)
 *
          IF( K.GT.2 ) THEN
 *
@@ -234,14 +234,14 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(i:n, i) into WORK
+*        Copy H(i:M, i) into WORK
 *
          CALL SCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
-*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
+*           Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
+*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
 *
             ALPHA = -A( K-1, J )
             CALL SAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
@@ -253,8 +253,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
-*            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
+*           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
+*            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( K, J )
@@ -262,7 +262,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = ISAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -277,14 +277,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1, I1+1:N) with A(I1+1:N, I2)
+*              Swap A(I1, I1+1:M) with A(I1+1:M, I2)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL SSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
      $                              A( J1+I1, I2 ), 1 )
 *
-*              Swap A(I1, I2+1:N) with A(I2, I2+1:N)
+*              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
 *
                CALL SSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
      $                           A( J1+I2-1, I2+1 ), LDA )
@@ -324,14 +324,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J:N, J),
+*              Copy A(J+1:M, J+1) into H(J:M, J),
 *
                CALL SCOPY( M-J, A( K+1, J+1 ), LDA,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( K, J+1 ).NE.ZERO ) THEN
                ALPHA = ONE / A( K, J+1 )
@@ -367,8 +367,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
-*         where H(J:N, J) has been initialized to be A(J:N, J)
+*        H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
+*         where H(J:M, J) has been initialized to be A(J:M, J)
 *
          IF( K.GT.2 ) THEN
 *
@@ -384,13 +384,13 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(J:N, J) into WORK
+*        Copy H(J:M, J) into WORK
 *
          CALL SCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
+*           Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 *
             ALPHA = -A( J, K-1 )
@@ -403,8 +403,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L((J+1):N, J)
-*            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
+*           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
+*            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( J, K )
@@ -412,7 +412,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = ISAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -427,14 +427,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1+1:N, I1) with A(I2, I1+1:N)
+*              Swap A(I1+1:M, I1) with A(I2, I1+1:M)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL SSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
      $                              A( I2, J1+I1 ), LDA )
 *
-*              Swap A(I2+1:N, I1) with A(I2+1:N, I2)
+*              Swap A(I2+1:M, I1) with A(I2+1:M, I2)
 *
                CALL SSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
      $                           A( I2+1, J1+I2-1 ), 1 )
@@ -473,14 +473,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J+1:N, J),
+*              Copy A(J+1:M, J+1) into H(J+1:M, J),
 *
                CALL SCOPY( M-J, A( J+1, K+1 ), 1,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( J+1, K ).NE.ZERO ) THEN
                ALPHA = ONE / A( J+1, K )
diff --git a/SRC/ssb2st_kernels.f b/SRC/ssb2st_kernels.f
index 7183c9ad..4d046d8f 100644
--- a/SRC/ssb2st_kernels.f
+++ b/SRC/ssb2st_kernels.f
@@ -47,45 +47,90 @@
 *  Arguments:
 *  ==========
 *
-*> @param[in] n
-*>          The order of the matrix A.
-*>
-*> @param[in] nb
-*>          The size of the band.
-*>
-*> @param[in, out] A
-*>          A pointer to the matrix A.
-*>
-*> @param[in] lda
-*>          The leading dimension of the matrix A.
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*> \endverbatim
 *>
-*> @param[out] V
-*>          REAL array, dimension 2*n if eigenvalues only are
-*>          requested or to be queried for vectors.
+*> \param[in] WANTZ
+*> \verbatim
+*>          WANTZ is LOGICAL which indicate if Eigenvalue are requested or both
+*>          Eigenvalue/Eigenvectors.
+*> \endverbatim
 *>
-*> @param[out] TAU
-*>          REAL array, dimension (2*n).
-*>          The scalar factors of the Householder reflectors are stored
-*>          in this array.
+*> \param[in] TTYPE
+*> \verbatim
+*>          TTYPE is INTEGER
+*> \endverbatim
 *>
-*> @param[in] st
+*> \param[in] ST
+*> \verbatim
+*>          ST is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] ed
+*> \param[in] ED
+*> \verbatim
+*>          ED is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] sweep
+*> \param[in] SWEEP
+*> \verbatim
+*>          SWEEP is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] Vblksiz
-*>          internal parameter for indices.
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER. The order of the matrix A.
+*> \endverbatim
 *>
-*> @param[in] wantz
-*>          logical which indicate if Eigenvalue are requested or both
-*>          Eigenvalue/Eigenvectors.
+*> \param[in] NB
+*> \verbatim
+*>          NB is INTEGER. The size of the band.
+*> \endverbatim
+*>
+*> \param[in] IB
+*> \verbatim
+*>          IB is INTEGER.
+*> \endverbatim
+*>
+*> \param[in, out] A
+*> \verbatim
+*>          A is REAL array. A pointer to the matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER. The leading dimension of the matrix A.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*>          V is REAL array, dimension 2*n if eigenvalues only are
+*>          requested or to be queried for vectors.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is REAL array, dimension (2*n).
+*>          The scalar factors of the Householder reflectors are stored
+*>          in this array.
+*> \endverbatim
+*>
+*> \param[in] LDVT
+*> \verbatim
+*>          LDVT is INTEGER.
+*> \endverbatim
+*>
+*> \param[in] WORK
+*> \verbatim
+*>          WORK is REAL array. Workspace of size nb.
+*> \endverbatim
+*> @param[in] n
+*>          The order of the matrix A.
 *>
-*> @param[in] work
-*>          Workspace of size nb.
 *>
 *> \par Further Details:
 *  =====================
diff --git a/SRC/zhb2st_kernels.f b/SRC/zhb2st_kernels.f
index e4114b5f..862eb9d5 100644
--- a/SRC/zhb2st_kernels.f
+++ b/SRC/zhb2st_kernels.f
@@ -47,45 +47,87 @@
 *  Arguments:
 *  ==========
 *
-*> @param[in] n
-*>          The order of the matrix A.
-*>
-*> @param[in] nb
-*>          The size of the band.
-*>
-*> @param[in, out] A
-*>          A pointer to the matrix A.
-*>
-*> @param[in] lda
-*>          The leading dimension of the matrix A.
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*> \endverbatim
 *>
-*> @param[out] V
-*>          COMPLEX*16 array, dimension 2*n if eigenvalues only are
-*>          requested or to be queried for vectors.
+*> \param[in] WANTZ
+*> \verbatim
+*>          WANTZ is LOGICAL which indicate if Eigenvalue are requested or both
+*>          Eigenvalue/Eigenvectors.
+*> \endverbatim
 *>
-*> @param[out] TAU
-*>          COMPLEX*16 array, dimension (2*n).
-*>          The scalar factors of the Householder reflectors are stored
-*>          in this array.
+*> \param[in] TTYPE
+*> \verbatim
+*>          TTYPE is INTEGER
+*> \endverbatim
 *>
-*> @param[in] st
+*> \param[in] ST
+*> \verbatim
+*>          ST is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] ed
+*> \param[in] ED
+*> \verbatim
+*>          ED is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] sweep
+*> \param[in] SWEEP
+*> \verbatim
+*>          SWEEP is INTEGER
 *>          internal parameter for indices.
+*> \endverbatim
 *>
-*> @param[in] Vblksiz
-*>          internal parameter for indices.
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER. The order of the matrix A.
+*> \endverbatim
 *>
-*> @param[in] wantz
-*>          logical which indicate if Eigenvalue are requested or both
-*>          Eigenvalue/Eigenvectors.
+*> \param[in] NB
+*> \verbatim
+*>          NB is INTEGER. The size of the band.
+*> \endverbatim
+*>
+*> \param[in] IB
+*> \verbatim
+*>          IB is INTEGER.
+*> \endverbatim
+*>
+*> \param[in, out] A
+*> \verbatim
+*>          A is COMPLEX*16 array. A pointer to the matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER. The leading dimension of the matrix A.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*>          V is COMPLEX*16 array, dimension 2*n if eigenvalues only are
+*>          requested or to be queried for vectors.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is COMPLEX*16 array, dimension (2*n).
+*>          The scalar factors of the Householder reflectors are stored
+*>          in this array.
+*> \endverbatim
+*>
+*> \param[in] LDVT
+*> \verbatim
+*>          LDVT is INTEGER.
+*> \endverbatim
 *>
-*> @param[in] work
-*>          Workspace of size nb.
+*> \param[in] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array. Workspace of size nb.
+*> \endverbatim
 *>
 *> \par Further Details:
 *  =====================
diff --git a/SRC/zlasyf_aa.f b/SRC/zlasyf_aa.f
index 7ac4ff3d..dff71042 100644
--- a/SRC/zlasyf_aa.f
+++ b/SRC/zlasyf_aa.f
@@ -99,12 +99,12 @@
 *> \param[in] LDA
 *> \verbatim
 *>          LDA is INTEGER
-*>          The leading dimension of the array A.  LDA >= max(1,N).
+*>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim
 *>
 *> \param[out] IPIV
 *> \verbatim
-*>          IPIV is INTEGER array, dimension (N)
+*>          IPIV is INTEGER array, dimension (M)
 *>          Details of the row and column interchanges,
 *>          the row and column k were interchanged with the row and
 *>          column IPIV(k).
@@ -217,8 +217,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
-*         where H(J:N, J) has been initialized to be A(J, J:N)
+*        H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
+*         where H(J:M, J) has been initialized to be A(J, J:M)
 *
          IF( K.GT.2 ) THEN
 *
@@ -234,14 +234,14 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(i:n, i) into WORK
+*        Copy H(i:M, i) into WORK
 *
          CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
-*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
+*           Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
+*            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
 *
             ALPHA = -A( K-1, J )
             CALL ZAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
@@ -253,8 +253,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
-*            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
+*           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
+*            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( K, J )
@@ -262,7 +262,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -277,14 +277,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1, I1+1:N) with A(I1+1:N, I2)
+*              Swap A(I1, I1+1:M) with A(I1+1:M, I2)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL ZSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
      $                              A( J1+I1, I2 ), 1 )
 *
-*              Swap A(I1, I2+1:N) with A(I2, I2+1:N)
+*              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
 *
                CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
      $                           A( J1+I2-1, I2+1 ), LDA )
@@ -324,14 +324,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J:N, J),
+*              Copy A(J+1:M, J+1) into H(J:M, J),
 *
                CALL ZCOPY( M-J, A( K+1, J+1 ), LDA,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( K, J+1 ).NE.ZERO ) THEN
                ALPHA = ONE / A( K, J+1 )
@@ -367,8 +367,8 @@
 *
          K = J1+J-1
 *
-*        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
-*         where H(J:N, J) has been initialized to be A(J:N, J)
+*        H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
+*         where H(J:M, J) has been initialized to be A(J:M, J)
 *
          IF( K.GT.2 ) THEN
 *
@@ -384,13 +384,13 @@
      $                  ONE, H( J, J ), 1 )
          END IF
 *
-*        Copy H(J:N, J) into WORK
+*        Copy H(J:M, J) into WORK
 *
          CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 *
          IF( J.GT.K1 ) THEN
 *
-*           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
+*           Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 *
             ALPHA = -A( J, K-1 )
@@ -403,8 +403,8 @@
 *
          IF( J.LT.M ) THEN
 *
-*           Compute WORK(2:N) = T(J, J) L((J+1):N, J)
-*            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
+*           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
+*            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
 *
             IF( K.GT.1 ) THEN
                ALPHA = -A( J, K )
@@ -412,7 +412,7 @@
      $                                 WORK( 2 ), 1 )
             ENDIF
 *
-*           Find max(|WORK(2:n)|)
+*           Find max(|WORK(2:M)|)
 *
             I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )
@@ -427,14 +427,14 @@
                WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV
 *
-*              Swap A(I1+1:N, I1) with A(I2, I1+1:N)
+*              Swap A(I1+1:M, I1) with A(I2, I1+1:M)
 *
                I1 = I1+J-1
                I2 = I2+J-1
                CALL ZSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
      $                              A( I2, J1+I1 ), LDA )
 *
-*              Swap A(I2+1:N, I1) with A(I2+1:N, I2)
+*              Swap A(I2+1:M, I1) with A(I2+1:M, I2)
 *
                CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
      $                           A( I2+1, J1+I2-1 ), 1 )
@@ -473,14 +473,14 @@
 *
             IF( J.LT.NB ) THEN
 *
-*              Copy A(J+1:N, J+1) into H(J+1:N, J),
+*              Copy A(J+1:M, J+1) into H(J+1:M, J),
 *
                CALL ZCOPY( M-J, A( J+1, K+1 ), 1,
      $                          H( J+1, J+1 ), 1 )
             END IF
 *
-*           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
-*            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
+*           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
+*            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *
             IF( A( J+1, K ).NE.ZERO ) THEN
                ALPHA = ONE / A( J+1, K )