Replaced banded/packed BLAS2 stubs with f2c code.

Details:
- Retired the blas2blis wrappers that simply called abort with a "not yet
  implemented" message. This includes all of the level-2 banded and packed
  routines.
- Replaced the aforementioned with the corresponding netlib implementations
  having been run through f2c (with some customization).
- Added directories named 'attic' to build/gen-make-frags/ignore_list.

build/gen-make-frags/ignore_list

@@ -1,3 +1,4 @@
+attic
 broken
 old
 other

frame/compat/bl2_blas.h

@@ -65,22 +65,22 @@
 
 // packed
 
-#include "bla_hpmv.h"
-#include "bla_hpr.h"
-#include "bla_hpr2.h"
-#include "bla_spmv.h"
-#include "bla_spr.h"
-#include "bla_spr2.h"
-#include "bla_tpmv.h"
-#include "bla_tpsv.h"
+//#include "bla_hpmv.h"
+//#include "bla_hpr.h"
+//#include "bla_hpr2.h"
+//#include "bla_spmv.h"
+//#include "bla_spr.h"
+//#include "bla_spr2.h"
+//#include "bla_tpmv.h"
+//#include "bla_tpsv.h"
 
 // banded
 
-#include "bla_gbmv.h"
-#include "bla_hbmv.h"
-#include "bla_sbmv.h"
-#include "bla_tbmv.h"
-#include "bla_tbsv.h"
+//#include "bla_gbmv.h"
+//#include "bla_hbmv.h"
+//#include "bla_sbmv.h"
+//#include "bla_tbmv.h"
+//#include "bla_tbsv.h"
 
 
 // -- Level-3 BLAS --

frame/compat/f2c/bla_hbmv.c

