Adds a modified f2c-generated C implmentation for BLAS.

This adds an optional implementation for the BLAS library that does
not require the use of a FORTRAN compiler. It can be enabled with
EIGEN_USE_F2C_BLAS.

The C implementation uses the standard gfortran calling convention
and does not require the use of -ff2c when compiled with gfortran.

40 files changed, 6647 insertions, 14 deletions

diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt
index a9bc05137..2bc956a64 100644
--- a/blas/CMakeLists.txt
+++ b/blas/CMakeLists.txt
@@ -16,21 +16,31 @@ add_custom_target(blas)
 
 set(EigenBlas_SRCS single.cpp double.cpp complex_single.cpp complex_double.cpp xerbla.cpp)
 
-if(EIGEN_Fortran_COMPILER_WORKS)
-
-set(EigenBlas_SRCS ${EigenBlas_SRCS}
-    complexdots.f
-    srotm.f srotmg.f drotm.f drotmg.f
-    lsame.f  dspmv.f ssbmv.f
-    chbmv.f  sspmv.f
-    zhbmv.f  chpmv.f dsbmv.f
-    zhpmv.f
-    dtbmv.f stbmv.f ctbmv.f ztbmv.f
-)
+if(EIGEN_USE_F2C_BLAS)
+  set(EigenBlas_SRCS ${EigenBlas_SRCS}
+    f2c/complexdots.c
+    f2c/srotm.c f2c/srotmg.c f2c/drotm.c f2c/drotmg.c
+    f2c/lsame.c  f2c/dspmv.c f2c/ssbmv.c
+    f2c/chbmv.c  f2c/sspmv.c
+    f2c/zhbmv.c  f2c/chpmv.c f2c/dsbmv.c
+    f2c/zhpmv.c
+    f2c/dtbmv.c f2c/stbmv.c f2c/ctbmv.c f2c/ztbmv.c
+    f2c/d_cnjg.c f2c/r_cnjg.c
+    )
 else()
-
-message(WARNING " No fortran compiler has been detected, the blas build will be incomplete.")
-
+  if (EIGEN_Fortran_COMPILER_WORKS)
+    set(EigenBlas_SRCS ${EigenBlas_SRCS}
+      fortran/complexdots.f
+      fortran/srotm.f fortran/srotmg.f fortran/drotm.f fortran/drotmg.f
+      fortran/lsame.f  fortran/dspmv.f fortran/ssbmv.f
+      fortran/chbmv.f  fortran/sspmv.f
+      fortran/zhbmv.f  fortran/chpmv.f fortran/dsbmv.f
+      fortran/zhpmv.f
+      fortran/dtbmv.f fortran/stbmv.f fortran/ctbmv.f fortran/ztbmv.f
+      )
+  else()
+    message(WARNING " No Fortran compiler has been detected, the blas build will be incomplete. Define EIGEN_USE_F2C_BLAS to build BLAS without Fortran")
+  endif()
 endif()
 
 add_library(eigen_blas_static ${EigenBlas_SRCS})
