| /* 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_ */ |
| |