@@ -0,0 +1,986 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* chbmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(c,hbmv)(char *uplo, integer *n, integer *k, complex *
+	alpha, complex *a, integer *lda, complex *x, integer *incx, complex *
+	beta, complex *y, integer *incy)
+{
+    /* 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 info;
+    complex temp1, temp2;
+    integer i__, j, l;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1 * 1;
+    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_ */
+
+/* zhbmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(z,hbmv)(char *uplo, integer *n, integer *k, doublecomplex 
+	*alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *
+	incx, doublecomplex *beta, doublecomplex *y, integer *incy)
+{
+    /* 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 info;
+    doublecomplex temp1, temp2;
+    integer i__, j, l;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1 * 1;
+    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_ */
+
+#endif
+

frame/compat/f2c/bla_hpmv.c

@@ -0,0 +1,890 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* chpmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(c,hpmv)(char *uplo, integer *n, complex *alpha, complex *
+	ap, complex *x, integer *incx, complex *beta, complex *y, integer *
+	incy)
+{
+    /* 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 info;
+    complex temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     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_ */
+
+/* zhpmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(z,hpmv)(char *uplo, integer *n, doublecomplex *alpha, 
+	doublecomplex *ap, doublecomplex *x, integer *incx, doublecomplex *
+	beta, doublecomplex *y, integer *incy)
+{
+    /* 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 info;
+    doublecomplex temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     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_ */
+
+#endif
+

frame/compat/f2c/bla_hpr.c

@@ -0,0 +1,700 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* chpr.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(c,hpr)(char *uplo, integer *n, real *alpha, complex *x, 
+	integer *incx, complex *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2, i__3, i__4, i__5;
+    real r__1;
+    complex q__1, q__2;
+
+    /* Builtin functions */
+    void r_cnjg(complex *, complex *);
+
+    /* Local variables */
+    integer info;
+    complex temp;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, jx, kx = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  CHPR    performs the hermitian rank 1 operation */
+
+/*     A := alpha*x*conjg( x' ) + A, */
+
+/*  where alpha is a real scalar, x is an n element vector and A is an */
+/*  n by n hermitian matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set, they are assumed to be zero, and on exit they */
+/*           are set to zero. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --x;
+
+    /* 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 = 5;
+    }
+    if (info != 0) {
+	xerbla_("CHPR  ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || *alpha == 0.f) {
+	return 0;
+    }
+
+/*     Set the start point in X if the increment is not unity. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	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) {
+		    r_cnjg(&q__2, &x[j]);
+		    q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    k = kk;
+		    i__2 = j - 1;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				q__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
+				q__2.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+			++k;
+/* L10: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = j;
+		    q__1.r = x[i__4].r * temp.r - x[i__4].i * temp.i, q__1.i =
+			     x[i__4].r * temp.i + x[i__4].i * temp.r;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
+		    r_cnjg(&q__2, &x[jx]);
+		    q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    ix = kx;
+		    i__2 = kk + j - 2;
+		    for (k = kk; k <= i__2; ++k) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				q__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
+				q__2.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+			ix += *incx;
+/* L30: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = jx;
+		    q__1.r = x[i__4].r * temp.r - x[i__4].i * temp.i, q__1.i =
+			     x[i__4].r * temp.i + x[i__4].i * temp.r;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		jx += *incx;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	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) {
+		    r_cnjg(&q__2, &x[j]);
+		    q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = j;
+		    q__1.r = temp.r * x[i__4].r - temp.i * x[i__4].i, q__1.i =
+			     temp.r * x[i__4].i + temp.i * x[i__4].r;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		    k = kk + 1;
+		    i__2 = *n;
+		    for (i__ = j + 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				q__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
+				q__2.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+			++k;
+/* L50: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
+		    r_cnjg(&q__2, &x[jx]);
+		    q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = jx;
+		    q__1.r = temp.r * x[i__4].r - temp.i * x[i__4].i, q__1.i =
+			     temp.r * x[i__4].i + temp.i * x[i__4].r;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		    ix = jx;
+		    i__2 = kk + *n - j;
+		    for (k = kk + 1; k <= i__2; ++k) {
+			ix += *incx;
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				q__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
+				q__2.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+/* L70: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		jx += *incx;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of CHPR  . */
+
+} /* chpr_ */
+
+/* zhpr.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(z,hpr)(char *uplo, integer *n, doublereal *alpha, 
+	doublecomplex *x, integer *incx, doublecomplex *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2, i__3, i__4, i__5;
+    doublereal d__1;
+    doublecomplex z__1, z__2;
+
+    /* Builtin functions */
+    void d_cnjg(doublecomplex *, doublecomplex *);
+
+    /* Local variables */
+    integer info;
+    doublecomplex temp;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, jx, kx = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  ZHPR    performs the hermitian rank 1 operation */
+
+/*     A := alpha*x*conjg( x' ) + A, */
+
+/*  where alpha is a real scalar, x is an n element vector and A is an */
+/*  n by n hermitian matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set, they are assumed to be zero, and on exit they */
+/*           are set to zero. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --x;
+
+    /* 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 = 5;
+    }
+    if (info != 0) {
+	xerbla_("ZHPR  ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || *alpha == 0.) {
+	return 0;
+    }
+
+/*     Set the start point in X if the increment is not unity. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	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.) {
+		    d_cnjg(&z__2, &x[j]);
+		    z__1.r = *alpha * z__2.r, z__1.i = *alpha * z__2.i;
+		    temp.r = z__1.r, temp.i = z__1.i;
+		    k = kk;
+		    i__2 = j - 1;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				z__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
+				z__2.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+			++k;
+/* L10: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = j;
+		    z__1.r = x[i__4].r * temp.r - x[i__4].i * temp.i, z__1.i =
+			     x[i__4].r * temp.i + x[i__4].i * temp.r;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		if (x[i__2].r != 0. || x[i__2].i != 0.) {
+		    d_cnjg(&z__2, &x[jx]);
+		    z__1.r = *alpha * z__2.r, z__1.i = *alpha * z__2.i;
+		    temp.r = z__1.r, temp.i = z__1.i;
+		    ix = kx;
+		    i__2 = kk + j - 2;
+		    for (k = kk; k <= i__2; ++k) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				z__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
+				z__2.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+			ix += *incx;
+/* L30: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = jx;
+		    z__1.r = x[i__4].r * temp.r - x[i__4].i * temp.i, z__1.i =
+			     x[i__4].r * temp.i + x[i__4].i * temp.r;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		jx += *incx;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	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.) {
+		    d_cnjg(&z__2, &x[j]);
+		    z__1.r = *alpha * z__2.r, z__1.i = *alpha * z__2.i;
+		    temp.r = z__1.r, temp.i = z__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = j;
+		    z__1.r = temp.r * x[i__4].r - temp.i * x[i__4].i, z__1.i =
+			     temp.r * x[i__4].i + temp.i * x[i__4].r;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		    k = kk + 1;
+		    i__2 = *n;
+		    for (i__ = j + 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				z__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
+				z__2.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+			++k;
+/* L50: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		if (x[i__2].r != 0. || x[i__2].i != 0.) {
+		    d_cnjg(&z__2, &x[jx]);
+		    z__1.r = *alpha * z__2.r, z__1.i = *alpha * z__2.i;
+		    temp.r = z__1.r, temp.i = z__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = jx;
+		    z__1.r = temp.r * x[i__4].r - temp.i * x[i__4].i, z__1.i =
+			     temp.r * x[i__4].i + temp.i * x[i__4].r;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		    ix = jx;
+		    i__2 = kk + *n - j;
+		    for (k = kk + 1; k <= i__2; ++k) {
+			ix += *incx;
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
+				z__2.i = x[i__5].r * temp.i + x[i__5].i * 
+				temp.r;
+			z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
+				z__2.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+/* L70: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		jx += *incx;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of ZHPR  . */
+
+} /* zhpr_ */
+
+#endif
+

frame/compat/f2c/bla_hpr2.c

@@ -0,0 +1,917 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* chpr2.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(c,hpr2)(char *uplo, integer *n, complex *alpha, complex *
+	x, integer *incx, complex *y, integer *incy, complex *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2, i__3, i__4, i__5, i__6;
+    real r__1;
+    complex q__1, q__2, q__3, q__4;
+
+    /* Builtin functions */
+    void r_cnjg(complex *, complex *);
+
+    /* Local variables */
+    integer info;
+    complex temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx = 0, jy = 0, kx = 0, ky = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  CHPR2  performs the hermitian rank 2 operation */
+
+/*     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, */
+
+/*  where alpha is a scalar, x and y are n element vectors and A is an */
+/*  n by n hermitian matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. */
+/*           Unchanged on exit. */
+
+/*  INCY   - INTEGER. */
+/*           On entry, INCY specifies the increment for the elements of */
+/*           Y. INCY must not be zero. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set, they are assumed to be zero, and on exit they */
+/*           are set to zero. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --y;
+    --x;
+
+    /* 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 = 5;
+    } else if (*incy == 0) {
+	info = 7;
+    }
+    if (info != 0) {
+	xerbla_("CHPR2 ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f)) {
+	return 0;
+    }
+
+/*     Set up the start points in X and Y if the increments are not both */
+/*     unity. */
+
+    if (*incx != 1 || *incy != 1) {
+	if (*incx > 0) {
+	    kx = 1;
+	} else {
+	    kx = 1 - (*n - 1) * *incx;
+	}
+	if (*incy > 0) {
+	    ky = 1;
+	} else {
+	    ky = 1 - (*n - 1) * *incy;
+	}
+	jx = kx;
+	jy = ky;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		i__3 = j;
+		if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
+			|| y[i__3].i != 0.f)) {
+		    r_cnjg(&q__2, &y[j]);
+		    q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
+			    alpha->r * q__2.i + alpha->i * q__2.r;
+		    temp1.r = q__1.r, temp1.i = q__1.i;
+		    i__2 = j;
+		    q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    r_cnjg(&q__1, &q__2);
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    k = kk;
+		    i__2 = j - 1;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
+				q__3.i;
+			i__6 = i__;
+			q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+			++k;
+/* L10: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = j;
+		    q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = j;
+		    q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		i__3 = jy;
+		if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
+			|| y[i__3].i != 0.f)) {
+		    r_cnjg(&q__2, &y[jy]);
+		    q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
+			    alpha->r * q__2.i + alpha->i * q__2.r;
+		    temp1.r = q__1.r, temp1.i = q__1.i;
+		    i__2 = jx;
+		    q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    r_cnjg(&q__1, &q__2);
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    ix = kx;
+		    iy = ky;
+		    i__2 = kk + j - 2;
+		    for (k = kk; k <= i__2; ++k) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
+				q__3.i;
+			i__6 = iy;
+			q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+			ix += *incx;
+			iy += *incy;
+/* L30: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = jx;
+		    q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = jy;
+		    q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		i__3 = j;
+		if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
+			|| y[i__3].i != 0.f)) {
+		    r_cnjg(&q__2, &y[j]);
+		    q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
+			    alpha->r * q__2.i + alpha->i * q__2.r;
+		    temp1.r = q__1.r, temp1.i = q__1.i;
+		    i__2 = j;
+		    q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    r_cnjg(&q__1, &q__2);
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = j;
+		    q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = j;
+		    q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		    k = kk + 1;
+		    i__2 = *n;
+		    for (i__ = j + 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
+				q__3.i;
+			i__6 = i__;
+			q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+			++k;
+/* L50: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		i__3 = jy;
+		if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
+			|| y[i__3].i != 0.f)) {
+		    r_cnjg(&q__2, &y[jy]);
+		    q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
+			    alpha->r * q__2.i + alpha->i * q__2.r;
+		    temp1.r = q__1.r, temp1.i = q__1.i;
+		    i__2 = jx;
+		    q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    r_cnjg(&q__1, &q__2);
+		    temp2.r = q__1.r, temp2.i = q__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = jx;
+		    q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = jy;
+		    q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+		    r__1 = ap[i__3].r + q__1.r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		    ix = jx;
+		    iy = jy;
+		    i__2 = kk + *n - j;
+		    for (k = kk + 1; k <= i__2; ++k) {
+			ix += *incx;
+			iy += *incy;
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
+				q__3.i;
+			i__6 = iy;
+			q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
+			ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
+/* L70: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    r__1 = ap[i__3].r;
+		    ap[i__2].r = r__1, ap[i__2].i = 0.f;
+		}
+		jx += *incx;
+		jy += *incy;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of CHPR2 . */
+
+} /* chpr2_ */
+
+/* zhpr2.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(z,hpr2)(char *uplo, integer *n, doublecomplex *alpha, 
+	doublecomplex *x, integer *incx, doublecomplex *y, integer *incy, 
+	doublecomplex *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2, i__3, i__4, i__5, i__6;
+    doublereal d__1;
+    doublecomplex z__1, z__2, z__3, z__4;
+
+    /* Builtin functions */
+    void d_cnjg(doublecomplex *, doublecomplex *);
+
+    /* Local variables */
+    integer info;
+    doublecomplex temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx = 0, jy = 0, kx = 0, ky = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  ZHPR2  performs the hermitian rank 2 operation */
+
+/*     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, */
+
+/*  where alpha is a scalar, x and y are n element vectors and A is an */
+/*  n by n hermitian matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. */
+/*           Unchanged on exit. */
+
+/*  INCY   - INTEGER. */
+/*           On entry, INCY specifies the increment for the elements of */
+/*           Y. INCY must not be zero. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+/*           Note that the imaginary parts of the diagonal elements need */
+/*           not be set, they are assumed to be zero, and on exit they */
+/*           are set to zero. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --y;
+    --x;
+
+    /* 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 = 5;
+    } else if (*incy == 0) {
+	info = 7;
+    }
+    if (info != 0) {
+	xerbla_("ZHPR2 ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || (alpha->r == 0. && alpha->i == 0.)) {
+	return 0;
+    }
+
+/*     Set up the start points in X and Y if the increments are not both */
+/*     unity. */
+
+    if (*incx != 1 || *incy != 1) {
+	if (*incx > 0) {
+	    kx = 1;
+	} else {
+	    kx = 1 - (*n - 1) * *incx;
+	}
+	if (*incy > 0) {
+	    ky = 1;
+	} else {
+	    ky = 1 - (*n - 1) * *incy;
+	}
+	jx = kx;
+	jy = ky;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		i__3 = j;
+		if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. || 
+			y[i__3].i != 0.)) {
+		    d_cnjg(&z__2, &y[j]);
+		    z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i = 
+			    alpha->r * z__2.i + alpha->i * z__2.r;
+		    temp1.r = z__1.r, temp1.i = z__1.i;
+		    i__2 = j;
+		    z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    d_cnjg(&z__1, &z__2);
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    k = kk;
+		    i__2 = j - 1;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				z__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			z__2.r = ap[i__4].r + z__3.r, z__2.i = ap[i__4].i + 
+				z__3.i;
+			i__6 = i__;
+			z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				z__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+			++k;
+/* L10: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = j;
+		    z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    z__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = j;
+		    z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    z__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		i__3 = jy;
+		if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. || 
+			y[i__3].i != 0.)) {
+		    d_cnjg(&z__2, &y[jy]);
+		    z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i = 
+			    alpha->r * z__2.i + alpha->i * z__2.r;
+		    temp1.r = z__1.r, temp1.i = z__1.i;
+		    i__2 = jx;
+		    z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    d_cnjg(&z__1, &z__2);
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    ix = kx;
+		    iy = ky;
+		    i__2 = kk + j - 2;
+		    for (k = kk; k <= i__2; ++k) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				z__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			z__2.r = ap[i__4].r + z__3.r, z__2.i = ap[i__4].i + 
+				z__3.i;
+			i__6 = iy;
+			z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				z__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+			ix += *incx;
+			iy += *incy;
+/* L30: */
+		    }
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    i__4 = jx;
+		    z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    z__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = jy;
+		    z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    z__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		} else {
+		    i__2 = kk + j - 1;
+		    i__3 = kk + j - 1;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = j;
+		i__3 = j;
+		if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. || 
+			y[i__3].i != 0.)) {
+		    d_cnjg(&z__2, &y[j]);
+		    z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i = 
+			    alpha->r * z__2.i + alpha->i * z__2.r;
+		    temp1.r = z__1.r, temp1.i = z__1.i;
+		    i__2 = j;
+		    z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    d_cnjg(&z__1, &z__2);
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = j;
+		    z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    z__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = j;
+		    z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    z__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		    k = kk + 1;
+		    i__2 = *n;
+		    for (i__ = j + 1; i__ <= i__2; ++i__) {
+			i__3 = k;
+			i__4 = k;
+			i__5 = i__;
+			z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				z__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			z__2.r = ap[i__4].r + z__3.r, z__2.i = ap[i__4].i + 
+				z__3.i;
+			i__6 = i__;
+			z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				z__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+			++k;
+/* L50: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = jx;
+		i__3 = jy;
+		if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. || 
+			y[i__3].i != 0.)) {
+		    d_cnjg(&z__2, &y[jy]);
+		    z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i = 
+			    alpha->r * z__2.i + alpha->i * z__2.r;
+		    temp1.r = z__1.r, temp1.i = z__1.i;
+		    i__2 = jx;
+		    z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
+			    z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
+			    .r;
+		    d_cnjg(&z__1, &z__2);
+		    temp2.r = z__1.r, temp2.i = z__1.i;
+		    i__2 = kk;
+		    i__3 = kk;
+		    i__4 = jx;
+		    z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
+			    z__2.i = x[i__4].r * temp1.i + x[i__4].i * 
+			    temp1.r;
+		    i__5 = jy;
+		    z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
+			    z__3.i = y[i__5].r * temp2.i + y[i__5].i * 
+			    temp2.r;
+		    z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
+		    d__1 = ap[i__3].r + z__1.r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		    ix = jx;
+		    iy = jy;
+		    i__2 = kk + *n - j;
+		    for (k = kk + 1; k <= i__2; ++k) {
+			ix += *incx;
+			iy += *incy;
+			i__3 = k;
+			i__4 = k;
+			i__5 = ix;
+			z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
+				z__3.i = x[i__5].r * temp1.i + x[i__5].i * 
+				temp1.r;
+			z__2.r = ap[i__4].r + z__3.r, z__2.i = ap[i__4].i + 
+				z__3.i;
+			i__6 = iy;
+			z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
+				z__4.i = y[i__6].r * temp2.i + y[i__6].i * 
+				temp2.r;
+			z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
+			ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
+/* L70: */
+		    }
+		} else {
+		    i__2 = kk;
+		    i__3 = kk;
+		    d__1 = ap[i__3].r;
+		    ap[i__2].r = d__1, ap[i__2].i = 0.;
+		}
+		jx += *incx;
+		jy += *incy;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of ZHPR2 . */
+
+} /* zhpr2_ */
+
+#endif
+