diff --git a/blas/f2c/chbmv.c b/blas/f2c/chbmv.c
new file mode 100644
index 000000000..f218fe3f5
--- /dev/null
+++ b/blas/f2c/chbmv.c
@@ -0,0 +1,487 @@
+/* chbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int chbmv_(char *uplo, integer *n, integer *k, complex *
+	alpha, complex *a, integer *lda, complex *x, integer *incx, complex *
+	beta, complex *y, integer *incy, ftnlen uplo_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
+    real r__1;
+    complex q__1, q__2, q__3, q__4;
+
+    /* Builtin functions */
+    void r_cnjg(complex *, complex *);
+
+    /* Local variables */
+    integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
+    complex temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  CHBMV  performs the matrix-vector  operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n hermitian band matrix, with k super-diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the band matrix A is being supplied as */
+/*           follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  being supplied. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  being supplied. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry, K specifies the number of super-diagonals of the */
+/*           matrix A. K must satisfy  0 .le. K. */
+/*           Unchanged on exit. */
+
+/*  ALPHA  - COMPLEX         . */
+/*           On entry, ALPHA specifies the scalar alpha. */
+/*           Unchanged on exit. */
+
+/*  A      - COMPLEX          array of DIMENSION ( LDA, n ). */
+/*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the hermitian matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer the upper */
+/*           triangular part of a hermitian band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the hermitian matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer the lower */
+/*           triangular part of a hermitian band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set and are assumed to be zero. */
+/*           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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - COMPLEX          array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           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. */
+/*           Unchanged on exit. */
+
+/*  Y      - COMPLEX          array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCY ) ). */
+/*           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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+    --y;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*k < 0) {
+	info = 3;
+    } else if (*lda < *k + 1) {
+	info = 6;
+    } else if (*incx == 0) {
+	info = 8;
+    } else if (*incy == 0) {
+	info = 11;
+    }
+    if (info != 0) {
+	xerbla_("CHBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f && (beta->r == 1.f && 
+                                                           beta->i == 0.f))) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array A */
+/*     are accessed sequentially with one pass through A. */
+
+/*     First form  y := beta*y. */
+
+    if (beta->r != 1.f || beta->i != 0.f) {
+	if (*incy == 1) {
+	    if (beta->r == 0.f && beta->i == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    y[i__2].r = 0.f, y[i__2].i = 0.f;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    i__3 = i__;
+		    q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (beta->r == 0.f && beta->i == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    y[i__2].r = 0.f, y[i__2].i = 0.f;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    i__3 = iy;
+		    q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (alpha->r == 0.f && alpha->i == 0.f) {
+	return 0;
+    }
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when upper triangle of A is stored. */
+
+	kplus1 = *k + 1;
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__2 = 1, i__3 = j - *k;
+		i__4 = j - 1;
+		for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+		    i__2 = i__;
+		    i__3 = i__;
+		    i__5 = l + i__ + j * a_dim1;
+		    q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
+		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+		    i__2 = i__;
+		    q__2.r = q__3.r * x[i__2].r - q__3.i * x[i__2].i, q__2.i =
+			     q__3.r * x[i__2].i + q__3.i * x[i__2].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+/* L50: */
+		}
+		i__4 = j;
+		i__2 = j;
+		i__3 = kplus1 + j * a_dim1;
+		r__1 = a[i__3].r;
+		q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
+		q__2.r = y[i__2].r + q__3.r, q__2.i = y[i__2].i + q__3.i;
+		q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+		y[i__4].r = q__1.r, y[i__4].i = q__1.i;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__4 = jx;
+		q__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, q__1.i =
+			 alpha->r * x[i__4].i + alpha->i * x[i__4].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		ix = kx;
+		iy = ky;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__4 = 1, i__2 = j - *k;
+		i__3 = j - 1;
+		for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+		    i__4 = iy;
+		    i__2 = iy;
+		    i__5 = l + i__ + j * a_dim1;
+		    q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i;
+		    y[i__4].r = q__1.r, y[i__4].i = q__1.i;
+		    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+		    i__4 = ix;
+		    q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i =
+			     q__3.r * x[i__4].i + q__3.i * x[i__4].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		i__3 = jy;
+		i__4 = jy;
+		i__2 = kplus1 + j * a_dim1;
+		r__1 = a[i__2].r;
+		q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
+		q__2.r = y[i__4].r + q__3.r, q__2.i = y[i__4].i + q__3.i;
+		q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+		y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		jx += *incx;
+		jy += *incy;
+		if (j > *k) {
+		    kx += *incx;
+		    ky += *incy;
+		}
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when lower triangle of A is stored. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__3 = j;
+		q__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, q__1.i =
+			 alpha->r * x[i__3].i + alpha->i * x[i__3].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		i__3 = j;
+		i__4 = j;
+		i__2 = j * a_dim1 + 1;
+		r__1 = a[i__2].r;
+		q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
+		q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		l = 1 - j;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    i__4 = i__;
+		    i__2 = i__;
+		    i__5 = l + i__ + j * a_dim1;
+		    q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i;
+		    y[i__4].r = q__1.r, y[i__4].i = q__1.i;
+		    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+		    i__4 = i__;
+		    q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i =
+			     q__3.r * x[i__4].i + q__3.i * x[i__4].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+/* L90: */
+		}
+		i__3 = j;
+		i__4 = j;
+		q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__3 = jx;
+		q__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, q__1.i =
+			 alpha->r * x[i__3].i + alpha->i * x[i__3].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		i__3 = jy;
+		i__4 = jy;
+		i__2 = j * a_dim1 + 1;
+		r__1 = a[i__2].r;
+		q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
+		q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		l = 1 - j;
+		ix = jx;
+		iy = jy;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    ix += *incx;
+		    iy += *incy;
+		    i__4 = iy;
+		    i__2 = iy;
+		    i__5 = l + i__ + j * a_dim1;
+		    q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i;
+		    y[i__4].r = q__1.r, y[i__4].i = q__1.i;
+		    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+		    i__4 = ix;
+		    q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i =
+			     q__3.r * x[i__4].i + q__3.i * x[i__4].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+/* L110: */
+		}
+		i__3 = jy;
+		i__4 = jy;
+		q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		jx += *incx;
+		jy += *incy;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of CHBMV . */
+
+} /* chbmv_ */
+
diff --git a/blas/f2c/chpmv.c b/blas/f2c/chpmv.c
new file mode 100644
index 000000000..65bab1c7f
--- /dev/null
+++ b/blas/f2c/chpmv.c
@@ -0,0 +1,438 @@
+/* chpmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int chpmv_(char *uplo, integer *n, complex *alpha, complex *
+	ap, complex *x, integer *incx, complex *beta, complex *y, integer *
+	incy, ftnlen uplo_len)
+{
+    /* System generated locals */
+    integer i__1, i__2, i__3, i__4, i__5;
+    real r__1;
+    complex q__1, q__2, q__3, q__4;
+
+    /* Builtin functions */
+    void r_cnjg(complex *, complex *);
+
+    /* Local variables */
+    integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
+    complex temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  CHPMV  performs the matrix-vector operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n hermitian matrix, supplied in packed form. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the matrix A is supplied in the packed */
+/*           array AP as follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  supplied in AP. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  supplied in AP. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  ALPHA  - COMPLEX         . */
+/*           On entry, ALPHA specifies the scalar alpha. */
+/*           Unchanged on exit. */
+
+/*  AP     - COMPLEX          array of DIMENSION at least */
+/*           ( ( n*( n + 1 ) )/2 ). */
+/*           Before entry with UPLO = 'U' or 'u', the array AP must */
+/*           contain the upper triangular part of the hermitian matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
+/*           and a( 2, 2 ) respectively, and so on. */
+/*           Before entry with UPLO = 'L' or 'l', the array AP must */
+/*           contain the lower triangular part of the hermitian matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
+/*           and a( 3, 1 ) respectively, and so on. */
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set and are assumed to be zero. */
+/*           Unchanged on exit. */
+
+/*  X      - COMPLEX          array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element 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 + ( n - 1 )*abs( INCY ) ). */
+/*           Before entry, the incremented array Y must contain the n */
+/*           element 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --y;
+    --x;
+    --ap;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*incx == 0) {
+	info = 6;
+    } else if (*incy == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("CHPMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f && (beta->r == 1.f && 
+                                                           beta->i == 0.f))) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+/*     First form  y := beta*y. */
+
+    if (beta->r != 1.f || beta->i != 0.f) {
+	if (*incy == 1) {
+	    if (beta->r == 0.f && beta->i == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    y[i__2].r = 0.f, y[i__2].i = 0.f;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    i__3 = i__;
+		    q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (beta->r == 0.f && beta->i == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    y[i__2].r = 0.f, y[i__2].i = 0.f;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    i__3 = iy;
+		    q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (alpha->r == 0.f && alpha->i == 0.f) {
+	return 0;
+    }
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when AP contains the upper triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		k = kk;
+		i__2 = j - 1;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    i__3 = i__;
+		    i__4 = i__;
+		    i__5 = k;
+		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		    r_cnjg(&q__3, &ap[k]);
+		    i__3 = i__;
+		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
+			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    ++k;
+/* L50: */
+		}
+		i__2 = j;
+		i__3 = j;
+		i__4 = kk + j - 1;
+		r__1 = ap[i__4].r;
+		q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
+		q__2.r = y[i__3].r + q__3.r, q__2.i = y[i__3].i + q__3.i;
+		q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		kk += j;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		ix = kx;
+		iy = ky;
+		i__2 = kk + j - 2;
+		for (k = kk; k <= i__2; ++k) {
+		    i__3 = iy;
+		    i__4 = iy;
+		    i__5 = k;
+		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		    r_cnjg(&q__3, &ap[k]);
+		    i__3 = ix;
+		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
+			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		i__2 = jy;
+		i__3 = jy;
+		i__4 = kk + j - 1;
+		r__1 = ap[i__4].r;
+		q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
+		q__2.r = y[i__3].r + q__3.r, q__2.i = y[i__3].i + q__3.i;
+		q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when AP contains the lower triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		i__2 = j;
+		i__3 = j;
+		i__4 = kk;
+		r__1 = ap[i__4].r;
+		q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
+		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
+		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		k = kk + 1;
+		i__2 = *n;
+		for (i__ = j + 1; i__ <= i__2; ++i__) {
+		    i__3 = i__;
+		    i__4 = i__;
+		    i__5 = k;
+		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		    r_cnjg(&q__3, &ap[k]);
+		    i__3 = i__;
+		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
+			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    ++k;
+/* L90: */
+		}
+		i__2 = j;
+		i__3 = j;
+		q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
+		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		kk += *n - j + 1;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = q__1.r, temp1.i = q__1.i;
+		temp2.r = 0.f, temp2.i = 0.f;
+		i__2 = jy;
+		i__3 = jy;
+		i__4 = kk;
+		r__1 = ap[i__4].r;
+		q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
+		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
+		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		ix = jx;
+		iy = jy;
+		i__2 = kk + *n - j;
+		for (k = kk + 1; k <= i__2; ++k) {
+		    ix += *incx;
+		    iy += *incy;
+		    i__3 = iy;
+		    i__4 = iy;
+		    i__5 = k;
+		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
+		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
+		    r_cnjg(&q__3, &ap[k]);
+		    i__3 = ix;
+		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
+			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
+		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+/* L110: */
+		}
+		i__2 = jy;
+		i__3 = jy;
+		q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
+		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
+		jx += *incx;
+		jy += *incy;
+		kk += *n - j + 1;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of CHPMV . */
+
+} /* chpmv_ */
+
diff --git a/blas/f2c/complexdots.c b/blas/f2c/complexdots.c
new file mode 100644
index 000000000..a856a231c
--- /dev/null
+++ b/blas/f2c/complexdots.c
@@ -0,0 +1,84 @@
+/* This file has been modified to use the standard gfortran calling
+   convention, rather than the f2c calling convention.
+
+   It does not require -ff2c when compiled with gfortran.
+*/
+
+/* complexdots.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+complex cdotc_(integer *n, complex *cx, integer 
+	*incx, complex *cy, integer *incy)
+{
+    complex res;
+    extern /* Subroutine */ int cdotcw_(integer *, complex *, integer *, 
+	    complex *, integer *, complex *);
+
+    /* Parameter adjustments */
+    --cy;
+    --cx;
+
+    /* Function Body */
+    cdotcw_(n, &cx[1], incx, &cy[1], incy, &res);
+    return res;
+} /* cdotc_ */
+
+complex cdotu_(integer *n, complex *cx, integer 
+	*incx, complex *cy, integer *incy)
+{
+    complex res;
+    extern /* Subroutine */ int cdotuw_(integer *, complex *, integer *, 
+	    complex *, integer *, complex *);
+
+    /* Parameter adjustments */
+    --cy;
+    --cx;
+
+    /* Function Body */
+    cdotuw_(n, &cx[1], incx, &cy[1], incy, &res);
+    return res;
+} /* cdotu_ */
+
+doublecomplex zdotc_(integer *n, doublecomplex *cx, integer *incx, 
+                     doublecomplex *cy, integer *incy)
+{
+    doublecomplex res;
+    extern /* Subroutine */ int zdotcw_(integer *, doublecomplex *, integer *,
+	     doublecomplex *, integer *, doublecomplex *);
+
+    /* Parameter adjustments */
+    --cy;
+    --cx;
+
+    /* Function Body */
+    zdotcw_(n, &cx[1], incx, &cy[1], incy, &res);
+    return res;
+} /* zdotc_ */
+
+doublecomplex zdotu_(integer *n, doublecomplex *cx, integer *incx, 
+                     doublecomplex *cy, integer *incy)
+{
+    doublecomplex res;
+    extern /* Subroutine */ int zdotuw_(integer *, doublecomplex *, integer *,
+	     doublecomplex *, integer *, doublecomplex *);
+
+    /* Parameter adjustments */
+    --cy;
+    --cx;
+
+    /* Function Body */
+    zdotuw_(n, &cx[1], incx, &cy[1], incy, &res);
+    return res;
+} /* zdotu_ */
+
diff --git a/blas/f2c/ctbmv.c b/blas/f2c/ctbmv.c
new file mode 100644
index 000000000..790fd581f
--- /dev/null
+++ b/blas/f2c/ctbmv.c
@@ -0,0 +1,647 @@
+/* ctbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int ctbmv_(char *uplo, char *trans, char *diag, integer *n, 
+	integer *k, complex *a, integer *lda, complex *x, integer *incx, 
+	ftnlen uplo_len, ftnlen trans_len, ftnlen diag_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
+    complex q__1, q__2, q__3;
+
+    /* Builtin functions */
+    void r_cnjg(complex *, complex *);
+
+    /* Local variables */
+    integer i__, j, l, ix, jx, kx, info;
+    complex temp;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+    logical noconj, nounit;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  CTBMV  performs one of the matrix-vector operations */
+
+/*     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*x, */
+
+/*  where x is an n element vector and  A is an n by n unit, or non-unit, */
+/*  upper or lower triangular band matrix, with ( k + 1 ) diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the matrix is an upper or */
+/*           lower triangular matrix as follows: */
+
+/*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
+
+/*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
+
+/*           Unchanged on exit. */
+
+/*  TRANS  - CHARACTER*1. */
+/*           On entry, TRANS specifies the operation to be performed as */
+/*           follows: */
+
+/*              TRANS = 'N' or 'n'   x := A*x. */
+
+/*              TRANS = 'T' or 't'   x := A'*x. */
+
+/*              TRANS = 'C' or 'c'   x := conjg( A' )*x. */
+
+/*           Unchanged on exit. */
+
+/*  DIAG   - CHARACTER*1. */
+/*           On entry, DIAG specifies whether or not A is unit */
+/*           triangular as follows: */
+
+/*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
+
+/*              DIAG = 'N' or 'n'   A is not assumed to be unit */
+/*                                  triangular. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry with UPLO = 'U' or 'u', K specifies the number of */
+/*           super-diagonals of the matrix A. */
+/*           On entry with UPLO = 'L' or 'l', K specifies the number of */
+/*           sub-diagonals of the matrix A. */
+/*           K must satisfy  0 .le. K. */
+/*           Unchanged on exit. */
+
+/*  A      - COMPLEX          array of DIMENSION ( LDA, n ). */
+/*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer an upper */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer a lower */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Note that when DIAG = 'U' or 'u' the elements of the array A */
+/*           corresponding to the diagonal elements of the matrix are not */
+/*           referenced, but are assumed to be unity. */
+/*           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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - COMPLEX          array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element vector x. On exit, X is overwritten with the */
+/*           tranformed vector x. */
+
+/*  INCX   - INTEGER. */
+/*           On entry, INCX specifies the increment for the elements of */
+/*           X. INCX 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, 
+	    "T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
+	    ftnlen)1)) {
+	info = 2;
+    } else if (! lsame_(diag, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(diag, 
+	    "N", (ftnlen)1, (ftnlen)1)) {
+	info = 3;
+    } else if (*n < 0) {
+	info = 4;
+    } else if (*k < 0) {
+	info = 5;
+    } else if (*lda < *k + 1) {
+	info = 7;
+    } else if (*incx == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("CTBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0) {
+	return 0;
+    }
+
+    noconj = lsame_(trans, "T", (ftnlen)1, (ftnlen)1);
+    nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
+
+/*     Set up the start point in X if the increment is not unity. This */
+/*     will be  ( N - 1 )*INCX   too small for descending loops. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of A are */
+/*     accessed sequentially with one pass through A. */
+
+    if (lsame_(trans, "N", (ftnlen)1, (ftnlen)1)) {
+
+/*         Form  x := A*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = j;
+		    if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
+			i__2 = j;
+			temp.r = x[i__2].r, temp.i = x[i__2].i;
+			l = kplus1 - j;
+/* Computing MAX */
+			i__2 = 1, i__3 = j - *k;
+			i__4 = j - 1;
+			for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+			    i__2 = i__;
+			    i__3 = i__;
+			    i__5 = l + i__ + j * a_dim1;
+			    q__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i, 
+				    q__2.i = temp.r * a[i__5].i + temp.i * a[
+				    i__5].r;
+			    q__1.r = x[i__3].r + q__2.r, q__1.i = x[i__3].i + 
+				    q__2.i;
+			    x[i__2].r = q__1.r, x[i__2].i = q__1.i;
+/* L10: */
+			}
+			if (nounit) {
+			    i__4 = j;
+			    i__2 = j;
+			    i__3 = kplus1 + j * a_dim1;
+			    q__1.r = x[i__2].r * a[i__3].r - x[i__2].i * a[
+				    i__3].i, q__1.i = x[i__2].r * a[i__3].i + 
+				    x[i__2].i * a[i__3].r;
+			    x[i__4].r = q__1.r, x[i__4].i = q__1.i;
+			}
+		    }
+/* L20: */
+		}
+	    } else {
+		jx = kx;
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    i__4 = jx;
+		    if (x[i__4].r != 0.f || x[i__4].i != 0.f) {
+			i__4 = jx;
+			temp.r = x[i__4].r, temp.i = x[i__4].i;
+			ix = kx;
+			l = kplus1 - j;
+/* Computing MAX */
+			i__4 = 1, i__2 = j - *k;
+			i__3 = j - 1;
+			for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+			    i__4 = ix;
+			    i__2 = ix;
+			    i__5 = l + i__ + j * a_dim1;
+			    q__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i, 
+				    q__2.i = temp.r * a[i__5].i + temp.i * a[
+				    i__5].r;
+			    q__1.r = x[i__2].r + q__2.r, q__1.i = x[i__2].i + 
+				    q__2.i;
+			    x[i__4].r = q__1.r, x[i__4].i = q__1.i;
+			    ix += *incx;
+/* L30: */
+			}
+			if (nounit) {
+			    i__3 = jx;
+			    i__4 = jx;
+			    i__2 = kplus1 + j * a_dim1;
+			    q__1.r = x[i__4].r * a[i__2].r - x[i__4].i * a[
+				    i__2].i, q__1.i = x[i__4].r * a[i__2].i + 
+				    x[i__4].i * a[i__2].r;
+			    x[i__3].r = q__1.r, x[i__3].i = q__1.i;
+			}
+		    }
+		    jx += *incx;
+		    if (j > *k) {
+			kx += *incx;
+		    }
+/* L40: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    i__1 = j;
+		    if (x[i__1].r != 0.f || x[i__1].i != 0.f) {
+			i__1 = j;
+			temp.r = x[i__1].r, temp.i = x[i__1].i;
+			l = 1 - j;
+/* Computing MIN */
+			i__1 = *n, i__3 = j + *k;
+			i__4 = j + 1;
+			for (i__ = min(i__1,i__3); i__ >= i__4; --i__) {
+			    i__1 = i__;
+			    i__3 = i__;
+			    i__2 = l + i__ + j * a_dim1;
+			    q__2.r = temp.r * a[i__2].r - temp.i * a[i__2].i, 
+				    q__2.i = temp.r * a[i__2].i + temp.i * a[
+				    i__2].r;
+			    q__1.r = x[i__3].r + q__2.r, q__1.i = x[i__3].i + 
+				    q__2.i;
+			    x[i__1].r = q__1.r, x[i__1].i = q__1.i;
+/* L50: */
+			}
+			if (nounit) {
+			    i__4 = j;
+			    i__1 = j;
+			    i__3 = j * a_dim1 + 1;
+			    q__1.r = x[i__1].r * a[i__3].r - x[i__1].i * a[
+				    i__3].i, q__1.i = x[i__1].r * a[i__3].i + 
+				    x[i__1].i * a[i__3].r;
+			    x[i__4].r = q__1.r, x[i__4].i = q__1.i;
+			}
+		    }
+/* L60: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    i__4 = jx;
+		    if (x[i__4].r != 0.f || x[i__4].i != 0.f) {
+			i__4 = jx;
+			temp.r = x[i__4].r, temp.i = x[i__4].i;
+			ix = kx;
+			l = 1 - j;
+/* Computing MIN */
+			i__4 = *n, i__1 = j + *k;
+			i__3 = j + 1;
+			for (i__ = min(i__4,i__1); i__ >= i__3; --i__) {
+			    i__4 = ix;
+			    i__1 = ix;
+			    i__2 = l + i__ + j * a_dim1;
+			    q__2.r = temp.r * a[i__2].r - temp.i * a[i__2].i, 
+				    q__2.i = temp.r * a[i__2].i + temp.i * a[
+				    i__2].r;
+			    q__1.r = x[i__1].r + q__2.r, q__1.i = x[i__1].i + 
+				    q__2.i;
+			    x[i__4].r = q__1.r, x[i__4].i = q__1.i;
+			    ix -= *incx;
+/* L70: */
+			}
+			if (nounit) {
+			    i__3 = jx;
+			    i__4 = jx;
+			    i__1 = j * a_dim1 + 1;
+			    q__1.r = x[i__4].r * a[i__1].r - x[i__4].i * a[
+				    i__1].i, q__1.i = x[i__4].r * a[i__1].i + 
+				    x[i__4].i * a[i__1].r;
+			    x[i__3].r = q__1.r, x[i__3].i = q__1.i;
+			}
+		    }
+		    jx -= *incx;
+		    if (*n - j >= *k) {
+			kx -= *incx;
+		    }
+/* L80: */
+		}
+	    }
+	}
+    } else {
+
+/*        Form  x := A'*x  or  x := conjg( A' )*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    i__3 = j;
+		    temp.r = x[i__3].r, temp.i = x[i__3].i;
+		    l = kplus1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__3 = kplus1 + j * a_dim1;
+			    q__1.r = temp.r * a[i__3].r - temp.i * a[i__3].i, 
+				    q__1.i = temp.r * a[i__3].i + temp.i * a[
+				    i__3].r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    i__4 = l + i__ + j * a_dim1;
+			    i__1 = i__;
+			    q__2.r = a[i__4].r * x[i__1].r - a[i__4].i * x[
+				    i__1].i, q__2.i = a[i__4].r * x[i__1].i + 
+				    a[i__4].i * x[i__1].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+/* L90: */
+			}
+		    } else {
+			if (nounit) {
+			    r_cnjg(&q__2, &a[kplus1 + j * a_dim1]);
+			    q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
+				    q__1.i = temp.r * q__2.i + temp.i * 
+				    q__2.r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+			    i__4 = i__;
+			    q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, 
+				    q__2.i = q__3.r * x[i__4].i + q__3.i * x[
+				    i__4].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+/* L100: */
+			}
+		    }
+		    i__3 = j;
+		    x[i__3].r = temp.r, x[i__3].i = temp.i;
+/* L110: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    i__3 = jx;
+		    temp.r = x[i__3].r, temp.i = x[i__3].i;
+		    kx -= *incx;
+		    ix = kx;
+		    l = kplus1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__3 = kplus1 + j * a_dim1;
+			    q__1.r = temp.r * a[i__3].r - temp.i * a[i__3].i, 
+				    q__1.i = temp.r * a[i__3].i + temp.i * a[
+				    i__3].r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    i__4 = l + i__ + j * a_dim1;
+			    i__1 = ix;
+			    q__2.r = a[i__4].r * x[i__1].r - a[i__4].i * x[
+				    i__1].i, q__2.i = a[i__4].r * x[i__1].i + 
+				    a[i__4].i * x[i__1].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			    ix -= *incx;
+/* L120: */
+			}
+		    } else {
+			if (nounit) {
+			    r_cnjg(&q__2, &a[kplus1 + j * a_dim1]);
+			    q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
+				    q__1.i = temp.r * q__2.i + temp.i * 
+				    q__2.r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+			    i__4 = ix;
+			    q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, 
+				    q__2.i = q__3.r * x[i__4].i + q__3.i * x[
+				    i__4].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			    ix -= *incx;
+/* L130: */
+			}
+		    }
+		    i__3 = jx;
+		    x[i__3].r = temp.r, x[i__3].i = temp.i;
+		    jx -= *incx;
+/* L140: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    i__4 = j;
+		    temp.r = x[i__4].r, temp.i = x[i__4].i;
+		    l = 1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__4 = j * a_dim1 + 1;
+			    q__1.r = temp.r * a[i__4].r - temp.i * a[i__4].i, 
+				    q__1.i = temp.r * a[i__4].i + temp.i * a[
+				    i__4].r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    i__1 = l + i__ + j * a_dim1;
+			    i__2 = i__;
+			    q__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
+				    i__2].i, q__2.i = a[i__1].r * x[i__2].i + 
+				    a[i__1].i * x[i__2].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+/* L150: */
+			}
+		    } else {
+			if (nounit) {
+			    r_cnjg(&q__2, &a[j * a_dim1 + 1]);
+			    q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
+				    q__1.i = temp.r * q__2.i + temp.i * 
+				    q__2.r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+			    i__1 = i__;
+			    q__2.r = q__3.r * x[i__1].r - q__3.i * x[i__1].i, 
+				    q__2.i = q__3.r * x[i__1].i + q__3.i * x[
+				    i__1].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+/* L160: */
+			}
+		    }
+		    i__4 = j;
+		    x[i__4].r = temp.r, x[i__4].i = temp.i;
+/* L170: */
+		}
+	    } else {
+		jx = kx;
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    i__4 = jx;
+		    temp.r = x[i__4].r, temp.i = x[i__4].i;
+		    kx += *incx;
+		    ix = kx;
+		    l = 1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__4 = j * a_dim1 + 1;
+			    q__1.r = temp.r * a[i__4].r - temp.i * a[i__4].i, 
+				    q__1.i = temp.r * a[i__4].i + temp.i * a[
+				    i__4].r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    i__1 = l + i__ + j * a_dim1;
+			    i__2 = ix;
+			    q__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
+				    i__2].i, q__2.i = a[i__1].r * x[i__2].i + 
+				    a[i__1].i * x[i__2].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			    ix += *incx;
+/* L180: */
+			}
+		    } else {
+			if (nounit) {
+			    r_cnjg(&q__2, &a[j * a_dim1 + 1]);
+			    q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
+				    q__1.i = temp.r * q__2.i + temp.i * 
+				    q__2.r;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
+			    i__1 = ix;
+			    q__2.r = q__3.r * x[i__1].r - q__3.i * x[i__1].i, 
+				    q__2.i = q__3.r * x[i__1].i + q__3.i * x[
+				    i__1].r;
+			    q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
+				    q__2.i;
+			    temp.r = q__1.r, temp.i = q__1.i;
+			    ix += *incx;
+/* L190: */
+			}
+		    }
+		    i__4 = jx;
+		    x[i__4].r = temp.r, x[i__4].i = temp.i;
+		    jx += *incx;
+/* L200: */
+		}
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of CTBMV . */
+
+} /* ctbmv_ */
+
diff --git a/blas/f2c/d_cnjg.c b/blas/f2c/d_cnjg.c
new file mode 100644
index 000000000..623090c6b
--- /dev/null
+++ b/blas/f2c/d_cnjg.c
@@ -0,0 +1,6 @@
+#include "datatypes.h"    
+
+void d_cnjg(doublecomplex *r, doublecomplex *z) {
+    r->r = z->r;
+    r->i = -(z->i);
+}
diff --git a/blas/f2c/datatypes.h b/blas/f2c/datatypes.h
new file mode 100644
index 000000000..63232b246
--- /dev/null
+++ b/blas/f2c/datatypes.h
@@ -0,0 +1,24 @@
+/* This contains a limited subset of the typedefs exposed by f2c
+   for use by the Eigen BLAS C-only implementation.
+*/
+
+#ifndef __EIGEN_DATATYPES_H__
+#define __EIGEN_DATATYPES_H__
+
+typedef int integer;
+typedef unsigned int uinteger;
+typedef float real;
+typedef double doublereal;
+typedef struct { real r, i; } complex;
+typedef struct { doublereal r, i; } doublecomplex;
+typedef int ftnlen;
+typedef int logical;
+
+#define abs(x) ((x) >= 0 ? (x) : -(x))
+#define dabs(x) (doublereal)abs(x)
+#define min(a,b) ((a) <= (b) ? (a) : (b))
+#define max(a,b) ((a) >= (b) ? (a) : (b))
+#define dmin(a,b) (doublereal)min(a,b)
+#define dmax(a,b) (doublereal)max(a,b)
+
+#endif
diff --git a/blas/f2c/drotm.c b/blas/f2c/drotm.c
new file mode 100644
index 000000000..17a779b74
--- /dev/null
+++ b/blas/f2c/drotm.c
@@ -0,0 +1,215 @@
+/* drotm.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int drotm_(integer *n, doublereal *dx, integer *incx, 
+	doublereal *dy, integer *incy, doublereal *dparam)
+{
+    /* Initialized data */
+
+    static doublereal zero = 0.;
+    static doublereal two = 2.;
+
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer i__;
+    doublereal w, z__;
+    integer kx, ky;
+    doublereal dh11, dh12, dh21, dh22, dflag;
+    integer nsteps;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*     APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX */
+
+/*     (DX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN */
+/*     (DY**T) */
+
+/*     DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE */
+/*     LX = (-INCX)*N, AND SIMILARLY FOR SY USING LY AND INCY. */
+/*     WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. */
+
+/*     DFLAG=-1.D0     DFLAG=0.D0        DFLAG=1.D0     DFLAG=-2.D0 */
+
+/*       (DH11  DH12)    (1.D0  DH12)    (DH11  1.D0)    (1.D0  0.D0) */
+/*     H=(          )    (          )    (          )    (          ) */
+/*       (DH21  DH22),   (DH21  1.D0),   (-1.D0 DH22),   (0.D0  1.D0). */
+/*     SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  N      (input) INTEGER */
+/*         number of elements in input vector(s) */
+
+/*  DX     (input/output) DOUBLE PRECISION array, dimension N */
+/*         double precision vector with N elements */
+
+/*  INCX   (input) INTEGER */
+/*         storage spacing between elements of DX */
+
+/*  DY     (input/output) DOUBLE PRECISION array, dimension N */
+/*         double precision vector with N elements */
+
+/*  INCY   (input) INTEGER */
+/*         storage spacing between elements of DY */
+
+/*  DPARAM (input/output)  DOUBLE PRECISION array, dimension 5 */
+/*     DPARAM(1)=DFLAG */
+/*     DPARAM(2)=DH11 */
+/*     DPARAM(3)=DH21 */
+/*     DPARAM(4)=DH12 */
+/*     DPARAM(5)=DH22 */
+
+/*  ===================================================================== */
+
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. Data statements .. */
+    /* Parameter adjustments */
+    --dparam;
+    --dy;
+    --dx;
+
+    /* Function Body */
+/*     .. */
+
+    dflag = dparam[1];
+    if (*n <= 0 || dflag + two == zero) {
+	goto L140;
+    }
+    if (! (*incx == *incy && *incx > 0)) {
+	goto L70;
+    }
+
+    nsteps = *n * *incx;
+    if (dflag < 0.) {
+	goto L50;
+    } else if (dflag == 0) {
+	goto L10;
+    } else {
+	goto L30;
+    }
+L10:
+    dh12 = dparam[4];
+    dh21 = dparam[3];
+    i__1 = nsteps;
+    i__2 = *incx;
+    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+	w = dx[i__];
+	z__ = dy[i__];
+	dx[i__] = w + z__ * dh12;
+	dy[i__] = w * dh21 + z__;
+/* L20: */
+    }
+    goto L140;
+L30:
+    dh11 = dparam[2];
+    dh22 = dparam[5];
+    i__2 = nsteps;
+    i__1 = *incx;
+    for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
+	w = dx[i__];
+	z__ = dy[i__];
+	dx[i__] = w * dh11 + z__;
+	dy[i__] = -w + dh22 * z__;
+/* L40: */
+    }
+    goto L140;
+L50:
+    dh11 = dparam[2];
+    dh12 = dparam[4];
+    dh21 = dparam[3];
+    dh22 = dparam[5];
+    i__1 = nsteps;
+    i__2 = *incx;
+    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+	w = dx[i__];
+	z__ = dy[i__];
+	dx[i__] = w * dh11 + z__ * dh12;
+	dy[i__] = w * dh21 + z__ * dh22;
+/* L60: */
+    }
+    goto L140;
+L70:
+    kx = 1;
+    ky = 1;
+    if (*incx < 0) {
+	kx = (1 - *n) * *incx + 1;
+    }
+    if (*incy < 0) {
+	ky = (1 - *n) * *incy + 1;
+    }
+
+    if (dflag < 0.) {
+	goto L120;
+    } else if (dflag == 0) {
+	goto L80;
+    } else {
+	goto L100;
+    }
+L80:
+    dh12 = dparam[4];
+    dh21 = dparam[3];
+    i__2 = *n;
+    for (i__ = 1; i__ <= i__2; ++i__) {
+	w = dx[kx];
+	z__ = dy[ky];
+	dx[kx] = w + z__ * dh12;
+	dy[ky] = w * dh21 + z__;
+	kx += *incx;
+	ky += *incy;
+/* L90: */
+    }
+    goto L140;
+L100:
+    dh11 = dparam[2];
+    dh22 = dparam[5];
+    i__2 = *n;
+    for (i__ = 1; i__ <= i__2; ++i__) {
+	w = dx[kx];
+	z__ = dy[ky];
+	dx[kx] = w * dh11 + z__;
+	dy[ky] = -w + dh22 * z__;
+	kx += *incx;
+	ky += *incy;
+/* L110: */
+    }
+    goto L140;
+L120:
+    dh11 = dparam[2];
+    dh12 = dparam[4];
+    dh21 = dparam[3];
+    dh22 = dparam[5];
+    i__2 = *n;
+    for (i__ = 1; i__ <= i__2; ++i__) {
+	w = dx[kx];
+	z__ = dy[ky];
+	dx[kx] = w * dh11 + z__ * dh12;
+	dy[ky] = w * dh21 + z__ * dh22;
+	kx += *incx;
+	ky += *incy;
+/* L130: */
+    }
+L140:
+    return 0;
+} /* drotm_ */
+
diff --git a/blas/f2c/drotmg.c b/blas/f2c/drotmg.c
new file mode 100644
index 000000000..a63eb1083
--- /dev/null
+++ b/blas/f2c/drotmg.c
@@ -0,0 +1,293 @@
+/* drotmg.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int drotmg_(doublereal *dd1, doublereal *dd2, doublereal *
+	dx1, doublereal *dy1, doublereal *dparam)
+{
+    /* Initialized data */
+
+    static doublereal zero = 0.;
+    static doublereal one = 1.;
+    static doublereal two = 2.;
+    static doublereal gam = 4096.;
+    static doublereal gamsq = 16777216.;
+    static doublereal rgamsq = 5.9604645e-8;
+
+    /* Format strings */
+    static char fmt_120[] = "";
+    static char fmt_150[] = "";
+    static char fmt_180[] = "";
+    static char fmt_210[] = "";
+
+    /* System generated locals */
+    doublereal d__1;
+
+    /* Local variables */
+    doublereal du, dp1, dp2, dq1, dq2, dh11, dh12, dh21, dh22;
+    integer igo;
+    doublereal dflag, dtemp;
+
+    /* Assigned format variables */
+    static char *igo_fmt;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*     CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS */
+/*     THE SECOND COMPONENT OF THE 2-VECTOR  (DSQRT(DD1)*DX1,DSQRT(DD2)* */
+/*     DY2)**T. */
+/*     WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. */
+
+/*     DFLAG=-1.D0     DFLAG=0.D0        DFLAG=1.D0     DFLAG=-2.D0 */
+
+/*       (DH11  DH12)    (1.D0  DH12)    (DH11  1.D0)    (1.D0  0.D0) */
+/*     H=(          )    (          )    (          )    (          ) */
+/*       (DH21  DH22),   (DH21  1.D0),   (-1.D0 DH22),   (0.D0  1.D0). */
+/*     LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22 */
+/*     RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE */
+/*     VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.) */
+
+/*     THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE */
+/*     INEXACT.  THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE */
+/*     OF DD1 AND DD2.  ALL ACTUAL SCALING OF DATA IS DONE USING GAM. */
+
+
+/*  Arguments */
+/*  ========= */
+
+/*  DD1    (input/output) DOUBLE PRECISION */
+
+/*  DD2    (input/output) DOUBLE PRECISION */
+
+/*  DX1    (input/output) DOUBLE PRECISION */
+
+/*  DY1    (input) DOUBLE PRECISION */
+
+/*  DPARAM (input/output)  DOUBLE PRECISION array, dimension 5 */
+/*     DPARAM(1)=DFLAG */
+/*     DPARAM(2)=DH11 */
+/*     DPARAM(3)=DH21 */
+/*     DPARAM(4)=DH12 */
+/*     DPARAM(5)=DH22 */
+
+/*  ===================================================================== */
+
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Data statements .. */
+
+    /* Parameter adjustments */
+    --dparam;
+
+    /* Function Body */
+/*     .. */
+    if (! (*dd1 < zero)) {
+	goto L10;
+    }
+/*       GO ZERO-H-D-AND-DX1.. */
+    goto L60;
+L10:
+/*     CASE-DD1-NONNEGATIVE */
+    dp2 = *dd2 * *dy1;
+    if (! (dp2 == zero)) {
+	goto L20;
+    }
+    dflag = -two;
+    goto L260;
+/*     REGULAR-CASE.. */
+L20:
+    dp1 = *dd1 * *dx1;
+    dq2 = dp2 * *dy1;
+    dq1 = dp1 * *dx1;
+
+    if (! (abs(dq1) > abs(dq2))) {
+	goto L40;
+    }
+    dh21 = -(*dy1) / *dx1;
+    dh12 = dp2 / dp1;
+
+    du = one - dh12 * dh21;
+
+    if (! (du <= zero)) {
+	goto L30;
+    }
+/*         GO ZERO-H-D-AND-DX1.. */
+    goto L60;
+L30:
+    dflag = zero;
+    *dd1 /= du;
+    *dd2 /= du;
+    *dx1 *= du;
+/*         GO SCALE-CHECK.. */
+    goto L100;
+L40:
+    if (! (dq2 < zero)) {
+	goto L50;
+    }
+/*         GO ZERO-H-D-AND-DX1.. */
+    goto L60;
+L50:
+    dflag = one;
+    dh11 = dp1 / dp2;
+    dh22 = *dx1 / *dy1;
+    du = one + dh11 * dh22;
+    dtemp = *dd2 / du;
+    *dd2 = *dd1 / du;
+    *dd1 = dtemp;
+    *dx1 = *dy1 * du;
+/*         GO SCALE-CHECK */
+    goto L100;
+/*     PROCEDURE..ZERO-H-D-AND-DX1.. */
+L60:
+    dflag = -one;
+    dh11 = zero;
+    dh12 = zero;
+    dh21 = zero;
+    dh22 = zero;
+
+    *dd1 = zero;
+    *dd2 = zero;
+    *dx1 = zero;
+/*         RETURN.. */
+    goto L220;
+/*     PROCEDURE..FIX-H.. */
+L70:
+    if (! (dflag >= zero)) {
+	goto L90;
+    }
+
+    if (! (dflag == zero)) {
+	goto L80;
+    }
+    dh11 = one;
+    dh22 = one;
+    dflag = -one;
+    goto L90;
+L80:
+    dh21 = -one;
+    dh12 = one;
+    dflag = -one;
+L90:
+    switch (igo) {
+	case 0: goto L120;
+	case 1: goto L150;
+	case 2: goto L180;
+	case 3: goto L210;
+    }
+/*     PROCEDURE..SCALE-CHECK */
+L100:
+L110:
+    if (! (*dd1 <= rgamsq)) {
+	goto L130;
+    }
+    if (*dd1 == zero) {
+	goto L160;
+    }
+    igo = 0;
+    igo_fmt = fmt_120;
+/*              FIX-H.. */
+    goto L70;
+L120:
+/* Computing 2nd power */
+    d__1 = gam;
+    *dd1 *= d__1 * d__1;
+    *dx1 /= gam;
+    dh11 /= gam;
+    dh12 /= gam;
+    goto L110;
+L130:
+L140:
+    if (! (*dd1 >= gamsq)) {
+	goto L160;
+    }
+    igo = 1;
+    igo_fmt = fmt_150;
+/*              FIX-H.. */
+    goto L70;
+L150:
+/* Computing 2nd power */
+    d__1 = gam;
+    *dd1 /= d__1 * d__1;
+    *dx1 *= gam;
+    dh11 *= gam;
+    dh12 *= gam;
+    goto L140;
+L160:
+L170:
+    if (! (abs(*dd2) <= rgamsq)) {
+	goto L190;
+    }
+    if (*dd2 == zero) {
+	goto L220;
+    }
+    igo = 2;
+    igo_fmt = fmt_180;
+/*              FIX-H.. */
+    goto L70;
+L180:
+/* Computing 2nd power */
+    d__1 = gam;
+    *dd2 *= d__1 * d__1;
+    dh21 /= gam;
+    dh22 /= gam;
+    goto L170;
+L190:
+L200:
+    if (! (abs(*dd2) >= gamsq)) {
+	goto L220;
+    }
+    igo = 3;
+    igo_fmt = fmt_210;
+/*              FIX-H.. */
+    goto L70;
+L210:
+/* Computing 2nd power */
+    d__1 = gam;
+    *dd2 /= d__1 * d__1;
+    dh21 *= gam;
+    dh22 *= gam;
+    goto L200;
+L220:
+    if (dflag < 0.) {
+	goto L250;
+    } else if (dflag == 0) {
+	goto L230;
+    } else {
+	goto L240;
+    }
+L230:
+    dparam[3] = dh21;
+    dparam[4] = dh12;
+    goto L260;
+L240:
+    dparam[2] = dh11;
+    dparam[5] = dh22;
+    goto L260;
+L250:
+    dparam[2] = dh11;
+    dparam[3] = dh21;
+    dparam[4] = dh12;
+    dparam[5] = dh22;
+L260:
+    dparam[1] = dflag;
+    return 0;
+} /* drotmg_ */
+
diff --git a/blas/f2c/dsbmv.c b/blas/f2c/dsbmv.c
new file mode 100644
index 000000000..c6b4b21d6
--- /dev/null
+++ b/blas/f2c/dsbmv.c
@@ -0,0 +1,366 @@
+/* dsbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int dsbmv_(char *uplo, integer *n, integer *k, doublereal *
+	alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, 
+	doublereal *beta, doublereal *y, integer *incy, ftnlen uplo_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
+    doublereal temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  DSBMV  performs the matrix-vector  operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n symmetric band matrix, with k super-diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the band matrix A is being supplied as */
+/*           follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  being supplied. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  being supplied. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry, K specifies the number of super-diagonals of the */
+/*           matrix A. K must satisfy  0 .le. K. */
+/*           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 with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the symmetric matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer the upper */
+/*           triangular part of a symmetric band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the symmetric matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer the lower */
+/*           triangular part of a symmetric band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + 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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - DOUBLE PRECISION array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           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. */
+/*           Unchanged on exit. */
+
+/*  Y      - DOUBLE PRECISION array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCY ) ). */
+/*           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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+    --y;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*k < 0) {
+	info = 3;
+    } else if (*lda < *k + 1) {
+	info = 6;
+    } else if (*incx == 0) {
+	info = 8;
+    } else if (*incy == 0) {
+	info = 11;
+    }
+    if (info != 0) {
+	xerbla_("DSBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (*alpha == 0. && *beta == 1.)) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array A */
+/*     are accessed sequentially with one pass through A. */
+
+/*     First form  y := beta*y. */
+
+    if (*beta != 1.) {
+	if (*incy == 1) {
+	    if (*beta == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = 0.;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = *beta * y[i__];
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (*beta == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = 0.;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = *beta * y[iy];
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (*alpha == 0.) {
+	return 0;
+    }
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when upper triangle of A is stored. */
+
+	kplus1 = *k + 1;
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__2 = 1, i__3 = j - *k;
+		i__4 = j - 1;
+		for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+		    y[i__] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[i__];
+/* L50: */
+		}
+		y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.;
+		ix = kx;
+		iy = ky;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__4 = 1, i__2 = j - *k;
+		i__3 = j - 1;
+		for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+		    y[iy] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[ix];
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha * 
+			temp2;
+		jx += *incx;
+		jy += *incy;
+		if (j > *k) {
+		    kx += *incx;
+		    ky += *incy;
+		}
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when lower triangle of A is stored. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.;
+		y[j] += temp1 * a[j * a_dim1 + 1];
+		l = 1 - j;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    y[i__] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[i__];
+/* L90: */
+		}
+		y[j] += *alpha * temp2;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.;
+		y[jy] += temp1 * a[j * a_dim1 + 1];
+		l = 1 - j;
+		ix = jx;
+		iy = jy;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    ix += *incx;
+		    iy += *incy;
+		    y[iy] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[ix];
+/* L110: */
+		}
+		y[jy] += *alpha * temp2;
+		jx += *incx;
+		jy += *incy;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of DSBMV . */
+
+} /* dsbmv_ */
+
diff --git a/blas/f2c/dspmv.c b/blas/f2c/dspmv.c
new file mode 100644
index 000000000..0b4e92d5c
--- /dev/null
+++ b/blas/f2c/dspmv.c
@@ -0,0 +1,316 @@
+/* dspmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int dspmv_(char *uplo, integer *n, doublereal *alpha, 
+	doublereal *ap, doublereal *x, integer *incx, doublereal *beta, 
+	doublereal *y, integer *incy, ftnlen uplo_len)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
+    doublereal temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  DSPMV  performs the matrix-vector operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n symmetric matrix, supplied in packed form. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the matrix A is supplied in the packed */
+/*           array AP as follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  supplied in AP. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  supplied in AP. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  ALPHA  - DOUBLE PRECISION. */
+/*           On entry, ALPHA specifies the scalar alpha. */
+/*           Unchanged on exit. */
+
+/*  AP     - DOUBLE PRECISION array of DIMENSION at least */
+/*           ( ( n*( n + 1 ) )/2 ). */
+/*           Before entry with UPLO = 'U' or 'u', the array AP must */
+/*           contain the upper triangular part of the symmetric matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
+/*           and a( 2, 2 ) respectively, and so on. */
+/*           Before entry with UPLO = 'L' or 'l', the array AP must */
+/*           contain the lower triangular part of the symmetric matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
+/*           and a( 3, 1 ) respectively, and so on. */
+/*           Unchanged on exit. */
+
+/*  X      - DOUBLE PRECISION array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element 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 + ( n - 1 )*abs( INCY ) ). */
+/*           Before entry, the incremented array Y must contain the n */
+/*           element 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --y;
+    --x;
+    --ap;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*incx == 0) {
+	info = 6;
+    } else if (*incy == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("DSPMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (*alpha == 0. && *beta == 1.)) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+/*     First form  y := beta*y. */
+
+    if (*beta != 1.) {
+	if (*incy == 1) {
+	    if (*beta == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = 0.;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = *beta * y[i__];
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (*beta == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = 0.;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = *beta * y[iy];
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (*alpha == 0.) {
+	return 0;
+    }
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when AP contains the upper triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.;
+		k = kk;
+		i__2 = j - 1;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    y[i__] += temp1 * ap[k];
+		    temp2 += ap[k] * x[i__];
+		    ++k;
+/* L50: */
+		}
+		y[j] = y[j] + temp1 * ap[kk + j - 1] + *alpha * temp2;
+		kk += j;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.;
+		ix = kx;
+		iy = ky;
+		i__2 = kk + j - 2;
+		for (k = kk; k <= i__2; ++k) {
+		    y[iy] += temp1 * ap[k];
+		    temp2 += ap[k] * x[ix];
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		y[jy] = y[jy] + temp1 * ap[kk + j - 1] + *alpha * temp2;
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when AP contains the lower triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.;
+		y[j] += temp1 * ap[kk];
+		k = kk + 1;
+		i__2 = *n;
+		for (i__ = j + 1; i__ <= i__2; ++i__) {
+		    y[i__] += temp1 * ap[k];
+		    temp2 += ap[k] * x[i__];
+		    ++k;
+/* L90: */
+		}
+		y[j] += *alpha * temp2;
+		kk += *n - j + 1;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.;
+		y[jy] += temp1 * ap[kk];
+		ix = jx;
+		iy = jy;
+		i__2 = kk + *n - j;
+		for (k = kk + 1; k <= i__2; ++k) {
+		    ix += *incx;
+		    iy += *incy;
+		    y[iy] += temp1 * ap[k];
+		    temp2 += ap[k] * x[ix];
+/* L110: */
+		}
+		y[jy] += *alpha * temp2;
+		jx += *incx;
+		jy += *incy;
+		kk += *n - j + 1;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of DSPMV . */
+
+} /* dspmv_ */
+
diff --git a/blas/f2c/dtbmv.c b/blas/f2c/dtbmv.c
new file mode 100644
index 000000000..fdf73ebb5
--- /dev/null
+++ b/blas/f2c/dtbmv.c
@@ -0,0 +1,428 @@
+/* dtbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int dtbmv_(char *uplo, char *trans, char *diag, integer *n, 
+	integer *k, doublereal *a, integer *lda, doublereal *x, integer *incx,
+	 ftnlen uplo_len, ftnlen trans_len, ftnlen diag_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    integer i__, j, l, ix, jx, kx, info;
+    doublereal temp;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+    logical nounit;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  DTBMV  performs one of the matrix-vector operations */
+
+/*     x := A*x,   or   x := A'*x, */
+
+/*  where x is an n element vector and  A is an n by n unit, or non-unit, */
+/*  upper or lower triangular band matrix, with ( k + 1 ) diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the matrix is an upper or */
+/*           lower triangular matrix as follows: */
+
+/*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
+
+/*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
+
+/*           Unchanged on exit. */
+
+/*  TRANS  - CHARACTER*1. */
+/*           On entry, TRANS specifies the operation to be performed as */
+/*           follows: */
+
+/*              TRANS = 'N' or 'n'   x := A*x. */
+
+/*              TRANS = 'T' or 't'   x := A'*x. */
+
+/*              TRANS = 'C' or 'c'   x := A'*x. */
+
+/*           Unchanged on exit. */
+
+/*  DIAG   - CHARACTER*1. */
+/*           On entry, DIAG specifies whether or not A is unit */
+/*           triangular as follows: */
+
+/*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
+
+/*              DIAG = 'N' or 'n'   A is not assumed to be unit */
+/*                                  triangular. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry with UPLO = 'U' or 'u', K specifies the number of */
+/*           super-diagonals of the matrix A. */
+/*           On entry with UPLO = 'L' or 'l', K specifies the number of */
+/*           sub-diagonals of the matrix A. */
+/*           K must satisfy  0 .le. K. */
+/*           Unchanged on exit. */
+
+/*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
+/*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer an upper */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer a lower */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Note that when DIAG = 'U' or 'u' the elements of the array A */
+/*           corresponding to the diagonal elements of the matrix are not */
+/*           referenced, but are assumed to be unity. */
+/*           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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - DOUBLE PRECISION array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element vector x. On exit, X is overwritten with the */
+/*           tranformed vector x. */
+
+/*  INCX   - INTEGER. */
+/*           On entry, INCX specifies the increment for the elements of */
+/*           X. INCX 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, 
+	    "T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
+	    ftnlen)1)) {
+	info = 2;
+    } else if (! lsame_(diag, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(diag, 
+	    "N", (ftnlen)1, (ftnlen)1)) {
+	info = 3;
+    } else if (*n < 0) {
+	info = 4;
+    } else if (*k < 0) {
+	info = 5;
+    } else if (*lda < *k + 1) {
+	info = 7;
+    } else if (*incx == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("DTBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0) {
+	return 0;
+    }
+
+    nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
+
+/*     Set up the start point in X if the increment is not unity. This */
+/*     will be  ( N - 1 )*INCX   too small for descending loops. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of A are */
+/*     accessed sequentially with one pass through A. */
+
+    if (lsame_(trans, "N", (ftnlen)1, (ftnlen)1)) {
+
+/*         Form  x := A*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    if (x[j] != 0.) {
+			temp = x[j];
+			l = kplus1 - j;
+/* Computing MAX */
+			i__2 = 1, i__3 = j - *k;
+			i__4 = j - 1;
+			for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+			    x[i__] += temp * a[l + i__ + j * a_dim1];
+/* L10: */
+			}
+			if (nounit) {
+			    x[j] *= a[kplus1 + j * a_dim1];
+			}
+		    }
+/* L20: */
+		}
+	    } else {
+		jx = kx;
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    if (x[jx] != 0.) {
+			temp = x[jx];
+			ix = kx;
+			l = kplus1 - j;
+/* Computing MAX */
+			i__4 = 1, i__2 = j - *k;
+			i__3 = j - 1;
+			for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+			    x[ix] += temp * a[l + i__ + j * a_dim1];
+			    ix += *incx;
+/* L30: */
+			}
+			if (nounit) {
+			    x[jx] *= a[kplus1 + j * a_dim1];
+			}
+		    }
+		    jx += *incx;
+		    if (j > *k) {
+			kx += *incx;
+		    }
+/* L40: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    if (x[j] != 0.) {
+			temp = x[j];
+			l = 1 - j;
+/* Computing MIN */
+			i__1 = *n, i__3 = j + *k;
+			i__4 = j + 1;
+			for (i__ = min(i__1,i__3); i__ >= i__4; --i__) {
+			    x[i__] += temp * a[l + i__ + j * a_dim1];
+/* L50: */
+			}
+			if (nounit) {
+			    x[j] *= a[j * a_dim1 + 1];
+			}
+		    }
+/* L60: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    if (x[jx] != 0.) {
+			temp = x[jx];
+			ix = kx;
+			l = 1 - j;
+/* Computing MIN */
+			i__4 = *n, i__1 = j + *k;
+			i__3 = j + 1;
+			for (i__ = min(i__4,i__1); i__ >= i__3; --i__) {
+			    x[ix] += temp * a[l + i__ + j * a_dim1];
+			    ix -= *incx;
+/* L70: */
+			}
+			if (nounit) {
+			    x[jx] *= a[j * a_dim1 + 1];
+			}
+		    }
+		    jx -= *incx;
+		    if (*n - j >= *k) {
+			kx -= *incx;
+		    }
+/* L80: */
+		}
+	    }
+	}
+    } else {
+
+/*        Form  x := A'*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    temp = x[j];
+		    l = kplus1 - j;
+		    if (nounit) {
+			temp *= a[kplus1 + j * a_dim1];
+		    }
+/* Computing MAX */
+		    i__4 = 1, i__1 = j - *k;
+		    i__3 = max(i__4,i__1);
+		    for (i__ = j - 1; i__ >= i__3; --i__) {
+			temp += a[l + i__ + j * a_dim1] * x[i__];
+/* L90: */
+		    }
+		    x[j] = temp;
+/* L100: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    temp = x[jx];
+		    kx -= *incx;
+		    ix = kx;
+		    l = kplus1 - j;
+		    if (nounit) {
+			temp *= a[kplus1 + j * a_dim1];
+		    }
+/* Computing MAX */
+		    i__4 = 1, i__1 = j - *k;
+		    i__3 = max(i__4,i__1);
+		    for (i__ = j - 1; i__ >= i__3; --i__) {
+			temp += a[l + i__ + j * a_dim1] * x[ix];
+			ix -= *incx;
+/* L110: */
+		    }
+		    x[jx] = temp;
+		    jx -= *incx;
+/* L120: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    temp = x[j];
+		    l = 1 - j;
+		    if (nounit) {
+			temp *= a[j * a_dim1 + 1];
+		    }
+/* Computing MIN */
+		    i__1 = *n, i__2 = j + *k;
+		    i__4 = min(i__1,i__2);
+		    for (i__ = j + 1; i__ <= i__4; ++i__) {
+			temp += a[l + i__ + j * a_dim1] * x[i__];
+/* L130: */
+		    }
+		    x[j] = temp;
+/* L140: */
+		}
+	    } else {
+		jx = kx;
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    temp = x[jx];
+		    kx += *incx;
+		    ix = kx;
+		    l = 1 - j;
+		    if (nounit) {
+			temp *= a[j * a_dim1 + 1];
+		    }
+/* Computing MIN */
+		    i__1 = *n, i__2 = j + *k;
+		    i__4 = min(i__1,i__2);
+		    for (i__ = j + 1; i__ <= i__4; ++i__) {
+			temp += a[l + i__ + j * a_dim1] * x[ix];
+			ix += *incx;
+/* L150: */
+		    }
+		    x[jx] = temp;
+		    jx += *incx;
+/* L160: */
+		}
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of DTBMV . */
+
+} /* dtbmv_ */
+
diff --git a/blas/f2c/lsame.c b/blas/f2c/lsame.c
new file mode 100644
index 000000000..46324d916
--- /dev/null
+++ b/blas/f2c/lsame.c
@@ -0,0 +1,117 @@
+/* lsame.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+logical lsame_(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
+{
+    /* System generated locals */
+    logical ret_val;
+
+    /* Local variables */
+    integer inta, intb, zcode;
+
+
+/*  -- LAPACK auxiliary routine (version 3.1) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/*     November 2006 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  LSAME returns .TRUE. if CA is the same letter as CB regardless of */
+/*  case. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  CA      (input) CHARACTER*1 */
+
+/*  CB      (input) CHARACTER*1 */
+/*          CA and CB specify the single characters to be compared. */
+
+/* ===================================================================== */
+
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+
+/*     Test if the characters are equal */
+
+    ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
+    if (ret_val) {
+	return ret_val;
+    }
+
+/*     Now test for equivalence if both characters are alphabetic. */
+
+    zcode = 'Z';
+
+/*     Use 'Z' rather than 'A' so that ASCII can be detected on Prime */
+/*     machines, on which ICHAR returns a value with bit 8 set. */
+/*     ICHAR('A') on Prime machines returns 193 which is the same as */
+/*     ICHAR('A') on an EBCDIC machine. */
+
+    inta = *(unsigned char *)ca;
+    intb = *(unsigned char *)cb;
+
+    if (zcode == 90 || zcode == 122) {
+
+/*        ASCII is assumed - ZCODE is the ASCII code of either lower or */
+/*        upper case 'Z'. */
+
+	if (inta >= 97 && inta <= 122) {
+	    inta += -32;
+	}
+	if (intb >= 97 && intb <= 122) {
+	    intb += -32;
+	}
+
+    } else if (zcode == 233 || zcode == 169) {
+
+/*        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or */
+/*        upper case 'Z'. */
+
+	if ((inta >= 129 && inta <= 137) || (inta >= 145 && inta <= 153) || 
+            (inta >= 162 && inta <= 169)) {
+	    inta += 64;
+	}
+	if ((intb >= 129 && intb <= 137) || (intb >= 145 && intb <= 153) || 
+            (intb >= 162 && intb <= 169)) {
+	    intb += 64;
+	}
+
+    } else if (zcode == 218 || zcode == 250) {
+
+/*        ASCII is assumed, on Prime machines - ZCODE is the ASCII code */
+/*        plus 128 of either lower or upper case 'Z'. */
+
+	if (inta >= 225 && inta <= 250) {
+	    inta += -32;
+	}
+	if (intb >= 225 && intb <= 250) {
+	    intb += -32;
+	}
+    }
+    ret_val = inta == intb;
+
+/*     RETURN */
+
+/*     End of LSAME */
+
+    return ret_val;
+} /* lsame_ */
+
diff --git a/blas/f2c/r_cnjg.c b/blas/f2c/r_cnjg.c
new file mode 100644
index 000000000..c08182f88
--- /dev/null
+++ b/blas/f2c/r_cnjg.c
@@ -0,0 +1,6 @@
+#include "datatypes.h"    
+
+void r_cnjg(complex *r, complex *z) {
+    r->r = z->r;
+    r->i = -(z->i);
+}
diff --git a/blas/f2c/srotm.c b/blas/f2c/srotm.c
new file mode 100644
index 000000000..bd5944a99
--- /dev/null
+++ b/blas/f2c/srotm.c
@@ -0,0 +1,216 @@
+/* srotm.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int srotm_(integer *n, real *sx, integer *incx, real *sy, 
+	integer *incy, real *sparam)
+{
+    /* Initialized data */
+
+    static real zero = 0.f;
+    static real two = 2.f;
+
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer i__;
+    real w, z__;
+    integer kx, ky;
+    real sh11, sh12, sh21, sh22, sflag;
+    integer nsteps;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*     APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX */
+
+/*     (SX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF SX ARE IN */
+/*     (DX**T) */
+
+/*     SX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE */
+/*     LX = (-INCX)*N, AND SIMILARLY FOR SY USING USING LY AND INCY. */
+/*     WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS.. */
+
+/*     SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0 */
+
+/*       (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0) */
+/*     H=(          )    (          )    (          )    (          ) */
+/*       (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0). */
+/*     SEE  SROTMG FOR A DESCRIPTION OF DATA STORAGE IN SPARAM. */
+
+
+/*  Arguments */
+/*  ========= */
+
+/*  N      (input) INTEGER */
+/*         number of elements in input vector(s) */
+
+/*  SX     (input/output) REAL array, dimension N */
+/*         double precision vector with N elements */
+
+/*  INCX   (input) INTEGER */
+/*         storage spacing between elements of SX */
+
+/*  SY     (input/output) REAL array, dimension N */
+/*         double precision vector with N elements */
+
+/*  INCY   (input) INTEGER */
+/*         storage spacing between elements of SY */
+
+/*  SPARAM (input/output)  REAL array, dimension 5 */
+/*     SPARAM(1)=SFLAG */
+/*     SPARAM(2)=SH11 */
+/*     SPARAM(3)=SH21 */
+/*     SPARAM(4)=SH12 */
+/*     SPARAM(5)=SH22 */
+
+/*  ===================================================================== */
+
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. Data statements .. */
+    /* Parameter adjustments */
+    --sparam;
+    --sy;
+    --sx;
+
+    /* Function Body */
+/*     .. */
+
+    sflag = sparam[1];
+    if (*n <= 0 || sflag + two == zero) {
+	goto L140;
+    }
+    if (! (*incx == *incy && *incx > 0)) {
+	goto L70;
+    }
+
+    nsteps = *n * *incx;
+    if (sflag < 0.f) {
+	goto L50;
+    } else if (sflag == 0) {
+	goto L10;
+    } else {
+	goto L30;
+    }
+L10:
+    sh12 = sparam[4];
+    sh21 = sparam[3];
+    i__1 = nsteps;
+    i__2 = *incx;
+    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+	w = sx[i__];
+	z__ = sy[i__];
+	sx[i__] = w + z__ * sh12;
+	sy[i__] = w * sh21 + z__;
+/* L20: */
+    }
+    goto L140;
+L30:
+    sh11 = sparam[2];
+    sh22 = sparam[5];
+    i__2 = nsteps;
+    i__1 = *incx;
+    for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
+	w = sx[i__];
+	z__ = sy[i__];
+	sx[i__] = w * sh11 + z__;
+	sy[i__] = -w + sh22 * z__;
+/* L40: */
+    }
+    goto L140;
+L50:
+    sh11 = sparam[2];
+    sh12 = sparam[4];
+    sh21 = sparam[3];
+    sh22 = sparam[5];
+    i__1 = nsteps;
+    i__2 = *incx;
+    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+	w = sx[i__];
+	z__ = sy[i__];
+	sx[i__] = w * sh11 + z__ * sh12;
+	sy[i__] = w * sh21 + z__ * sh22;
+/* L60: */
+    }
+    goto L140;
+L70:
+    kx = 1;
+    ky = 1;
+    if (*incx < 0) {
+	kx = (1 - *n) * *incx + 1;
+    }
+    if (*incy < 0) {
+	ky = (1 - *n) * *incy + 1;
+    }
+
+    if (sflag < 0.f) {
+	goto L120;
+    } else if (sflag == 0) {
+	goto L80;
+    } else {
+	goto L100;
+    }
+L80:
+    sh12 = sparam[4];
+    sh21 = sparam[3];
+    i__2 = *n;
+    for (i__ = 1; i__ <= i__2; ++i__) {
+	w = sx[kx];
+	z__ = sy[ky];
+	sx[kx] = w + z__ * sh12;
+	sy[ky] = w * sh21 + z__;
+	kx += *incx;
+	ky += *incy;
+/* L90: */
+    }
+    goto L140;
+L100:
+    sh11 = sparam[2];
+    sh22 = sparam[5];
+    i__2 = *n;
+    for (i__ = 1; i__ <= i__2; ++i__) {
+	w = sx[kx];
+	z__ = sy[ky];
+	sx[kx] = w * sh11 + z__;
+	sy[ky] = -w + sh22 * z__;
+	kx += *incx;
+	ky += *incy;
+/* L110: */
+    }
+    goto L140;
+L120:
+    sh11 = sparam[2];
+    sh12 = sparam[4];
+    sh21 = sparam[3];
+    sh22 = sparam[5];
+    i__2 = *n;
+    for (i__ = 1; i__ <= i__2; ++i__) {
+	w = sx[kx];
+	z__ = sy[ky];
+	sx[kx] = w * sh11 + z__ * sh12;
+	sy[ky] = w * sh21 + z__ * sh22;
+	kx += *incx;
+	ky += *incy;
+/* L130: */
+    }
+L140:
+    return 0;
+} /* srotm_ */
+
diff --git a/blas/f2c/srotmg.c b/blas/f2c/srotmg.c
new file mode 100644
index 000000000..75f789fe2
--- /dev/null
+++ b/blas/f2c/srotmg.c
@@ -0,0 +1,295 @@
+/* srotmg.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int srotmg_(real *sd1, real *sd2, real *sx1, real *sy1, real 
+	*sparam)
+{
+    /* Initialized data */
+
+    static real zero = 0.f;
+    static real one = 1.f;
+    static real two = 2.f;
+    static real gam = 4096.f;
+    static real gamsq = 16777200.f;
+    static real rgamsq = 5.96046e-8f;
+
+    /* Format strings */
+    static char fmt_120[] = "";
+    static char fmt_150[] = "";
+    static char fmt_180[] = "";
+    static char fmt_210[] = "";
+
+    /* System generated locals */
+    real r__1;
+
+    /* Local variables */
+    real su, sp1, sp2, sq1, sq2, sh11, sh12, sh21, sh22;
+    integer igo;
+    real sflag, stemp;
+
+    /* Assigned format variables */
+    static char *igo_fmt;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*     CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS */
+/*     THE SECOND COMPONENT OF THE 2-VECTOR  (SQRT(SD1)*SX1,SQRT(SD2)* */
+/*     SY2)**T. */
+/*     WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS.. */
+
+/*     SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0 */
+
+/*       (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0) */
+/*     H=(          )    (          )    (          )    (          ) */
+/*       (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0). */
+/*     LOCATIONS 2-4 OF SPARAM CONTAIN SH11,SH21,SH12, AND SH22 */
+/*     RESPECTIVELY. (VALUES OF 1.E0, -1.E0, OR 0.E0 IMPLIED BY THE */
+/*     VALUE OF SPARAM(1) ARE NOT STORED IN SPARAM.) */
+
+/*     THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE */
+/*     INEXACT.  THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE */
+/*     OF SD1 AND SD2.  ALL ACTUAL SCALING OF DATA IS DONE USING GAM. */
+
+
+/*  Arguments */
+/*  ========= */
+
+
+/*  SD1    (input/output) REAL */
+
+/*  SD2    (input/output) REAL */
+
+/*  SX1    (input/output) REAL */
+
+/*  SY1    (input) REAL */
+
+
+/*  SPARAM (input/output)  REAL array, dimension 5 */
+/*     SPARAM(1)=SFLAG */
+/*     SPARAM(2)=SH11 */
+/*     SPARAM(3)=SH21 */
+/*     SPARAM(4)=SH12 */
+/*     SPARAM(5)=SH22 */
+
+/*  ===================================================================== */
+
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Data statements .. */
+
+    /* Parameter adjustments */
+    --sparam;
+
+    /* Function Body */
+/*     .. */
+    if (! (*sd1 < zero)) {
+	goto L10;
+    }
+/*       GO ZERO-H-D-AND-SX1.. */
+    goto L60;
+L10:
+/*     CASE-SD1-NONNEGATIVE */
+    sp2 = *sd2 * *sy1;
+    if (! (sp2 == zero)) {
+	goto L20;
+    }
+    sflag = -two;
+    goto L260;
+/*     REGULAR-CASE.. */
+L20:
+    sp1 = *sd1 * *sx1;
+    sq2 = sp2 * *sy1;
+    sq1 = sp1 * *sx1;
+
+    if (! (dabs(sq1) > dabs(sq2))) {
+	goto L40;
+    }
+    sh21 = -(*sy1) / *sx1;
+    sh12 = sp2 / sp1;
+
+    su = one - sh12 * sh21;
+
+    if (! (su <= zero)) {
+	goto L30;
+    }
+/*         GO ZERO-H-D-AND-SX1.. */
+    goto L60;
+L30:
+    sflag = zero;
+    *sd1 /= su;
+    *sd2 /= su;
+    *sx1 *= su;
+/*         GO SCALE-CHECK.. */
+    goto L100;
+L40:
+    if (! (sq2 < zero)) {
+	goto L50;
+    }
+/*         GO ZERO-H-D-AND-SX1.. */
+    goto L60;
+L50:
+    sflag = one;
+    sh11 = sp1 / sp2;
+    sh22 = *sx1 / *sy1;
+    su = one + sh11 * sh22;
+    stemp = *sd2 / su;
+    *sd2 = *sd1 / su;
+    *sd1 = stemp;
+    *sx1 = *sy1 * su;
+/*         GO SCALE-CHECK */
+    goto L100;
+/*     PROCEDURE..ZERO-H-D-AND-SX1.. */
+L60:
+    sflag = -one;
+    sh11 = zero;
+    sh12 = zero;
+    sh21 = zero;
+    sh22 = zero;
+
+    *sd1 = zero;
+    *sd2 = zero;
+    *sx1 = zero;
+/*         RETURN.. */
+    goto L220;
+/*     PROCEDURE..FIX-H.. */
+L70:
+    if (! (sflag >= zero)) {
+	goto L90;
+    }
+
+    if (! (sflag == zero)) {
+	goto L80;
+    }
+    sh11 = one;
+    sh22 = one;
+    sflag = -one;
+    goto L90;
+L80:
+    sh21 = -one;
+    sh12 = one;
+    sflag = -one;
+L90:
+    switch (igo) {
+	case 0: goto L120;
+	case 1: goto L150;
+	case 2: goto L180;
+	case 3: goto L210;
+    }
+/*     PROCEDURE..SCALE-CHECK */
+L100:
+L110:
+    if (! (*sd1 <= rgamsq)) {
+	goto L130;
+    }
+    if (*sd1 == zero) {
+	goto L160;
+    }
+    igo = 0;
+    igo_fmt = fmt_120;
+/*              FIX-H.. */
+    goto L70;
+L120:
+/* Computing 2nd power */
+    r__1 = gam;
+    *sd1 *= r__1 * r__1;
+    *sx1 /= gam;
+    sh11 /= gam;
+    sh12 /= gam;
+    goto L110;
+L130:
+L140:
+    if (! (*sd1 >= gamsq)) {
+	goto L160;
+    }
+    igo = 1;
+    igo_fmt = fmt_150;
+/*              FIX-H.. */
+    goto L70;
+L150:
+/* Computing 2nd power */
+    r__1 = gam;
+    *sd1 /= r__1 * r__1;
+    *sx1 *= gam;
+    sh11 *= gam;
+    sh12 *= gam;
+    goto L140;
+L160:
+L170:
+    if (! (dabs(*sd2) <= rgamsq)) {
+	goto L190;
+    }
+    if (*sd2 == zero) {
+	goto L220;
+    }
+    igo = 2;
+    igo_fmt = fmt_180;
+/*              FIX-H.. */
+    goto L70;
+L180:
+/* Computing 2nd power */
+    r__1 = gam;
+    *sd2 *= r__1 * r__1;
+    sh21 /= gam;
+    sh22 /= gam;
+    goto L170;
+L190:
+L200:
+    if (! (dabs(*sd2) >= gamsq)) {
+	goto L220;
+    }
+    igo = 3;
+    igo_fmt = fmt_210;
+/*              FIX-H.. */
+    goto L70;
+L210:
+/* Computing 2nd power */
+    r__1 = gam;
+    *sd2 /= r__1 * r__1;
+    sh21 *= gam;
+    sh22 *= gam;
+    goto L200;
+L220:
+    if (sflag < 0.f) {
+	goto L250;
+    } else if (sflag == 0) {
+	goto L230;
+    } else {
+	goto L240;
+    }
+L230:
+    sparam[3] = sh21;
+    sparam[4] = sh12;
+    goto L260;
+L240:
+    sparam[2] = sh11;
+    sparam[5] = sh22;
+    goto L260;
+L250:
+    sparam[2] = sh11;
+    sparam[3] = sh21;
+    sparam[4] = sh12;
+    sparam[5] = sh22;
+L260:
+    sparam[1] = sflag;
+    return 0;
+} /* srotmg_ */
+
diff --git a/blas/f2c/ssbmv.c b/blas/f2c/ssbmv.c
new file mode 100644
index 000000000..8599325f2
--- /dev/null
+++ b/blas/f2c/ssbmv.c
@@ -0,0 +1,368 @@
+/* ssbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int ssbmv_(char *uplo, integer *n, integer *k, real *alpha, 
+	real *a, integer *lda, real *x, integer *incx, real *beta, real *y, 
+	integer *incy, ftnlen uplo_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
+    real temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  SSBMV  performs the matrix-vector  operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n symmetric band matrix, with k super-diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the band matrix A is being supplied as */
+/*           follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  being supplied. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  being supplied. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry, K specifies the number of super-diagonals of the */
+/*           matrix A. K must satisfy  0 .le. K. */
+/*           Unchanged on exit. */
+
+/*  ALPHA  - REAL            . */
+/*           On entry, ALPHA specifies the scalar alpha. */
+/*           Unchanged on exit. */
+
+/*  A      - REAL             array of DIMENSION ( LDA, n ). */
+/*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the symmetric matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer the upper */
+/*           triangular part of a symmetric band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the symmetric matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer the lower */
+/*           triangular part of a symmetric band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + 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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - REAL             array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           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. */
+/*           Unchanged on exit. */
+
+/*  Y      - REAL             array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCY ) ). */
+/*           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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+    --y;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*k < 0) {
+	info = 3;
+    } else if (*lda < *k + 1) {
+	info = 6;
+    } else if (*incx == 0) {
+	info = 8;
+    } else if (*incy == 0) {
+	info = 11;
+    }
+    if (info != 0) {
+	xerbla_("SSBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (*alpha == 0.f && *beta == 1.f)) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array A */
+/*     are accessed sequentially with one pass through A. */
+
+/*     First form  y := beta*y. */
+
+    if (*beta != 1.f) {
+	if (*incy == 1) {
+	    if (*beta == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = 0.f;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = *beta * y[i__];
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (*beta == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = 0.f;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = *beta * y[iy];
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (*alpha == 0.f) {
+	return 0;
+    }
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when upper triangle of A is stored. */
+
+	kplus1 = *k + 1;
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.f;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__2 = 1, i__3 = j - *k;
+		i__4 = j - 1;
+		for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+		    y[i__] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[i__];
+/* L50: */
+		}
+		y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.f;
+		ix = kx;
+		iy = ky;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__4 = 1, i__2 = j - *k;
+		i__3 = j - 1;
+		for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+		    y[iy] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[ix];
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha * 
+			temp2;
+		jx += *incx;
+		jy += *incy;
+		if (j > *k) {
+		    kx += *incx;
+		    ky += *incy;
+		}
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when lower triangle of A is stored. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.f;
+		y[j] += temp1 * a[j * a_dim1 + 1];
+		l = 1 - j;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    y[i__] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[i__];
+/* L90: */
+		}
+		y[j] += *alpha * temp2;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.f;
+		y[jy] += temp1 * a[j * a_dim1 + 1];
+		l = 1 - j;
+		ix = jx;
+		iy = jy;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    ix += *incx;
+		    iy += *incy;
+		    y[iy] += temp1 * a[l + i__ + j * a_dim1];
+		    temp2 += a[l + i__ + j * a_dim1] * x[ix];
+/* L110: */
+		}
+		y[jy] += *alpha * temp2;
+		jx += *incx;
+		jy += *incy;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of SSBMV . */
+
+} /* ssbmv_ */
+
diff --git a/blas/f2c/sspmv.c b/blas/f2c/sspmv.c
new file mode 100644
index 000000000..47858ec6c
--- /dev/null
+++ b/blas/f2c/sspmv.c
@@ -0,0 +1,316 @@
+/* sspmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int sspmv_(char *uplo, integer *n, real *alpha, real *ap, 
+	real *x, integer *incx, real *beta, real *y, integer *incy, ftnlen 
+	uplo_len)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
+    real temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  SSPMV  performs the matrix-vector operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n symmetric matrix, supplied in packed form. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the matrix A is supplied in the packed */
+/*           array AP as follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  supplied in AP. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  supplied in AP. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  ALPHA  - REAL            . */
+/*           On entry, ALPHA specifies the scalar alpha. */
+/*           Unchanged on exit. */
+
+/*  AP     - REAL             array of DIMENSION at least */
+/*           ( ( n*( n + 1 ) )/2 ). */
+/*           Before entry with UPLO = 'U' or 'u', the array AP must */
+/*           contain the upper triangular part of the symmetric matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
+/*           and a( 2, 2 ) respectively, and so on. */
+/*           Before entry with UPLO = 'L' or 'l', the array AP must */
+/*           contain the lower triangular part of the symmetric matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
+/*           and a( 3, 1 ) respectively, and so on. */
+/*           Unchanged on exit. */
+
+/*  X      - REAL             array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element 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 + ( n - 1 )*abs( INCY ) ). */
+/*           Before entry, the incremented array Y must contain the n */
+/*           element 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --y;
+    --x;
+    --ap;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*incx == 0) {
+	info = 6;
+    } else if (*incy == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("SSPMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (*alpha == 0.f && *beta == 1.f)) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+/*     First form  y := beta*y. */
+
+    if (*beta != 1.f) {
+	if (*incy == 1) {
+	    if (*beta == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = 0.f;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = *beta * y[i__];
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (*beta == 0.f) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = 0.f;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = *beta * y[iy];
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (*alpha == 0.f) {
+	return 0;
+    }
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when AP contains the upper triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.f;
+		k = kk;
+		i__2 = j - 1;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    y[i__] += temp1 * ap[k];
+		    temp2 += ap[k] * x[i__];
+		    ++k;
+/* L50: */
+		}
+		y[j] = y[j] + temp1 * ap[kk + j - 1] + *alpha * temp2;
+		kk += j;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.f;
+		ix = kx;
+		iy = ky;
+		i__2 = kk + j - 2;
+		for (k = kk; k <= i__2; ++k) {
+		    y[iy] += temp1 * ap[k];
+		    temp2 += ap[k] * x[ix];
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		y[jy] = y[jy] + temp1 * ap[kk + j - 1] + *alpha * temp2;
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when AP contains the lower triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[j];
+		temp2 = 0.f;
+		y[j] += temp1 * ap[kk];
+		k = kk + 1;
+		i__2 = *n;
+		for (i__ = j + 1; i__ <= i__2; ++i__) {
+		    y[i__] += temp1 * ap[k];
+		    temp2 += ap[k] * x[i__];
+		    ++k;
+/* L90: */
+		}
+		y[j] += *alpha * temp2;
+		kk += *n - j + 1;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp1 = *alpha * x[jx];
+		temp2 = 0.f;
+		y[jy] += temp1 * ap[kk];
+		ix = jx;
+		iy = jy;
+		i__2 = kk + *n - j;
+		for (k = kk + 1; k <= i__2; ++k) {
+		    ix += *incx;
+		    iy += *incy;
+		    y[iy] += temp1 * ap[k];
+		    temp2 += ap[k] * x[ix];
+/* L110: */
+		}
+		y[jy] += *alpha * temp2;
+		jx += *incx;
+		jy += *incy;
+		kk += *n - j + 1;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of SSPMV . */
+
+} /* sspmv_ */
+
diff --git a/blas/f2c/stbmv.c b/blas/f2c/stbmv.c
new file mode 100644
index 000000000..fcf9ce336
--- /dev/null
+++ b/blas/f2c/stbmv.c
@@ -0,0 +1,428 @@
+/* stbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int stbmv_(char *uplo, char *trans, char *diag, integer *n, 
+	integer *k, real *a, integer *lda, real *x, integer *incx, ftnlen 
+	uplo_len, ftnlen trans_len, ftnlen diag_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    integer i__, j, l, ix, jx, kx, info;
+    real temp;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+    logical nounit;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  STBMV  performs one of the matrix-vector operations */
+
+/*     x := A*x,   or   x := A'*x, */
+
+/*  where x is an n element vector and  A is an n by n unit, or non-unit, */
+/*  upper or lower triangular band matrix, with ( k + 1 ) diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the matrix is an upper or */
+/*           lower triangular matrix as follows: */
+
+/*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
+
+/*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
+
+/*           Unchanged on exit. */
+
+/*  TRANS  - CHARACTER*1. */
+/*           On entry, TRANS specifies the operation to be performed as */
+/*           follows: */
+
+/*              TRANS = 'N' or 'n'   x := A*x. */
+
+/*              TRANS = 'T' or 't'   x := A'*x. */
+
+/*              TRANS = 'C' or 'c'   x := A'*x. */
+
+/*           Unchanged on exit. */
+
+/*  DIAG   - CHARACTER*1. */
+/*           On entry, DIAG specifies whether or not A is unit */
+/*           triangular as follows: */
+
+/*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
+
+/*              DIAG = 'N' or 'n'   A is not assumed to be unit */
+/*                                  triangular. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry with UPLO = 'U' or 'u', K specifies the number of */
+/*           super-diagonals of the matrix A. */
+/*           On entry with UPLO = 'L' or 'l', K specifies the number of */
+/*           sub-diagonals of the matrix A. */
+/*           K must satisfy  0 .le. K. */
+/*           Unchanged on exit. */
+
+/*  A      - REAL             array of DIMENSION ( LDA, n ). */
+/*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer an upper */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer a lower */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Note that when DIAG = 'U' or 'u' the elements of the array A */
+/*           corresponding to the diagonal elements of the matrix are not */
+/*           referenced, but are assumed to be unity. */
+/*           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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - REAL             array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element vector x. On exit, X is overwritten with the */
+/*           tranformed vector x. */
+
+/*  INCX   - INTEGER. */
+/*           On entry, INCX specifies the increment for the elements of */
+/*           X. INCX 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, 
+	    "T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
+	    ftnlen)1)) {
+	info = 2;
+    } else if (! lsame_(diag, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(diag, 
+	    "N", (ftnlen)1, (ftnlen)1)) {
+	info = 3;
+    } else if (*n < 0) {
+	info = 4;
+    } else if (*k < 0) {
+	info = 5;
+    } else if (*lda < *k + 1) {
+	info = 7;
+    } else if (*incx == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("STBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0) {
+	return 0;
+    }
+
+    nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
+
+/*     Set up the start point in X if the increment is not unity. This */
+/*     will be  ( N - 1 )*INCX   too small for descending loops. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of A are */
+/*     accessed sequentially with one pass through A. */
+
+    if (lsame_(trans, "N", (ftnlen)1, (ftnlen)1)) {
+
+/*         Form  x := A*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    if (x[j] != 0.f) {
+			temp = x[j];
+			l = kplus1 - j;
+/* Computing MAX */
+			i__2 = 1, i__3 = j - *k;
+			i__4 = j - 1;
+			for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+			    x[i__] += temp * a[l + i__ + j * a_dim1];
+/* L10: */
+			}
+			if (nounit) {
+			    x[j] *= a[kplus1 + j * a_dim1];
+			}
+		    }
+/* L20: */
+		}
+	    } else {
+		jx = kx;
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    if (x[jx] != 0.f) {
+			temp = x[jx];
+			ix = kx;
+			l = kplus1 - j;
+/* Computing MAX */
+			i__4 = 1, i__2 = j - *k;
+			i__3 = j - 1;
+			for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+			    x[ix] += temp * a[l + i__ + j * a_dim1];
+			    ix += *incx;
+/* L30: */
+			}
+			if (nounit) {
+			    x[jx] *= a[kplus1 + j * a_dim1];
+			}
+		    }
+		    jx += *incx;
+		    if (j > *k) {
+			kx += *incx;
+		    }
+/* L40: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    if (x[j] != 0.f) {
+			temp = x[j];
+			l = 1 - j;
+/* Computing MIN */
+			i__1 = *n, i__3 = j + *k;
+			i__4 = j + 1;
+			for (i__ = min(i__1,i__3); i__ >= i__4; --i__) {
+			    x[i__] += temp * a[l + i__ + j * a_dim1];
+/* L50: */
+			}
+			if (nounit) {
+			    x[j] *= a[j * a_dim1 + 1];
+			}
+		    }
+/* L60: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    if (x[jx] != 0.f) {
+			temp = x[jx];
+			ix = kx;
+			l = 1 - j;
+/* Computing MIN */
+			i__4 = *n, i__1 = j + *k;
+			i__3 = j + 1;
+			for (i__ = min(i__4,i__1); i__ >= i__3; --i__) {
+			    x[ix] += temp * a[l + i__ + j * a_dim1];
+			    ix -= *incx;
+/* L70: */
+			}
+			if (nounit) {
+			    x[jx] *= a[j * a_dim1 + 1];
+			}
+		    }
+		    jx -= *incx;
+		    if (*n - j >= *k) {
+			kx -= *incx;
+		    }
+/* L80: */
+		}
+	    }
+	}
+    } else {
+
+/*        Form  x := A'*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    temp = x[j];
+		    l = kplus1 - j;
+		    if (nounit) {
+			temp *= a[kplus1 + j * a_dim1];
+		    }
+/* Computing MAX */
+		    i__4 = 1, i__1 = j - *k;
+		    i__3 = max(i__4,i__1);
+		    for (i__ = j - 1; i__ >= i__3; --i__) {
+			temp += a[l + i__ + j * a_dim1] * x[i__];
+/* L90: */
+		    }
+		    x[j] = temp;
+/* L100: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    temp = x[jx];
+		    kx -= *incx;
+		    ix = kx;
+		    l = kplus1 - j;
+		    if (nounit) {
+			temp *= a[kplus1 + j * a_dim1];
+		    }
+/* Computing MAX */
+		    i__4 = 1, i__1 = j - *k;
+		    i__3 = max(i__4,i__1);
+		    for (i__ = j - 1; i__ >= i__3; --i__) {
+			temp += a[l + i__ + j * a_dim1] * x[ix];
+			ix -= *incx;
+/* L110: */
+		    }
+		    x[jx] = temp;
+		    jx -= *incx;
+/* L120: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    temp = x[j];
+		    l = 1 - j;
+		    if (nounit) {
+			temp *= a[j * a_dim1 + 1];
+		    }
+/* Computing MIN */
+		    i__1 = *n, i__2 = j + *k;
+		    i__4 = min(i__1,i__2);
+		    for (i__ = j + 1; i__ <= i__4; ++i__) {
+			temp += a[l + i__ + j * a_dim1] * x[i__];
+/* L130: */
+		    }
+		    x[j] = temp;
+/* L140: */
+		}
+	    } else {
+		jx = kx;
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    temp = x[jx];
+		    kx += *incx;
+		    ix = kx;
+		    l = 1 - j;
+		    if (nounit) {
+			temp *= a[j * a_dim1 + 1];
+		    }
+/* Computing MIN */
+		    i__1 = *n, i__2 = j + *k;
+		    i__4 = min(i__1,i__2);
+		    for (i__ = j + 1; i__ <= i__4; ++i__) {
+			temp += a[l + i__ + j * a_dim1] * x[ix];
+			ix += *incx;
+/* L150: */
+		    }
+		    x[jx] = temp;
+		    jx += *incx;
+/* L160: */
+		}
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of STBMV . */
+
+} /* stbmv_ */
+
diff --git a/blas/f2c/zhbmv.c b/blas/f2c/zhbmv.c
new file mode 100644
index 000000000..42da13dbb
--- /dev/null
+++ b/blas/f2c/zhbmv.c
@@ -0,0 +1,488 @@
+/* zhbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int zhbmv_(char *uplo, integer *n, integer *k, doublecomplex 
+	*alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *
+	incx, doublecomplex *beta, doublecomplex *y, integer *incy, ftnlen 
+	uplo_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
+    doublereal d__1;
+    doublecomplex z__1, z__2, z__3, z__4;
+
+    /* Builtin functions */
+    void d_cnjg(doublecomplex *, doublecomplex *);
+
+    /* Local variables */
+    integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
+    doublecomplex temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  ZHBMV  performs the matrix-vector  operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n hermitian band matrix, with k super-diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the band matrix A is being supplied as */
+/*           follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  being supplied. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  being supplied. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry, K specifies the number of super-diagonals of the */
+/*           matrix A. K must satisfy  0 .le. K. */
+/*           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 with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the hermitian matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer the upper */
+/*           triangular part of a hermitian band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the hermitian matrix, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer the lower */
+/*           triangular part of a hermitian band matrix from conventional */
+/*           full matrix storage to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set and are assumed to be zero. */
+/*           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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - COMPLEX*16       array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           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. */
+/*           Unchanged on exit. */
+
+/*  Y      - COMPLEX*16       array of DIMENSION at least */
+/*           ( 1 + ( n - 1 )*abs( INCY ) ). */
+/*           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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+    --y;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*k < 0) {
+	info = 3;
+    } else if (*lda < *k + 1) {
+	info = 6;
+    } else if (*incx == 0) {
+	info = 8;
+    } else if (*incy == 0) {
+	info = 11;
+    }
+    if (info != 0) {
+	xerbla_("ZHBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (alpha->r == 0. && alpha->i == 0. && (beta->r == 1. && 
+                                                         beta->i == 0.))) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array A */
+/*     are accessed sequentially with one pass through A. */
+
+/*     First form  y := beta*y. */
+
+    if (beta->r != 1. || beta->i != 0.) {
+	if (*incy == 1) {
+	    if (beta->r == 0. && beta->i == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    y[i__2].r = 0., y[i__2].i = 0.;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    i__3 = i__;
+		    z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (beta->r == 0. && beta->i == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    y[i__2].r = 0., y[i__2].i = 0.;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    i__3 = iy;
+		    z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (alpha->r == 0. && alpha->i == 0.) {
+	return 0;
+    }
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when upper triangle of A is stored. */
+
+	kplus1 = *k + 1;
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__2 = 1, i__3 = j - *k;
+		i__4 = j - 1;
+		for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+		    i__2 = i__;
+		    i__3 = i__;
+		    i__5 = l + i__ + j * a_dim1;
+		    z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
+		    y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+		    i__2 = i__;
+		    z__2.r = z__3.r * x[i__2].r - z__3.i * x[i__2].i, z__2.i =
+			     z__3.r * x[i__2].i + z__3.i * x[i__2].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+/* L50: */
+		}
+		i__4 = j;
+		i__2 = j;
+		i__3 = kplus1 + j * a_dim1;
+		d__1 = a[i__3].r;
+		z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
+		z__2.r = y[i__2].r + z__3.r, z__2.i = y[i__2].i + z__3.i;
+		z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+		y[i__4].r = z__1.r, y[i__4].i = z__1.i;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__4 = jx;
+		z__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, z__1.i =
+			 alpha->r * x[i__4].i + alpha->i * x[i__4].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		ix = kx;
+		iy = ky;
+		l = kplus1 - j;
+/* Computing MAX */
+		i__4 = 1, i__2 = j - *k;
+		i__3 = j - 1;
+		for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+		    i__4 = iy;
+		    i__2 = iy;
+		    i__5 = l + i__ + j * a_dim1;
+		    z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
+		    y[i__4].r = z__1.r, y[i__4].i = z__1.i;
+		    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+		    i__4 = ix;
+		    z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
+			     z__3.r * x[i__4].i + z__3.i * x[i__4].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		i__3 = jy;
+		i__4 = jy;
+		i__2 = kplus1 + j * a_dim1;
+		d__1 = a[i__2].r;
+		z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
+		z__2.r = y[i__4].r + z__3.r, z__2.i = y[i__4].i + z__3.i;
+		z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+		y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		jx += *incx;
+		jy += *incy;
+		if (j > *k) {
+		    kx += *incx;
+		    ky += *incy;
+		}
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when lower triangle of A is stored. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__3 = j;
+		z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
+			 alpha->r * x[i__3].i + alpha->i * x[i__3].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		i__3 = j;
+		i__4 = j;
+		i__2 = j * a_dim1 + 1;
+		d__1 = a[i__2].r;
+		z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
+		z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		l = 1 - j;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    i__4 = i__;
+		    i__2 = i__;
+		    i__5 = l + i__ + j * a_dim1;
+		    z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
+		    y[i__4].r = z__1.r, y[i__4].i = z__1.i;
+		    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+		    i__4 = i__;
+		    z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
+			     z__3.r * x[i__4].i + z__3.i * x[i__4].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+/* L90: */
+		}
+		i__3 = j;
+		i__4 = j;
+		z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__3 = jx;
+		z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
+			 alpha->r * x[i__3].i + alpha->i * x[i__3].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		i__3 = jy;
+		i__4 = jy;
+		i__2 = j * a_dim1 + 1;
+		d__1 = a[i__2].r;
+		z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
+		z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		l = 1 - j;
+		ix = jx;
+		iy = jy;
+/* Computing MIN */
+		i__4 = *n, i__2 = j + *k;
+		i__3 = min(i__4,i__2);
+		for (i__ = j + 1; i__ <= i__3; ++i__) {
+		    ix += *incx;
+		    iy += *incy;
+		    i__4 = iy;
+		    i__2 = iy;
+		    i__5 = l + i__ + j * a_dim1;
+		    z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
+			    z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
+			    .r;
+		    z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
+		    y[i__4].r = z__1.r, y[i__4].i = z__1.i;
+		    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+		    i__4 = ix;
+		    z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
+			     z__3.r * x[i__4].i + z__3.i * x[i__4].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+/* L110: */
+		}
+		i__3 = jy;
+		i__4 = jy;
+		z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		jx += *incx;
+		jy += *incy;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of ZHBMV . */
+
+} /* zhbmv_ */
+
diff --git a/blas/f2c/zhpmv.c b/blas/f2c/zhpmv.c
new file mode 100644
index 000000000..fbe2f42b3
--- /dev/null
+++ b/blas/f2c/zhpmv.c
@@ -0,0 +1,438 @@
+/* zhpmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int zhpmv_(char *uplo, integer *n, doublecomplex *alpha, 
+	doublecomplex *ap, doublecomplex *x, integer *incx, doublecomplex *
+	beta, doublecomplex *y, integer *incy, ftnlen uplo_len)
+{
+    /* System generated locals */
+    integer i__1, i__2, i__3, i__4, i__5;
+    doublereal d__1;
+    doublecomplex z__1, z__2, z__3, z__4;
+
+    /* Builtin functions */
+    void d_cnjg(doublecomplex *, doublecomplex *);
+
+    /* Local variables */
+    integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
+    doublecomplex temp1, temp2;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  ZHPMV  performs the matrix-vector operation */
+
+/*     y := alpha*A*x + beta*y, */
+
+/*  where alpha and beta are scalars, x and y are n element vectors and */
+/*  A is an n by n hermitian matrix, supplied in packed form. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the upper or lower */
+/*           triangular part of the matrix A is supplied in the packed */
+/*           array AP as follows: */
+
+/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
+/*                                  supplied in AP. */
+
+/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
+/*                                  supplied in AP. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  ALPHA  - COMPLEX*16      . */
+/*           On entry, ALPHA specifies the scalar alpha. */
+/*           Unchanged on exit. */
+
+/*  AP     - COMPLEX*16       array of DIMENSION at least */
+/*           ( ( n*( n + 1 ) )/2 ). */
+/*           Before entry with UPLO = 'U' or 'u', the array AP must */
+/*           contain the upper triangular part of the hermitian matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
+/*           and a( 2, 2 ) respectively, and so on. */
+/*           Before entry with UPLO = 'L' or 'l', the array AP must */
+/*           contain the lower triangular part of the hermitian matrix */
+/*           packed sequentially, column by column, so that AP( 1 ) */
+/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
+/*           and a( 3, 1 ) respectively, and so on. */
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set and are assumed to be zero. */
+/*           Unchanged on exit. */
+
+/*  X      - COMPLEX*16       array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element 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 + ( n - 1 )*abs( INCY ) ). */
+/*           Before entry, the incremented array Y must contain the n */
+/*           element 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --y;
+    --x;
+    --ap;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*incx == 0) {
+	info = 6;
+    } else if (*incy == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("ZHPMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (alpha->r == 0. && alpha->i == 0. && (beta->r == 1. && 
+                                                         beta->i == 0.))) {
+	return 0;
+    }
+
+/*     Set up the start points in  X  and  Y. */
+
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (*n - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (*n - 1) * *incy;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+/*     First form  y := beta*y. */
+
+    if (beta->r != 1. || beta->i != 0.) {
+	if (*incy == 1) {
+	    if (beta->r == 0. && beta->i == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    y[i__2].r = 0., y[i__2].i = 0.;
+/* L10: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = i__;
+		    i__3 = i__;
+		    z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+/* L20: */
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (beta->r == 0. && beta->i == 0.) {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    y[i__2].r = 0., y[i__2].i = 0.;
+		    iy += *incy;
+/* L30: */
+		}
+	    } else {
+		i__1 = *n;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    i__2 = iy;
+		    i__3 = iy;
+		    z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
+			    z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
+			    .r;
+		    y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		    iy += *incy;
+/* L40: */
+		}
+	    }
+	}
+    }
+    if (alpha->r == 0. && alpha->i == 0.) {
+	return 0;
+    }
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  y  when AP contains the upper triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		k = kk;
+		i__2 = j - 1;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    i__3 = i__;
+		    i__4 = i__;
+		    i__5 = k;
+		    z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		    y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		    d_cnjg(&z__3, &ap[k]);
+		    i__3 = i__;
+		    z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i =
+			     z__3.r * x[i__3].i + z__3.i * x[i__3].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    ++k;
+/* L50: */
+		}
+		i__2 = j;
+		i__3 = j;
+		i__4 = kk + j - 1;
+		d__1 = ap[i__4].r;
+		z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
+		z__2.r = y[i__3].r + z__3.r, z__2.i = y[i__3].i + z__3.i;
+		z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+		y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		kk += j;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		ix = kx;
+		iy = ky;
+		i__2 = kk + j - 2;
+		for (k = kk; k <= i__2; ++k) {
+		    i__3 = iy;
+		    i__4 = iy;
+		    i__5 = k;
+		    z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		    y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		    d_cnjg(&z__3, &ap[k]);
+		    i__3 = ix;
+		    z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i =
+			     z__3.r * x[i__3].i + z__3.i * x[i__3].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    ix += *incx;
+		    iy += *incy;
+/* L70: */
+		}
+		i__2 = jy;
+		i__3 = jy;
+		i__4 = kk + j - 1;
+		d__1 = ap[i__4].r;
+		z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
+		z__2.r = y[i__3].r + z__3.r, z__2.i = y[i__3].i + z__3.i;
+		z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+		y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L80: */
+	    }
+	}
+    } else {
+
+/*        Form  y  when AP contains the lower triangle. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		i__2 = j;
+		i__3 = j;
+		i__4 = kk;
+		d__1 = ap[i__4].r;
+		z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
+		z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
+		y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		k = kk + 1;
+		i__2 = *n;
+		for (i__ = j + 1; i__ <= i__2; ++i__) {
+		    i__3 = i__;
+		    i__4 = i__;
+		    i__5 = k;
+		    z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		    y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		    d_cnjg(&z__3, &ap[k]);
+		    i__3 = i__;
+		    z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i =
+			     z__3.r * x[i__3].i + z__3.i * x[i__3].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    ++k;
+/* L90: */
+		}
+		i__2 = j;
+		i__3 = j;
+		z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
+		y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		kk += *n - j + 1;
+/* L100: */
+	    }
+	} else {
+	    jx = kx;
+	    jy = ky;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
+			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
+		temp1.r = z__1.r, temp1.i = z__1.i;
+		temp2.r = 0., temp2.i = 0.;
+		i__2 = jy;
+		i__3 = jy;
+		i__4 = kk;
+		d__1 = ap[i__4].r;
+		z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
+		z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
+		y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		ix = jx;
+		iy = jy;
+		i__2 = kk + *n - j;
+		for (k = kk + 1; k <= i__2; ++k) {
+		    ix += *incx;
+		    iy += *incy;
+		    i__3 = iy;
+		    i__4 = iy;
+		    i__5 = k;
+		    z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 
+			    z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
+			    .r;
+		    z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
+		    y[i__3].r = z__1.r, y[i__3].i = z__1.i;
+		    d_cnjg(&z__3, &ap[k]);
+		    i__3 = ix;
+		    z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i =
+			     z__3.r * x[i__3].i + z__3.i * x[i__3].r;
+		    z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+/* L110: */
+		}
+		i__2 = jy;
+		i__3 = jy;
+		z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 
+			alpha->r * temp2.i + alpha->i * temp2.r;
+		z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
+		y[i__2].r = z__1.r, y[i__2].i = z__1.i;
+		jx += *incx;
+		jy += *incy;
+		kk += *n - j + 1;
+/* L120: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of ZHPMV . */
+
+} /* zhpmv_ */
+
diff --git a/blas/f2c/ztbmv.c b/blas/f2c/ztbmv.c
new file mode 100644
index 000000000..4cdcd7f88
--- /dev/null
+++ b/blas/f2c/ztbmv.c
@@ -0,0 +1,647 @@
+/* ztbmv.f -- translated by f2c (version 20100827).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "datatypes.h"
+
+/* Subroutine */ int ztbmv_(char *uplo, char *trans, char *diag, integer *n, 
+	integer *k, doublecomplex *a, integer *lda, doublecomplex *x, integer 
+	*incx, ftnlen uplo_len, ftnlen trans_len, ftnlen diag_len)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
+    doublecomplex z__1, z__2, z__3;
+
+    /* Builtin functions */
+    void d_cnjg(doublecomplex *, doublecomplex *);
+
+    /* Local variables */
+    integer i__, j, l, ix, jx, kx, info;
+    doublecomplex temp;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+    logical noconj, nounit;
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  ZTBMV  performs one of the matrix-vector operations */
+
+/*     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*x, */
+
+/*  where x is an n element vector and  A is an n by n unit, or non-unit, */
+/*  upper or lower triangular band matrix, with ( k + 1 ) diagonals. */
+
+/*  Arguments */
+/*  ========== */
+
+/*  UPLO   - CHARACTER*1. */
+/*           On entry, UPLO specifies whether the matrix is an upper or */
+/*           lower triangular matrix as follows: */
+
+/*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
+
+/*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
+
+/*           Unchanged on exit. */
+
+/*  TRANS  - CHARACTER*1. */
+/*           On entry, TRANS specifies the operation to be performed as */
+/*           follows: */
+
+/*              TRANS = 'N' or 'n'   x := A*x. */
+
+/*              TRANS = 'T' or 't'   x := A'*x. */
+
+/*              TRANS = 'C' or 'c'   x := conjg( A' )*x. */
+
+/*           Unchanged on exit. */
+
+/*  DIAG   - CHARACTER*1. */
+/*           On entry, DIAG specifies whether or not A is unit */
+/*           triangular as follows: */
+
+/*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
+
+/*              DIAG = 'N' or 'n'   A is not assumed to be unit */
+/*                                  triangular. */
+
+/*           Unchanged on exit. */
+
+/*  N      - INTEGER. */
+/*           On entry, N specifies the order of the matrix A. */
+/*           N must be at least zero. */
+/*           Unchanged on exit. */
+
+/*  K      - INTEGER. */
+/*           On entry with UPLO = 'U' or 'u', K specifies the number of */
+/*           super-diagonals of the matrix A. */
+/*           On entry with UPLO = 'L' or 'l', K specifies the number of */
+/*           sub-diagonals of the matrix A. */
+/*           K must satisfy  0 .le. K. */
+/*           Unchanged on exit. */
+
+/*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ). */
+/*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the upper triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row */
+/*           ( k + 1 ) of the array, the first super-diagonal starting at */
+/*           position 2 in row k, and so on. The top left k by k triangle */
+/*           of the array A is not referenced. */
+/*           The following program segment will transfer an upper */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = K + 1 - J */
+/*                    DO 10, I = MAX( 1, J - K ), J */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
+/*           by n part of the array A must contain the lower triangular */
+/*           band part of the matrix of coefficients, supplied column by */
+/*           column, with the leading diagonal of the matrix in row 1 of */
+/*           the array, the first sub-diagonal starting at position 1 in */
+/*           row 2, and so on. The bottom right k by k triangle of the */
+/*           array A is not referenced. */
+/*           The following program segment will transfer a lower */
+/*           triangular band matrix from conventional full matrix storage */
+/*           to band storage: */
+
+/*                 DO 20, J = 1, N */
+/*                    M = 1 - J */
+/*                    DO 10, I = J, MIN( N, J + K ) */
+/*                       A( M + I, J ) = matrix( I, J ) */
+/*              10    CONTINUE */
+/*              20 CONTINUE */
+
+/*           Note that when DIAG = 'U' or 'u' the elements of the array A */
+/*           corresponding to the diagonal elements of the matrix are not */
+/*           referenced, but are assumed to be unity. */
+/*           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 */
+/*           ( k + 1 ). */
+/*           Unchanged on exit. */
+
+/*  X      - COMPLEX*16       array of dimension at least */
+/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
+/*           Before entry, the incremented array X must contain the n */
+/*           element vector x. On exit, X is overwritten with the */
+/*           tranformed vector x. */
+
+/*  INCX   - INTEGER. */
+/*           On entry, INCX specifies the increment for the elements of */
+/*           X. INCX 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 .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+
+    /* Function Body */
+    info = 0;
+    if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
+	    ftnlen)1, (ftnlen)1)) {
+	info = 1;
+    } else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, 
+	    "T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
+	    ftnlen)1)) {
+	info = 2;
+    } else if (! lsame_(diag, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(diag, 
+	    "N", (ftnlen)1, (ftnlen)1)) {
+	info = 3;
+    } else if (*n < 0) {
+	info = 4;
+    } else if (*k < 0) {
+	info = 5;
+    } else if (*lda < *k + 1) {
+	info = 7;
+    } else if (*incx == 0) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("ZTBMV ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0) {
+	return 0;
+    }
+
+    noconj = lsame_(trans, "T", (ftnlen)1, (ftnlen)1);
+    nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
+
+/*     Set up the start point in X if the increment is not unity. This */
+/*     will be  ( N - 1 )*INCX   too small for descending loops. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of A are */
+/*     accessed sequentially with one pass through A. */
+
+    if (lsame_(trans, "N", (ftnlen)1, (ftnlen)1)) {
+
+/*         Form  x := A*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = j;
+		    if (x[i__2].r != 0. || x[i__2].i != 0.) {
+			i__2 = j;
+			temp.r = x[i__2].r, temp.i = x[i__2].i;
+			l = kplus1 - j;
+/* Computing MAX */
+			i__2 = 1, i__3 = j - *k;
+			i__4 = j - 1;
+			for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
+			    i__2 = i__;
+			    i__3 = i__;
+			    i__5 = l + i__ + j * a_dim1;
+			    z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i, 
+				    z__2.i = temp.r * a[i__5].i + temp.i * a[
+				    i__5].r;
+			    z__1.r = x[i__3].r + z__2.r, z__1.i = x[i__3].i + 
+				    z__2.i;
+			    x[i__2].r = z__1.r, x[i__2].i = z__1.i;
+/* L10: */
+			}
+			if (nounit) {
+			    i__4 = j;
+			    i__2 = j;
+			    i__3 = kplus1 + j * a_dim1;
+			    z__1.r = x[i__2].r * a[i__3].r - x[i__2].i * a[
+				    i__3].i, z__1.i = x[i__2].r * a[i__3].i + 
+				    x[i__2].i * a[i__3].r;
+			    x[i__4].r = z__1.r, x[i__4].i = z__1.i;
+			}
+		    }
+/* L20: */
+		}
+	    } else {
+		jx = kx;
+		i__1 = *n;
+		for (j = 1; j <= i__1; ++j) {
+		    i__4 = jx;
+		    if (x[i__4].r != 0. || x[i__4].i != 0.) {
+			i__4 = jx;
+			temp.r = x[i__4].r, temp.i = x[i__4].i;
+			ix = kx;
+			l = kplus1 - j;
+/* Computing MAX */
+			i__4 = 1, i__2 = j - *k;
+			i__3 = j - 1;
+			for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
+			    i__4 = ix;
+			    i__2 = ix;
+			    i__5 = l + i__ + j * a_dim1;
+			    z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i, 
+				    z__2.i = temp.r * a[i__5].i + temp.i * a[
+				    i__5].r;
+			    z__1.r = x[i__2].r + z__2.r, z__1.i = x[i__2].i + 
+				    z__2.i;
+			    x[i__4].r = z__1.r, x[i__4].i = z__1.i;
+			    ix += *incx;
+/* L30: */
+			}
+			if (nounit) {
+			    i__3 = jx;
+			    i__4 = jx;
+			    i__2 = kplus1 + j * a_dim1;
+			    z__1.r = x[i__4].r * a[i__2].r - x[i__4].i * a[
+				    i__2].i, z__1.i = x[i__4].r * a[i__2].i + 
+				    x[i__4].i * a[i__2].r;
+			    x[i__3].r = z__1.r, x[i__3].i = z__1.i;
+			}
+		    }
+		    jx += *incx;
+		    if (j > *k) {
+			kx += *incx;
+		    }
+/* L40: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    i__1 = j;
+		    if (x[i__1].r != 0. || x[i__1].i != 0.) {
+			i__1 = j;
+			temp.r = x[i__1].r, temp.i = x[i__1].i;
+			l = 1 - j;
+/* Computing MIN */
+			i__1 = *n, i__3 = j + *k;
+			i__4 = j + 1;
+			for (i__ = min(i__1,i__3); i__ >= i__4; --i__) {
+			    i__1 = i__;
+			    i__3 = i__;
+			    i__2 = l + i__ + j * a_dim1;
+			    z__2.r = temp.r * a[i__2].r - temp.i * a[i__2].i, 
+				    z__2.i = temp.r * a[i__2].i + temp.i * a[
+				    i__2].r;
+			    z__1.r = x[i__3].r + z__2.r, z__1.i = x[i__3].i + 
+				    z__2.i;
+			    x[i__1].r = z__1.r, x[i__1].i = z__1.i;
+/* L50: */
+			}
+			if (nounit) {
+			    i__4 = j;
+			    i__1 = j;
+			    i__3 = j * a_dim1 + 1;
+			    z__1.r = x[i__1].r * a[i__3].r - x[i__1].i * a[
+				    i__3].i, z__1.i = x[i__1].r * a[i__3].i + 
+				    x[i__1].i * a[i__3].r;
+			    x[i__4].r = z__1.r, x[i__4].i = z__1.i;
+			}
+		    }
+/* L60: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    i__4 = jx;
+		    if (x[i__4].r != 0. || x[i__4].i != 0.) {
+			i__4 = jx;
+			temp.r = x[i__4].r, temp.i = x[i__4].i;
+			ix = kx;
+			l = 1 - j;
+/* Computing MIN */
+			i__4 = *n, i__1 = j + *k;
+			i__3 = j + 1;
+			for (i__ = min(i__4,i__1); i__ >= i__3; --i__) {
+			    i__4 = ix;
+			    i__1 = ix;
+			    i__2 = l + i__ + j * a_dim1;
+			    z__2.r = temp.r * a[i__2].r - temp.i * a[i__2].i, 
+				    z__2.i = temp.r * a[i__2].i + temp.i * a[
+				    i__2].r;
+			    z__1.r = x[i__1].r + z__2.r, z__1.i = x[i__1].i + 
+				    z__2.i;
+			    x[i__4].r = z__1.r, x[i__4].i = z__1.i;
+			    ix -= *incx;
+/* L70: */
+			}
+			if (nounit) {
+			    i__3 = jx;
+			    i__4 = jx;
+			    i__1 = j * a_dim1 + 1;
+			    z__1.r = x[i__4].r * a[i__1].r - x[i__4].i * a[
+				    i__1].i, z__1.i = x[i__4].r * a[i__1].i + 
+				    x[i__4].i * a[i__1].r;
+			    x[i__3].r = z__1.r, x[i__3].i = z__1.i;
+			}
+		    }
+		    jx -= *incx;
+		    if (*n - j >= *k) {
+			kx -= *incx;
+		    }
+/* L80: */
+		}
+	    }
+	}
+    } else {
+
+/*        Form  x := A'*x  or  x := conjg( A' )*x. */
+
+	if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+	    kplus1 = *k + 1;
+	    if (*incx == 1) {
+		for (j = *n; j >= 1; --j) {
+		    i__3 = j;
+		    temp.r = x[i__3].r, temp.i = x[i__3].i;
+		    l = kplus1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__3 = kplus1 + j * a_dim1;
+			    z__1.r = temp.r * a[i__3].r - temp.i * a[i__3].i, 
+				    z__1.i = temp.r * a[i__3].i + temp.i * a[
+				    i__3].r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    i__4 = l + i__ + j * a_dim1;
+			    i__1 = i__;
+			    z__2.r = a[i__4].r * x[i__1].r - a[i__4].i * x[
+				    i__1].i, z__2.i = a[i__4].r * x[i__1].i + 
+				    a[i__4].i * x[i__1].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+/* L90: */
+			}
+		    } else {
+			if (nounit) {
+			    d_cnjg(&z__2, &a[kplus1 + j * a_dim1]);
+			    z__1.r = temp.r * z__2.r - temp.i * z__2.i, 
+				    z__1.i = temp.r * z__2.i + temp.i * 
+				    z__2.r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+			    i__4 = i__;
+			    z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, 
+				    z__2.i = z__3.r * x[i__4].i + z__3.i * x[
+				    i__4].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+/* L100: */
+			}
+		    }
+		    i__3 = j;
+		    x[i__3].r = temp.r, x[i__3].i = temp.i;
+/* L110: */
+		}
+	    } else {
+		kx += (*n - 1) * *incx;
+		jx = kx;
+		for (j = *n; j >= 1; --j) {
+		    i__3 = jx;
+		    temp.r = x[i__3].r, temp.i = x[i__3].i;
+		    kx -= *incx;
+		    ix = kx;
+		    l = kplus1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__3 = kplus1 + j * a_dim1;
+			    z__1.r = temp.r * a[i__3].r - temp.i * a[i__3].i, 
+				    z__1.i = temp.r * a[i__3].i + temp.i * a[
+				    i__3].r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    i__4 = l + i__ + j * a_dim1;
+			    i__1 = ix;
+			    z__2.r = a[i__4].r * x[i__1].r - a[i__4].i * x[
+				    i__1].i, z__2.i = a[i__4].r * x[i__1].i + 
+				    a[i__4].i * x[i__1].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			    ix -= *incx;
+/* L120: */
+			}
+		    } else {
+			if (nounit) {
+			    d_cnjg(&z__2, &a[kplus1 + j * a_dim1]);
+			    z__1.r = temp.r * z__2.r - temp.i * z__2.i, 
+				    z__1.i = temp.r * z__2.i + temp.i * 
+				    z__2.r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MAX */
+			i__4 = 1, i__1 = j - *k;
+			i__3 = max(i__4,i__1);
+			for (i__ = j - 1; i__ >= i__3; --i__) {
+			    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+			    i__4 = ix;
+			    z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, 
+				    z__2.i = z__3.r * x[i__4].i + z__3.i * x[
+				    i__4].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			    ix -= *incx;
+/* L130: */
+			}
+		    }
+		    i__3 = jx;
+		    x[i__3].r = temp.r, x[i__3].i = temp.i;
+		    jx -= *incx;
+/* L140: */
+		}
+	    }
+	} else {
+	    if (*incx == 1) {
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    i__4 = j;
+		    temp.r = x[i__4].r, temp.i = x[i__4].i;
+		    l = 1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__4 = j * a_dim1 + 1;
+			    z__1.r = temp.r * a[i__4].r - temp.i * a[i__4].i, 
+				    z__1.i = temp.r * a[i__4].i + temp.i * a[
+				    i__4].r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    i__1 = l + i__ + j * a_dim1;
+			    i__2 = i__;
+			    z__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
+				    i__2].i, z__2.i = a[i__1].r * x[i__2].i + 
+				    a[i__1].i * x[i__2].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+/* L150: */
+			}
+		    } else {
+			if (nounit) {
+			    d_cnjg(&z__2, &a[j * a_dim1 + 1]);
+			    z__1.r = temp.r * z__2.r - temp.i * z__2.i, 
+				    z__1.i = temp.r * z__2.i + temp.i * 
+				    z__2.r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+			    i__1 = i__;
+			    z__2.r = z__3.r * x[i__1].r - z__3.i * x[i__1].i, 
+				    z__2.i = z__3.r * x[i__1].i + z__3.i * x[
+				    i__1].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+/* L160: */
+			}
+		    }
+		    i__4 = j;
+		    x[i__4].r = temp.r, x[i__4].i = temp.i;
+/* L170: */
+		}
+	    } else {
+		jx = kx;
+		i__3 = *n;
+		for (j = 1; j <= i__3; ++j) {
+		    i__4 = jx;
+		    temp.r = x[i__4].r, temp.i = x[i__4].i;
+		    kx += *incx;
+		    ix = kx;
+		    l = 1 - j;
+		    if (noconj) {
+			if (nounit) {
+			    i__4 = j * a_dim1 + 1;
+			    z__1.r = temp.r * a[i__4].r - temp.i * a[i__4].i, 
+				    z__1.i = temp.r * a[i__4].i + temp.i * a[
+				    i__4].r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    i__1 = l + i__ + j * a_dim1;
+			    i__2 = ix;
+			    z__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
+				    i__2].i, z__2.i = a[i__1].r * x[i__2].i + 
+				    a[i__1].i * x[i__2].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			    ix += *incx;
+/* L180: */
+			}
+		    } else {
+			if (nounit) {
+			    d_cnjg(&z__2, &a[j * a_dim1 + 1]);
+			    z__1.r = temp.r * z__2.r - temp.i * z__2.i, 
+				    z__1.i = temp.r * z__2.i + temp.i * 
+				    z__2.r;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			}
+/* Computing MIN */
+			i__1 = *n, i__2 = j + *k;
+			i__4 = min(i__1,i__2);
+			for (i__ = j + 1; i__ <= i__4; ++i__) {
+			    d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
+			    i__1 = ix;
+			    z__2.r = z__3.r * x[i__1].r - z__3.i * x[i__1].i, 
+				    z__2.i = z__3.r * x[i__1].i + z__3.i * x[
+				    i__1].r;
+			    z__1.r = temp.r + z__2.r, z__1.i = temp.i + 
+				    z__2.i;
+			    temp.r = z__1.r, temp.i = z__1.i;
+			    ix += *incx;
+/* L190: */
+			}
+		    }
+		    i__4 = jx;
+		    x[i__4].r = temp.r, x[i__4].i = temp.i;
+		    jx += *incx;
+/* L200: */
+		}
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of ZTBMV . */
+
+} /* ztbmv_ */
+
diff --git a/blas/chbmv.f b/blas/fortran/chbmv.f
index 1b1c330ea..1b1c330ea 100644
--- a/blas/chbmv.f
+++ b/blas/fortran/chbmv.f
diff --git a/blas/chpmv.f b/blas/fortran/chpmv.f
index 158be5a7b..158be5a7b 100644
--- a/blas/chpmv.f
+++ b/blas/fortran/chpmv.f
diff --git a/blas/complexdots.f b/blas/fortran/complexdots.f
index a7da51d16..a7da51d16 100644
--- a/blas/complexdots.f
+++ b/blas/fortran/complexdots.f
diff --git a/blas/ctbmv.f b/blas/fortran/ctbmv.f
index 5a879fa01..5a879fa01 100644
--- a/blas/ctbmv.f
+++ b/blas/fortran/ctbmv.f
diff --git a/blas/drotm.f b/blas/fortran/drotm.f
index 63a3b1134..63a3b1134 100644
--- a/blas/drotm.f
+++ b/blas/fortran/drotm.f
diff --git a/blas/drotmg.f b/blas/fortran/drotmg.f
index 3ae647b08..3ae647b08 100644
--- a/blas/drotmg.f
+++ b/blas/fortran/drotmg.f
diff --git a/blas/dsbmv.f b/blas/fortran/dsbmv.f
index 8c82d1fa1..8c82d1fa1 100644
--- a/blas/dsbmv.f
+++ b/blas/fortran/dsbmv.f
diff --git a/blas/dspmv.f b/blas/fortran/dspmv.f
index f6e121e76..f6e121e76 100644
--- a/blas/dspmv.f
+++ b/blas/fortran/dspmv.f
diff --git a/blas/dtbmv.f b/blas/fortran/dtbmv.f
index a87ffdeae..a87ffdeae 100644
--- a/blas/dtbmv.f
+++ b/blas/fortran/dtbmv.f
diff --git a/blas/lsame.f b/blas/fortran/lsame.f
index f53690268..f53690268 100644
--- a/blas/lsame.f
+++ b/blas/fortran/lsame.f
diff --git a/blas/srotm.f b/blas/fortran/srotm.f
index fc5a59333..fc5a59333 100644
--- a/blas/srotm.f
+++ b/blas/fortran/srotm.f
diff --git a/blas/srotmg.f b/blas/fortran/srotmg.f
index 7b3bd4272..7b3bd4272 100644
--- a/blas/srotmg.f
+++ b/blas/fortran/srotmg.f
diff --git a/blas/ssbmv.f b/blas/fortran/ssbmv.f
index 16893a295..16893a295 100644
--- a/blas/ssbmv.f
+++ b/blas/fortran/ssbmv.f
diff --git a/blas/sspmv.f b/blas/fortran/sspmv.f
index 0b8449824..0b8449824 100644
--- a/blas/sspmv.f
+++ b/blas/fortran/sspmv.f
diff --git a/blas/stbmv.f b/blas/fortran/stbmv.f
index c0b8f1136..c0b8f1136 100644
--- a/blas/stbmv.f
+++ b/blas/fortran/stbmv.f
diff --git a/blas/zhbmv.f b/blas/fortran/zhbmv.f
index bca0da5fc..bca0da5fc 100644
--- a/blas/zhbmv.f
+++ b/blas/fortran/zhbmv.f
diff --git a/blas/zhpmv.f b/blas/fortran/zhpmv.f
index b686108b3..b686108b3 100644
--- a/blas/zhpmv.f
+++ b/blas/fortran/zhpmv.f
diff --git a/blas/ztbmv.f b/blas/fortran/ztbmv.f
index 7c85c1b55..7c85c1b55 100644
--- a/blas/ztbmv.f
+++ b/blas/fortran/ztbmv.f