implement GBMV

6 files changed, 64 insertions, 1255 deletions

diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt
index c8d5c241f..5cbf90ddd 100644
--- a/blas/CMakeLists.txt
+++ b/blas/CMakeLists.txt
@@ -19,7 +19,7 @@ add_custom_target(blas)
 
 set(EigenBlas_SRCS single.cpp double.cpp complex_single.cpp complex_double.cpp xerbla.cpp
     srotm.f srotmg.f drotm.f drotmg.f
-    lsame.f cgbmv.f  chpr2.f  ctbsv.f  dspmv.f  dtbmv.f  dtpsv.f ssbmv.f  sspr.f   stpmv.f   zgbmv.f  zhpr2.f  ztbsv.f chbmv.f  chpr.f   ctpmv.f  dgbmv.f   dspr2.f dtbsv.f  sspmv.f  stbmv.f  stpsv.f   zhbmv.f  zhpr.f   ztpmv.f chpmv.f  ctbmv.f  ctpsv.f    dsbmv.f  dspr.f   dtpmv.f  sgbmv.f sspr2.f  stbsv.f  zhpmv.f  ztbmv.f  ztpsv.f
+    lsame.f   chpr2.f  ctbsv.f  dspmv.f  dtbmv.f  dtpsv.f ssbmv.f  sspr.f   stpmv.f     zhpr2.f  ztbsv.f chbmv.f  chpr.f   ctpmv.f     dspr2.f dtbsv.f  sspmv.f  stbmv.f  stpsv.f   zhbmv.f  zhpr.f   ztpmv.f chpmv.f  ctbmv.f  ctpsv.f    dsbmv.f  dspr.f   dtpmv.f   sspr2.f  stbsv.f  zhpmv.f  ztbmv.f  ztpsv.f
 )
 
 add_library(eigen_blas_static ${EigenBlas_SRCS})