frame/compat/f2c/bla_sbmv.c

@@ -0,0 +1,748 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* dsbmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(d,sbmv)(char *uplo, integer *n, integer *k, doublereal *
+	alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, 
+	doublereal *beta, doublereal *y, integer *incy)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    integer info;
+    doublereal temp1, temp2;
+    integer i__, j, l;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1 * 1;
+    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_ */
+
+/* ssbmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(s,sbmv)(char *uplo, integer *n, integer *k, real *alpha, 
+	real *a, integer *lda, real *x, integer *incx, real *beta, real *y, 
+	integer *incy)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    integer info;
+    real temp1, temp2;
+    integer i__, j, l;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kplus1, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1 * 1;
+    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_ */
+
+#endif
+

frame/compat/f2c/bla_spmv.c

@@ -0,0 +1,647 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* dspmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(d,spmv)(char *uplo, integer *n, doublereal *alpha, 
+	doublereal *ap, doublereal *x, integer *incx, doublereal *beta, 
+	doublereal *y, integer *incy)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer info;
+    doublereal temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     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_ */
+
+/* sspmv.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(s,spmv)(char *uplo, integer *n, real *alpha, real *ap, 
+	real *x, integer *incx, real *beta, real *y, integer *incy)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer info;
+    real temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx, jy, kx, ky;
+    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. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     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_ */
+
+#endif
+

