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