diff --git a/blas/cgbmv.f b/blas/cgbmv.f
deleted file mode 100644
index 2a837dba3..000000000
--- a/blas/cgbmv.f
+++ /dev/null
@@ -1,322 +0,0 @@
-      SUBROUTINE CGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-*     .. Scalar Arguments ..
-      COMPLEX ALPHA,BETA
-      INTEGER INCX,INCY,KL,KU,LDA,M,N
-      CHARACTER TRANS
-*     ..
-*     .. Array Arguments ..
-      COMPLEX A(LDA,*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  CGBMV  performs one of the matrix-vector operations
-*
-*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or
-*
-*     y := alpha*conjg( A' )*x + beta*y,
-*
-*  where alpha and beta are scalars, x and y are vectors and A is an
-*  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
-*
-*  Arguments
-*  ==========
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the operation to be performed as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
-*
-*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
-*
-*              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y.
-*
-*           Unchanged on exit.
-*
-*  M      - INTEGER.
-*           On entry, M specifies the number of rows of the matrix A.
-*           M must be at least zero.
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the number of columns of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  KL     - INTEGER.
-*           On entry, KL specifies the number of sub-diagonals of the
-*           matrix A. KL must satisfy  0 .le. KL.
-*           Unchanged on exit.
-*
-*  KU     - INTEGER.
-*           On entry, KU specifies the number of super-diagonals of the
-*           matrix A. KU must satisfy  0 .le. KU.
-*           Unchanged on exit.
-*
-*  ALPHA  - COMPLEX         .
-*           On entry, ALPHA specifies the scalar alpha.
-*           Unchanged on exit.
-*
-*  A      - COMPLEX          array of DIMENSION ( LDA, n ).
-*           Before entry, the leading ( kl + ku + 1 ) by n part of the
-*           array A must contain the matrix of coefficients, supplied
-*           column by column, with the leading diagonal of the matrix in
-*           row ( ku + 1 ) of the array, the first super-diagonal
-*           starting at position 2 in row ku, the first sub-diagonal
-*           starting at position 1 in row ( ku + 2 ), and so on.
-*           Elements in the array A that do not correspond to elements
-*           in the band matrix (such as the top left ku by ku triangle)
-*           are not referenced.
-*           The following program segment will transfer a band matrix
-*           from conventional full matrix storage to band storage:
-*
-*                 DO 20, J = 1, N
-*                    K = KU + 1 - J
-*                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
-*                       A( K + I, J ) = matrix( I, J )
-*              10    CONTINUE
-*              20 CONTINUE
-*
-*           Unchanged on exit.
-*
-*  LDA    - INTEGER.
-*           On entry, LDA specifies the first dimension of A as declared
-*           in the calling (sub) program. LDA must be at least
-*           ( kl + ku + 1 ).
-*           Unchanged on exit.
-*
-*  X      - COMPLEX          array of DIMENSION at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
-*           Before entry, the incremented array X must contain the
-*           vector x.
-*           Unchanged on exit.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  BETA   - COMPLEX         .
-*           On entry, BETA specifies the scalar beta. When BETA is
-*           supplied as zero then Y need not be set on input.
-*           Unchanged on exit.
-*
-*  Y      - COMPLEX          array of DIMENSION at least
-*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
-*           Before entry, the incremented array Y must contain the
-*           vector y. On exit, Y is overwritten by the updated vector y.
-*
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ONE
-      PARAMETER (ONE= (1.0E+0,0.0E+0))
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
-*     ..
-*     .. Local Scalars ..
-      COMPLEX TEMP
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
-      LOGICAL NOCONJ
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC CONJG,MAX,MIN
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +    .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 1
-      ELSE IF (M.LT.0) THEN
-          INFO = 2
-      ELSE IF (N.LT.0) THEN
-          INFO = 3
-      ELSE IF (KL.LT.0) THEN
-          INFO = 4
-      ELSE IF (KU.LT.0) THEN
-          INFO = 5
-      ELSE IF (LDA.LT. (KL+KU+1)) THEN
-          INFO = 8
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 10
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 13
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('CGBMV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
-     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
-*
-      NOCONJ = LSAME(TRANS,'T')
-*
-*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
-*     up the start points in  X  and  Y.
-*
-      IF (LSAME(TRANS,'N')) THEN
-          LENX = N
-          LENY = M
-      ELSE
-          LENX = M
-          LENY = N
-      END IF
-      IF (INCX.GT.0) THEN
-          KX = 1
-      ELSE
-          KX = 1 - (LENX-1)*INCX
-      END IF
-      IF (INCY.GT.0) THEN
-          KY = 1
-      ELSE
-          KY = 1 - (LENY-1)*INCY
-      END IF
-*
-*     Start the operations. In this version the elements of A are
-*     accessed sequentially with one pass through the band part of A.
-*
-*     First form  y := beta*y.
-*
-      IF (BETA.NE.ONE) THEN
-          IF (INCY.EQ.1) THEN
-              IF (BETA.EQ.ZERO) THEN
-                  DO 10 I = 1,LENY
-                      Y(I) = ZERO
-   10             CONTINUE
-              ELSE
-                  DO 20 I = 1,LENY
-                      Y(I) = BETA*Y(I)
-   20             CONTINUE
-              END IF
-          ELSE
-              IY = KY
-              IF (BETA.EQ.ZERO) THEN
-                  DO 30 I = 1,LENY
-                      Y(IY) = ZERO
-                      IY = IY + INCY
-   30             CONTINUE
-              ELSE
-                  DO 40 I = 1,LENY
-                      Y(IY) = BETA*Y(IY)
-                      IY = IY + INCY
-   40             CONTINUE
-              END IF
-          END IF
-      END IF
-      IF (ALPHA.EQ.ZERO) RETURN
-      KUP1 = KU + 1
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  y := alpha*A*x + y.
-*
-          JX = KX
-          IF (INCY.EQ.1) THEN
-              DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  IF (J.GT.KU) KY = KY + INCY
-   80         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y.
-*
-          JY = KY
-          IF (INCX.EQ.1) THEN
-              DO 110 J = 1,N
-                  TEMP = ZERO
-                  K = KUP1 - J
-                  IF (NOCONJ) THEN
-                      DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + A(K+I,J)*X(I)
-   90                 CONTINUE
-                  ELSE
-                      DO 100 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + CONJG(A(K+I,J))*X(I)
-  100                 CONTINUE
-                  END IF
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-  110         CONTINUE
-          ELSE
-              DO 140 J = 1,N
-                  TEMP = ZERO
-                  IX = KX
-                  K = KUP1 - J
-                  IF (NOCONJ) THEN
-                      DO 120 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + A(K+I,J)*X(IX)
-                          IX = IX + INCX
-  120                 CONTINUE
-                  ELSE
-                      DO 130 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + CONJG(A(K+I,J))*X(IX)
-                          IX = IX + INCX
-  130                 CONTINUE
-                  END IF
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-                  IF (J.GT.KU) KX = KX + INCX
-  140         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of CGBMV .
-*
-      END
diff --git a/blas/dgbmv.f b/blas/dgbmv.f
deleted file mode 100644
index c3dc64aa3..000000000
--- a/blas/dgbmv.f
+++ /dev/null
@@ -1,301 +0,0 @@
-      SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-*     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA,BETA
-      INTEGER INCX,INCY,KL,KU,LDA,M,N
-      CHARACTER TRANS
-*     ..
-*     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DGBMV  performs one of the matrix-vector operations
-*
-*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
-*
-*  where alpha and beta are scalars, x and y are vectors and A is an
-*  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
-*
-*  Arguments
-*  ==========
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the operation to be performed as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
-*
-*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
-*
-*              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
-*
-*           Unchanged on exit.
-*
-*  M      - INTEGER.
-*           On entry, M specifies the number of rows of the matrix A.
-*           M must be at least zero.
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the number of columns of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  KL     - INTEGER.
-*           On entry, KL specifies the number of sub-diagonals of the
-*           matrix A. KL must satisfy  0 .le. KL.
-*           Unchanged on exit.
-*
-*  KU     - INTEGER.
-*           On entry, KU specifies the number of super-diagonals of the
-*           matrix A. KU must satisfy  0 .le. KU.
-*           Unchanged on exit.
-*
-*  ALPHA  - DOUBLE PRECISION.
-*           On entry, ALPHA specifies the scalar alpha.
-*           Unchanged on exit.
-*
-*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
-*           Before entry, the leading ( kl + ku + 1 ) by n part of the
-*           array A must contain the matrix of coefficients, supplied
-*           column by column, with the leading diagonal of the matrix in
-*           row ( ku + 1 ) of the array, the first super-diagonal
-*           starting at position 2 in row ku, the first sub-diagonal
-*           starting at position 1 in row ( ku + 2 ), and so on.
-*           Elements in the array A that do not correspond to elements
-*           in the band matrix (such as the top left ku by ku triangle)
-*           are not referenced.
-*           The following program segment will transfer a band matrix
-*           from conventional full matrix storage to band storage:
-*
-*                 DO 20, J = 1, N
-*                    K = KU + 1 - J
-*                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
-*                       A( K + I, J ) = matrix( I, J )
-*              10    CONTINUE
-*              20 CONTINUE
-*
-*           Unchanged on exit.
-*
-*  LDA    - INTEGER.
-*           On entry, LDA specifies the first dimension of A as declared
-*           in the calling (sub) program. LDA must be at least
-*           ( kl + ku + 1 ).
-*           Unchanged on exit.
-*
-*  X      - DOUBLE PRECISION array of DIMENSION at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
-*           Before entry, the incremented array X must contain the
-*           vector x.
-*           Unchanged on exit.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  BETA   - DOUBLE PRECISION.
-*           On entry, BETA specifies the scalar beta. When BETA is
-*           supplied as zero then Y need not be set on input.
-*           Unchanged on exit.
-*
-*  Y      - DOUBLE PRECISION array of DIMENSION at least
-*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
-*           Before entry, the incremented array Y must contain the
-*           vector y. On exit, Y is overwritten by the updated vector y.
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
-      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
-*     ..
-*     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC MAX,MIN
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +    .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 1
-      ELSE IF (M.LT.0) THEN
-          INFO = 2
-      ELSE IF (N.LT.0) THEN
-          INFO = 3
-      ELSE IF (KL.LT.0) THEN
-          INFO = 4
-      ELSE IF (KU.LT.0) THEN
-          INFO = 5
-      ELSE IF (LDA.LT. (KL+KU+1)) THEN
-          INFO = 8
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 10
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 13
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('DGBMV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
-     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
-*
-*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
-*     up the start points in  X  and  Y.
-*
-      IF (LSAME(TRANS,'N')) THEN
-          LENX = N
-          LENY = M
-      ELSE
-          LENX = M
-          LENY = N
-      END IF
-      IF (INCX.GT.0) THEN
-          KX = 1
-      ELSE
-          KX = 1 - (LENX-1)*INCX
-      END IF
-      IF (INCY.GT.0) THEN
-          KY = 1
-      ELSE
-          KY = 1 - (LENY-1)*INCY
-      END IF
-*
-*     Start the operations. In this version the elements of A are
-*     accessed sequentially with one pass through the band part of A.
-*
-*     First form  y := beta*y.
-*
-      IF (BETA.NE.ONE) THEN
-          IF (INCY.EQ.1) THEN
-              IF (BETA.EQ.ZERO) THEN
-                  DO 10 I = 1,LENY
-                      Y(I) = ZERO
-   10             CONTINUE
-              ELSE
-                  DO 20 I = 1,LENY
-                      Y(I) = BETA*Y(I)
-   20             CONTINUE
-              END IF
-          ELSE
-              IY = KY
-              IF (BETA.EQ.ZERO) THEN
-                  DO 30 I = 1,LENY
-                      Y(IY) = ZERO
-                      IY = IY + INCY
-   30             CONTINUE
-              ELSE
-                  DO 40 I = 1,LENY
-                      Y(IY) = BETA*Y(IY)
-                      IY = IY + INCY
-   40             CONTINUE
-              END IF
-          END IF
-      END IF
-      IF (ALPHA.EQ.ZERO) RETURN
-      KUP1 = KU + 1
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  y := alpha*A*x + y.
-*
-          JX = KX
-          IF (INCY.EQ.1) THEN
-              DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  IF (J.GT.KU) KY = KY + INCY
-   80         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  y := alpha*A'*x + y.
-*
-          JY = KY
-          IF (INCX.EQ.1) THEN
-              DO 100 J = 1,N
-                  TEMP = ZERO
-                  K = KUP1 - J
-                  DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
-                      TEMP = TEMP + A(K+I,J)*X(I)
-   90             CONTINUE
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-  100         CONTINUE
-          ELSE
-              DO 120 J = 1,N
-                  TEMP = ZERO
-                  IX = KX
-                  K = KUP1 - J
-                  DO 110 I = MAX(1,J-KU),MIN(M,J+KL)
-                      TEMP = TEMP + A(K+I,J)*X(IX)
-                      IX = IX + INCX
-  110             CONTINUE
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-                  IF (J.GT.KU) KX = KX + INCX
-  120         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of DGBMV .
-*
-      END
diff --git a/blas/level2_impl.h b/blas/level2_impl.h
index 0e4649f57..8cbc2f424 100644
--- a/blas/level2_impl.h
+++ b/blas/level2_impl.h
@@ -202,21 +202,76 @@ int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar
   return 0;
 }
 
-/**  DGBMV  performs one of the matrix-vector operations
+/**  GBMV  performs one of the matrix-vector operations
   *
   *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
   *
   *  where alpha and beta are scalars, x and y are vectors and A is an
   *  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
   */
-// int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *alpha, RealScalar *a, int *lda,
-//                           RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
-// {
-//   return 1;
-// }
-
+int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda,
+                          RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
+{
+  Scalar* a = reinterpret_cast<Scalar*>(pa);
+  Scalar* x = reinterpret_cast<Scalar*>(px);
+  Scalar* y = reinterpret_cast<Scalar*>(py);
+  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
+  Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
+  int coeff_rows = *kl+*ku+1;
+  
+  int info = 0;
+       if(OP(*trans)==INVALID)                                        info = 1;
+  else if(*m<0)                                                       info = 2;
+  else if(*n<0)                                                       info = 3;
+  else if(*kl<0)                                                      info = 4;
+  else if(*ku<0)                                                      info = 5;
+  else if(*lda<coeff_rows)                                            info = 8;
+  else if(*incx==0)                                                   info = 10;
+  else if(*incy==0)                                                   info = 13;
+  if(info)
+    return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6);
+  
+  if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
+    return 0;
+  
+  int actual_m = *m;
+  int actual_n = *n;
+  if(OP(*trans)!=NOTR)
+    std::swap(actual_m,actual_n);
+  
+  Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
+  Scalar* actual_y = get_compact_vector(y,actual_m,*incy);
+  
+  if(beta!=Scalar(1))
+  {
+    if(beta==Scalar(0)) vector(actual_y, actual_m).setZero();
+    else                vector(actual_y, actual_m) *= beta;
+  }
+  
+  MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
+  
+  int nb = std::min(*n,(*m)+(*ku));
+  for(int j=0; j<nb; ++j)
+  {
+    int start = std::max(0,j - *ku);
+    int end = std::min((*m)-1,j + *kl);
+    int len = end - start + 1;
+    int offset = (*ku) - j + start;
+    if(OP(*trans)==NOTR)
+      vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
+    else if(OP(*trans)==TR)
+      actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value();
+    else
+      actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint()   * vector(actual_x+start,len) ).value();
+  }    
+  
+  if(actual_x!=x) delete[] actual_x;
+  if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
+  
+  return 0;
+}
 
-/**  DTBMV  performs one of the matrix-vector operations
+/**  TBMV  performs one of the matrix-vector operations
   *
   *     x := A*x,   or   x := A'*x,
   *
diff --git a/blas/sgbmv.f b/blas/sgbmv.f
deleted file mode 100644
index aaa8b1a17..000000000
--- a/blas/sgbmv.f
+++ /dev/null
@@ -1,301 +0,0 @@
-      SUBROUTINE SGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-*     .. Scalar Arguments ..
-      REAL ALPHA,BETA
-      INTEGER INCX,INCY,KL,KU,LDA,M,N
-      CHARACTER TRANS
-*     ..
-*     .. Array Arguments ..
-      REAL A(LDA,*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  SGBMV  performs one of the matrix-vector operations
-*
-*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
-*
-*  where alpha and beta are scalars, x and y are vectors and A is an
-*  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
-*
-*  Arguments
-*  ==========
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the operation to be performed as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
-*
-*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
-*
-*              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
-*
-*           Unchanged on exit.
-*
-*  M      - INTEGER.
-*           On entry, M specifies the number of rows of the matrix A.
-*           M must be at least zero.
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the number of columns of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  KL     - INTEGER.
-*           On entry, KL specifies the number of sub-diagonals of the
-*           matrix A. KL must satisfy  0 .le. KL.
-*           Unchanged on exit.
-*
-*  KU     - INTEGER.
-*           On entry, KU specifies the number of super-diagonals of the
-*           matrix A. KU must satisfy  0 .le. KU.
-*           Unchanged on exit.
-*
-*  ALPHA  - REAL            .
-*           On entry, ALPHA specifies the scalar alpha.
-*           Unchanged on exit.
-*
-*  A      - REAL             array of DIMENSION ( LDA, n ).
-*           Before entry, the leading ( kl + ku + 1 ) by n part of the
-*           array A must contain the matrix of coefficients, supplied
-*           column by column, with the leading diagonal of the matrix in
-*           row ( ku + 1 ) of the array, the first super-diagonal
-*           starting at position 2 in row ku, the first sub-diagonal
-*           starting at position 1 in row ( ku + 2 ), and so on.
-*           Elements in the array A that do not correspond to elements
-*           in the band matrix (such as the top left ku by ku triangle)
-*           are not referenced.
-*           The following program segment will transfer a band matrix
-*           from conventional full matrix storage to band storage:
-*
-*                 DO 20, J = 1, N
-*                    K = KU + 1 - J
-*                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
-*                       A( K + I, J ) = matrix( I, J )
-*              10    CONTINUE
-*              20 CONTINUE
-*
-*           Unchanged on exit.
-*
-*  LDA    - INTEGER.
-*           On entry, LDA specifies the first dimension of A as declared
-*           in the calling (sub) program. LDA must be at least
-*           ( kl + ku + 1 ).
-*           Unchanged on exit.
-*
-*  X      - REAL             array of DIMENSION at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
-*           Before entry, the incremented array X must contain the
-*           vector x.
-*           Unchanged on exit.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  BETA   - REAL            .
-*           On entry, BETA specifies the scalar beta. When BETA is
-*           supplied as zero then Y need not be set on input.
-*           Unchanged on exit.
-*
-*  Y      - REAL             array of DIMENSION at least
-*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
-*           Before entry, the incremented array Y must contain the
-*           vector y. On exit, Y is overwritten by the updated vector y.
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      REAL ONE,ZERO
-      PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
-*     ..
-*     .. Local Scalars ..
-      REAL TEMP
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC MAX,MIN
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +    .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 1
-      ELSE IF (M.LT.0) THEN
-          INFO = 2
-      ELSE IF (N.LT.0) THEN
-          INFO = 3
-      ELSE IF (KL.LT.0) THEN
-          INFO = 4
-      ELSE IF (KU.LT.0) THEN
-          INFO = 5
-      ELSE IF (LDA.LT. (KL+KU+1)) THEN
-          INFO = 8
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 10
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 13
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('SGBMV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
-     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
-*
-*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
-*     up the start points in  X  and  Y.
-*
-      IF (LSAME(TRANS,'N')) THEN
-          LENX = N
-          LENY = M
-      ELSE
-          LENX = M
-          LENY = N
-      END IF
-      IF (INCX.GT.0) THEN
-          KX = 1
-      ELSE
-          KX = 1 - (LENX-1)*INCX
-      END IF
-      IF (INCY.GT.0) THEN
-          KY = 1
-      ELSE
-          KY = 1 - (LENY-1)*INCY
-      END IF
-*
-*     Start the operations. In this version the elements of A are
-*     accessed sequentially with one pass through the band part of A.
-*
-*     First form  y := beta*y.
-*
-      IF (BETA.NE.ONE) THEN
-          IF (INCY.EQ.1) THEN
-              IF (BETA.EQ.ZERO) THEN
-                  DO 10 I = 1,LENY
-                      Y(I) = ZERO
-   10             CONTINUE
-              ELSE
-                  DO 20 I = 1,LENY
-                      Y(I) = BETA*Y(I)
-   20             CONTINUE
-              END IF
-          ELSE
-              IY = KY
-              IF (BETA.EQ.ZERO) THEN
-                  DO 30 I = 1,LENY
-                      Y(IY) = ZERO
-                      IY = IY + INCY
-   30             CONTINUE
-              ELSE
-                  DO 40 I = 1,LENY
-                      Y(IY) = BETA*Y(IY)
-                      IY = IY + INCY
-   40             CONTINUE
-              END IF
-          END IF
-      END IF
-      IF (ALPHA.EQ.ZERO) RETURN
-      KUP1 = KU + 1
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  y := alpha*A*x + y.
-*
-          JX = KX
-          IF (INCY.EQ.1) THEN
-              DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  IF (J.GT.KU) KY = KY + INCY
-   80         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  y := alpha*A'*x + y.
-*
-          JY = KY
-          IF (INCX.EQ.1) THEN
-              DO 100 J = 1,N
-                  TEMP = ZERO
-                  K = KUP1 - J
-                  DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
-                      TEMP = TEMP + A(K+I,J)*X(I)
-   90             CONTINUE
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-  100         CONTINUE
-          ELSE
-              DO 120 J = 1,N
-                  TEMP = ZERO
-                  IX = KX
-                  K = KUP1 - J
-                  DO 110 I = MAX(1,J-KU),MIN(M,J+KL)
-                      TEMP = TEMP + A(K+I,J)*X(IX)
-                      IX = IX + INCX
-  110             CONTINUE
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-                  IF (J.GT.KU) KX = KX + INCX
-  120         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of SGBMV .
-*
-      END
diff --git a/blas/zgbmv.f b/blas/zgbmv.f
deleted file mode 100644
index 5a2228d04..000000000
--- a/blas/zgbmv.f
+++ /dev/null
@@ -1,322 +0,0 @@
-      SUBROUTINE ZGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-*     .. Scalar Arguments ..
-      DOUBLE COMPLEX ALPHA,BETA
-      INTEGER INCX,INCY,KL,KU,LDA,M,N
-      CHARACTER TRANS
-*     ..
-*     .. Array Arguments ..
-      DOUBLE COMPLEX A(LDA,*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  ZGBMV  performs one of the matrix-vector operations
-*
-*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or
-*
-*     y := alpha*conjg( A' )*x + beta*y,
-*
-*  where alpha and beta are scalars, x and y are vectors and A is an
-*  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
-*
-*  Arguments
-*  ==========
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the operation to be performed as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
-*
-*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
-*
-*              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y.
-*
-*           Unchanged on exit.
-*
-*  M      - INTEGER.
-*           On entry, M specifies the number of rows of the matrix A.
-*           M must be at least zero.
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the number of columns of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  KL     - INTEGER.
-*           On entry, KL specifies the number of sub-diagonals of the
-*           matrix A. KL must satisfy  0 .le. KL.
-*           Unchanged on exit.
-*
-*  KU     - INTEGER.
-*           On entry, KU specifies the number of super-diagonals of the
-*           matrix A. KU must satisfy  0 .le. KU.
-*           Unchanged on exit.
-*
-*  ALPHA  - COMPLEX*16      .
-*           On entry, ALPHA specifies the scalar alpha.
-*           Unchanged on exit.
-*
-*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
-*           Before entry, the leading ( kl + ku + 1 ) by n part of the
-*           array A must contain the matrix of coefficients, supplied
-*           column by column, with the leading diagonal of the matrix in
-*           row ( ku + 1 ) of the array, the first super-diagonal
-*           starting at position 2 in row ku, the first sub-diagonal
-*           starting at position 1 in row ( ku + 2 ), and so on.
-*           Elements in the array A that do not correspond to elements
-*           in the band matrix (such as the top left ku by ku triangle)
-*           are not referenced.
-*           The following program segment will transfer a band matrix
-*           from conventional full matrix storage to band storage:
-*
-*                 DO 20, J = 1, N
-*                    K = KU + 1 - J
-*                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
-*                       A( K + I, J ) = matrix( I, J )
-*              10    CONTINUE
-*              20 CONTINUE
-*
-*           Unchanged on exit.
-*
-*  LDA    - INTEGER.
-*           On entry, LDA specifies the first dimension of A as declared
-*           in the calling (sub) program. LDA must be at least
-*           ( kl + ku + 1 ).
-*           Unchanged on exit.
-*
-*  X      - COMPLEX*16       array of DIMENSION at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
-*           Before entry, the incremented array X must contain the
-*           vector x.
-*           Unchanged on exit.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  BETA   - COMPLEX*16      .
-*           On entry, BETA specifies the scalar beta. When BETA is
-*           supplied as zero then Y need not be set on input.
-*           Unchanged on exit.
-*
-*  Y      - COMPLEX*16       array of DIMENSION at least
-*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
-*           and at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
-*           Before entry, the incremented array Y must contain the
-*           vector y. On exit, Y is overwritten by the updated vector y.
-*
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE COMPLEX ONE
-      PARAMETER (ONE= (1.0D+0,0.0D+0))
-      DOUBLE COMPLEX ZERO
-      PARAMETER (ZERO= (0.0D+0,0.0D+0))
-*     ..
-*     .. Local Scalars ..
-      DOUBLE COMPLEX TEMP
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
-      LOGICAL NOCONJ
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC DCONJG,MAX,MIN
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +    .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 1
-      ELSE IF (M.LT.0) THEN
-          INFO = 2
-      ELSE IF (N.LT.0) THEN
-          INFO = 3
-      ELSE IF (KL.LT.0) THEN
-          INFO = 4
-      ELSE IF (KU.LT.0) THEN
-          INFO = 5
-      ELSE IF (LDA.LT. (KL+KU+1)) THEN
-          INFO = 8
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 10
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 13
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('ZGBMV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
-     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
-*
-      NOCONJ = LSAME(TRANS,'T')
-*
-*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
-*     up the start points in  X  and  Y.
-*
-      IF (LSAME(TRANS,'N')) THEN
-          LENX = N
-          LENY = M
-      ELSE
-          LENX = M
-          LENY = N
-      END IF
-      IF (INCX.GT.0) THEN
-          KX = 1
-      ELSE
-          KX = 1 - (LENX-1)*INCX
-      END IF
-      IF (INCY.GT.0) THEN
-          KY = 1
-      ELSE
-          KY = 1 - (LENY-1)*INCY
-      END IF
-*
-*     Start the operations. In this version the elements of A are
-*     accessed sequentially with one pass through the band part of A.
-*
-*     First form  y := beta*y.
-*
-      IF (BETA.NE.ONE) THEN
-          IF (INCY.EQ.1) THEN
-              IF (BETA.EQ.ZERO) THEN
-                  DO 10 I = 1,LENY
-                      Y(I) = ZERO
-   10             CONTINUE
-              ELSE
-                  DO 20 I = 1,LENY
-                      Y(I) = BETA*Y(I)
-   20             CONTINUE
-              END IF
-          ELSE
-              IY = KY
-              IF (BETA.EQ.ZERO) THEN
-                  DO 30 I = 1,LENY
-                      Y(IY) = ZERO
-                      IY = IY + INCY
-   30             CONTINUE
-              ELSE
-                  DO 40 I = 1,LENY
-                      Y(IY) = BETA*Y(IY)
-                      IY = IY + INCY
-   40             CONTINUE
-              END IF
-          END IF
-      END IF
-      IF (ALPHA.EQ.ZERO) RETURN
-      KUP1 = KU + 1
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  y := alpha*A*x + y.
-*
-          JX = KX
-          IF (INCY.EQ.1) THEN
-              DO 60 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      K = KUP1 - J
-                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(I) = Y(I) + TEMP*A(K+I,J)
-   50                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF (X(JX).NE.ZERO) THEN
-                      TEMP = ALPHA*X(JX)
-                      IY = KY
-                      K = KUP1 - J
-                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
-                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  IF (J.GT.KU) KY = KY + INCY
-   80         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y.
-*
-          JY = KY
-          IF (INCX.EQ.1) THEN
-              DO 110 J = 1,N
-                  TEMP = ZERO
-                  K = KUP1 - J
-                  IF (NOCONJ) THEN
-                      DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + A(K+I,J)*X(I)
-   90                 CONTINUE
-                  ELSE
-                      DO 100 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + DCONJG(A(K+I,J))*X(I)
-  100                 CONTINUE
-                  END IF
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-  110         CONTINUE
-          ELSE
-              DO 140 J = 1,N
-                  TEMP = ZERO
-                  IX = KX
-                  K = KUP1 - J
-                  IF (NOCONJ) THEN
-                      DO 120 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + A(K+I,J)*X(IX)
-                          IX = IX + INCX
-  120                 CONTINUE
-                  ELSE
-                      DO 130 I = MAX(1,J-KU),MIN(M,J+KL)
-                          TEMP = TEMP + DCONJG(A(K+I,J))*X(IX)
-                          IX = IX + INCX
-  130                 CONTINUE
-                  END IF
-                  Y(JY) = Y(JY) + ALPHA*TEMP
-                  JY = JY + INCY
-                  IF (J.GT.KU) KX = KX + INCX
-  140         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of ZGBMV .
-*
-      END