frame/compat/f2c/bla_spr.c

@@ -0,0 +1,498 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* dspr.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(d,spr)(char *uplo, integer *n, doublereal *alpha, 
+	doublereal *x, integer *incx, doublereal *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer info;
+    doublereal temp;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, jx, kx = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  DSPR    performs the symmetric rank 1 operation */
+
+/*     A := alpha*x*x' + A, */
+
+/*  where alpha is a real scalar, x is an n element vector and A is an */
+/*  n by n symmetric matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --x;
+
+    /* 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 = 5;
+    }
+    if (info != 0) {
+	xerbla_("DSPR  ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || *alpha == 0.) {
+	return 0;
+    }
+
+/*     Set the start point in X if the increment is not unity. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	if (*incx == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0.) {
+		    temp = *alpha * x[j];
+		    k = kk;
+		    i__2 = j;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			ap[k] += x[i__] * temp;
+			++k;
+/* L10: */
+		    }
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0.) {
+		    temp = *alpha * x[jx];
+		    ix = kx;
+		    i__2 = kk + j - 1;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] += x[ix] * temp;
+			ix += *incx;
+/* L30: */
+		    }
+		}
+		jx += *incx;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	if (*incx == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0.) {
+		    temp = *alpha * x[j];
+		    k = kk;
+		    i__2 = *n;
+		    for (i__ = j; i__ <= i__2; ++i__) {
+			ap[k] += x[i__] * temp;
+			++k;
+/* L50: */
+		    }
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0.) {
+		    temp = *alpha * x[jx];
+		    ix = jx;
+		    i__2 = kk + *n - j;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] += x[ix] * temp;
+			ix += *incx;
+/* L70: */
+		    }
+		}
+		jx += *incx;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of DSPR  . */
+
+} /* dspr_ */
+
+/* sspr.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(s,spr)(char *uplo, integer *n, real *alpha, real *x, 
+	integer *incx, real *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer info;
+    real temp;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, jx, kx = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  SSPR    performs the symmetric rank 1 operation */
+
+/*     A := alpha*x*x' + A, */
+
+/*  where alpha is a real scalar, x is an n element vector and A is an */
+/*  n by n symmetric matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --x;
+
+    /* 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 = 5;
+    }
+    if (info != 0) {
+	xerbla_("SSPR  ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || *alpha == 0.f) {
+	return 0;
+    }
+
+/*     Set the start point in X if the increment is not unity. */
+
+    if (*incx <= 0) {
+	kx = 1 - (*n - 1) * *incx;
+    } else if (*incx != 1) {
+	kx = 1;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	if (*incx == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0.f) {
+		    temp = *alpha * x[j];
+		    k = kk;
+		    i__2 = j;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			ap[k] += x[i__] * temp;
+			++k;
+/* L10: */
+		    }
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0.f) {
+		    temp = *alpha * x[jx];
+		    ix = kx;
+		    i__2 = kk + j - 1;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] += x[ix] * temp;
+			ix += *incx;
+/* L30: */
+		    }
+		}
+		jx += *incx;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	if (*incx == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0.f) {
+		    temp = *alpha * x[j];
+		    k = kk;
+		    i__2 = *n;
+		    for (i__ = j; i__ <= i__2; ++i__) {
+			ap[k] += x[i__] * temp;
+			++k;
+/* L50: */
+		    }
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    jx = kx;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0.f) {
+		    temp = *alpha * x[jx];
+		    ix = jx;
+		    i__2 = kk + *n - j;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] += x[ix] * temp;
+			ix += *incx;
+/* L70: */
+		    }
+		}
+		jx += *incx;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of SSPR  . */
+
+} /* sspr_ */
+
+#endif
+

frame/compat/f2c/bla_spr2.c

@@ -0,0 +1,563 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+/* dspr2.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(d,spr2)(char *uplo, integer *n, doublereal *alpha, 
+	doublereal *x, integer *incx, doublereal *y, integer *incy, 
+	doublereal *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer info;
+    doublereal temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx = 0, jy = 0, kx = 0, ky = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  DSPR2  performs the symmetric rank 2 operation */
+
+/*     A := alpha*x*y' + alpha*y*x' + A, */
+
+/*  where alpha is a scalar, x and y are n element vectors and A is an */
+/*  n by n symmetric matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. */
+/*           Unchanged on exit. */
+
+/*  INCY   - INTEGER. */
+/*           On entry, INCY specifies the increment for the elements of */
+/*           Y. INCY must not be zero. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --y;
+    --x;
+
+    /* 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 = 5;
+    } else if (*incy == 0) {
+	info = 7;
+    }
+    if (info != 0) {
+	xerbla_("DSPR2 ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || *alpha == 0.) {
+	return 0;
+    }
+
+/*     Set up the start points in X and Y if the increments are not both */
+/*     unity. */
+
+    if (*incx != 1 || *incy != 1) {
+	if (*incx > 0) {
+	    kx = 1;
+	} else {
+	    kx = 1 - (*n - 1) * *incx;
+	}
+	if (*incy > 0) {
+	    ky = 1;
+	} else {
+	    ky = 1 - (*n - 1) * *incy;
+	}
+	jx = kx;
+	jy = ky;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0. || y[j] != 0.) {
+		    temp1 = *alpha * y[j];
+		    temp2 = *alpha * x[j];
+		    k = kk;
+		    i__2 = j;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
+			++k;
+/* L10: */
+		    }
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0. || y[jy] != 0.) {
+		    temp1 = *alpha * y[jy];
+		    temp2 = *alpha * x[jx];
+		    ix = kx;
+		    iy = ky;
+		    i__2 = kk + j - 1;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
+			ix += *incx;
+			iy += *incy;
+/* L30: */
+		    }
+		}
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0. || y[j] != 0.) {
+		    temp1 = *alpha * y[j];
+		    temp2 = *alpha * x[j];
+		    k = kk;
+		    i__2 = *n;
+		    for (i__ = j; i__ <= i__2; ++i__) {
+			ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
+			++k;
+/* L50: */
+		    }
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0. || y[jy] != 0.) {
+		    temp1 = *alpha * y[jy];
+		    temp2 = *alpha * x[jx];
+		    ix = jx;
+		    iy = jy;
+		    i__2 = kk + *n - j;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
+			ix += *incx;
+			iy += *incy;
+/* L70: */
+		    }
+		}
+		jx += *incx;
+		jy += *incy;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of DSPR2 . */
+
+} /* dspr2_ */
+
+/* sspr2.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Subroutine */ int PASTEF77(s,spr2)(char *uplo, integer *n, real *alpha, real *x, 
+	integer *incx, real *y, integer *incy, real *ap)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+
+    /* Local variables */
+    integer info;
+    real temp1, temp2;
+    integer i__, j, k;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    integer kk, ix, iy, jx = 0, jy = 0, kx = 0, ky = 0;
+    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
+
+/*     .. Scalar Arguments .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  SSPR2  performs the symmetric rank 2 operation */
+
+/*     A := alpha*x*y' + alpha*y*x' + A, */
+
+/*  where alpha is a scalar, x and y are n element vectors and A is an */
+/*  n by n symmetric matrix, supplied in packed form. */
+
+/*  Parameters */
+/*  ========== */
+
+/*  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. */
+
+/*  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. */
+
+/*  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. */
+/*           Unchanged on exit. */
+
+/*  INCY   - INTEGER. */
+/*           On entry, INCY specifies the increment for the elements of */
+/*           Y. INCY must not be zero. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the upper triangular part of the */
+/*           updated matrix. */
+/*           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. On exit, the array */
+/*           AP is overwritten by the lower triangular part of the */
+/*           updated matrix. */
+
+
+/*  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 .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Test the input parameters. */
+
+    /* Parameter adjustments */
+    --ap;
+    --y;
+    --x;
+
+    /* 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 = 5;
+    } else if (*incy == 0) {
+	info = 7;
+    }
+    if (info != 0) {
+	xerbla_("SSPR2 ", &info, (ftnlen)6);
+	return 0;
+    }
+
+/*     Quick return if possible. */
+
+    if (*n == 0 || *alpha == 0.f) {
+	return 0;
+    }
+
+/*     Set up the start points in X and Y if the increments are not both */
+/*     unity. */
+
+    if (*incx != 1 || *incy != 1) {
+	if (*incx > 0) {
+	    kx = 1;
+	} else {
+	    kx = 1 - (*n - 1) * *incx;
+	}
+	if (*incy > 0) {
+	    ky = 1;
+	} else {
+	    ky = 1 - (*n - 1) * *incy;
+	}
+	jx = kx;
+	jy = ky;
+    }
+
+/*     Start the operations. In this version the elements of the array AP */
+/*     are accessed sequentially with one pass through AP. */
+
+    kk = 1;
+    if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
+
+/*        Form  A  when upper triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0.f || y[j] != 0.f) {
+		    temp1 = *alpha * y[j];
+		    temp2 = *alpha * x[j];
+		    k = kk;
+		    i__2 = j;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
+			++k;
+/* L10: */
+		    }
+		}
+		kk += j;
+/* L20: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0.f || y[jy] != 0.f) {
+		    temp1 = *alpha * y[jy];
+		    temp2 = *alpha * x[jx];
+		    ix = kx;
+		    iy = ky;
+		    i__2 = kk + j - 1;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
+			ix += *incx;
+			iy += *incy;
+/* L30: */
+		    }
+		}
+		jx += *incx;
+		jy += *incy;
+		kk += j;
+/* L40: */
+	    }
+	}
+    } else {
+
+/*        Form  A  when lower triangle is stored in AP. */
+
+	if (*incx == 1 && *incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[j] != 0.f || y[j] != 0.f) {
+		    temp1 = *alpha * y[j];
+		    temp2 = *alpha * x[j];
+		    k = kk;
+		    i__2 = *n;
+		    for (i__ = j; i__ <= i__2; ++i__) {
+			ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
+			++k;
+/* L50: */
+		    }
+		}
+		kk = kk + *n - j + 1;
+/* L60: */
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (x[jx] != 0.f || y[jy] != 0.f) {
+		    temp1 = *alpha * y[jy];
+		    temp2 = *alpha * x[jx];
+		    ix = jx;
+		    iy = jy;
+		    i__2 = kk + *n - j;
+		    for (k = kk; k <= i__2; ++k) {
+			ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
+			ix += *incx;
+			iy += *incy;
+/* L70: */
+		    }
+		}
+		jx += *incx;
+		jy += *incy;
+		kk = kk + *n - j + 1;
+/* L80: */
+	    }
+	}
+    }
+
+    return 0;
+
+/*     End of SSPR2 . */
+
+} /* sspr2_ */
+
+#endif
+

frame/compat/f2c/bla_xerbla.c

@@ -0,0 +1,88 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+#include "stdio.h"
+
+/* xerbla.f -- translated by f2c (version 19991025).
+   You must link the resulting object file with the libraries:
+	-lf2c -lm   (in that order)
+*/
+
+/* Table of constant values */
+
+/* Subroutine */ int xerbla_(char *srname, integer *info, ftnlen srname_len)
+{
+/*  -- LAPACK auxiliary routine (preliminary version) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
+/*     Courant Institute, Argonne National Lab, and Rice University */
+/*     February 29, 1992 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  XERBLA  is an error handler for the LAPACK routines. */
+/*  It is called by an LAPACK routine if an input parameter has an */
+/*  invalid value.  A message is printed and execution stops. */
+
+/*  Installers may consider modifying the STOP statement in order to */
+/*  call system-specific exception-handling facilities. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  SRNAME  (input) CHARACTER*6 */
+/*          The name of the routine which called XERBLA. */
+
+/*  INFO    (input) INTEGER */
+/*          The position of the invalid parameter in the parameter list */
+/*          of the calling routine. */
+
+    printf("** On entry to %6s, parameter number %2i had an illegal value\n",
+        srname, (int)*info);
+
+/*     End of XERBLA */
+
+    return 0;
+} /* xerbla_ */
+
+#endif
+

frame/compat/f2c/util/bla_c_div.c

@@ -0,0 +1,54 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+void c_div(complex *cp, complex *ap, complex *bp)
+{
+	complex a = *ap;
+	complex b = *bp;
+	real    temp;
+
+	temp = b.r * b.r + b.i * b.i;
+
+	cp->r = ( a.r * b.r + a.i * b.i ) / temp;
+	cp->i = ( a.i * b.r - a.r * b.i ) / temp;
+}
+
+#endif
+

frame/compat/f2c/util/bla_d_cnjg.c

@@ -0,0 +1,48 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+void d_cnjg(doublecomplex *dest, doublecomplex *src)
+{
+	dest->r =   src->r ;
+	dest->i = -(src->i);
+}
+
+#endif
+

frame/compat/f2c/util/bla_d_imag.c

@@ -0,0 +1,47 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+double d_imag(doublecomplex *z)
+{
+    return(z->i);
+}
+
+#endif
+

frame/compat/f2c/util/bla_r_cnjg.c

@@ -0,0 +1,48 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+void r_cnjg(complex *dest, complex *src)
+{
+	dest->r =   src->r ;
+	dest->i = -(src->i);
+}
+
+#endif
+

frame/compat/f2c/util/bla_r_imag.c

@@ -0,0 +1,47 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+real r_imag(complex *z)
+{
+    return(z->i);
+}
+
+#endif
+

frame/compat/f2c/util/bla_z_div.c

@@ -0,0 +1,54 @@
+/*
+
+   BLIS    
+   An object-based framework for developing high-performance BLAS-like
+   libraries.
+
+   Copyright (C) 2013, The University of Texas
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions are
+   met:
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    - Neither the name of The University of Texas nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+*/
+
+#include "blis2.h"
+
+#ifdef BLIS_ENABLE_BLAS2BLIS
+
+#include "bl2_f2c.h"
+
+void z_div(doublecomplex *cp, doublecomplex *ap, doublecomplex *bp)
+{
+	doublecomplex a = *ap;
+	doublecomplex b = *bp;
+	double        temp;
+
+	temp = b.r * b.r + b.i * b.i;
+
+	cp->r = ( a.r * b.r + a.i * b.i ) / temp;
+	cp->i = ( a.i * b.r - a.r * b.i ) / temp;
+}
+
+#endif
+

version

@@ -1 +1 @@
-0.0.3
+0.0.3-1