/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:15 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dgemv.h" #include /* PARAMETER translations */ #define ONE 1.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dgemv( byte trans, long m, long n, double alpha, double *a, long lda, double x[], long incx, double beta, double y[], long incy) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny; double temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /*--D replaces "?": ?GEMV * .. Scalar Arguments .. */ /* .. * .. Array Arguments .. */ /* .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Arguments * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * 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 * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, 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. * */ info = 0; if ((!lsame( trans, 'N' ) && !lsame( trans, 'T' )) && !lsame( trans, 'C' )) { info = 1; } else if (m < 0) { info = 2; } else if (n < 0) { info = 3; } else if (lda < max( 1, m )) { info = 6; } else if (incx == 0) { info = 8; } else if (incy == 0) { info = 11; } if (info != 0) { xerbla( "DGEMV ", info ); return; } /* Quick return if possible. * */ if (((m == 0) || (n == 0)) || ((alpha == ZERO) && (beta == ONE))) return; /* Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * */ if (lsame( trans, 'N' )) { lenx = n; leny = m; } else { lenx = m; leny = n; } if (incx > 0) { kx = 1; } else { kx = 1 - (lenx - 1)*incx; } if (incy > 0) { ky = 1; } else { ky = 1 - (leny - 1)*incy; } /* Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * */ if (beta != ONE) { if (incy == 1) { if (beta == ZERO) { for (i = 1; i <= leny; i++) { Y[i] = ZERO; } } else { for (i = 1; i <= leny; i++) { Y[i] *= beta; } } } else { iy = ky; if (beta == ZERO) { for (i = 1; i <= leny; i++) { Y[iy] = ZERO; iy += incy; } } else { for (i = 1; i <= leny; i++) { Y[iy] *= beta; iy += incy; } } } } if (alpha == ZERO) return; if (lsame( trans, 'N' )) { /* Form y := alpha*A*x + y. * */ jx = kx; if (incy == 1) { for (j = 1; j <= n; j++) { if (X[jx] != ZERO) { temp = alpha*X[jx]; for (i = 1; i <= m; i++) { Y[i] += temp*A(j - 1,i - 1); } } jx += incx; } } else { for (j = 1; j <= n; j++) { if (X[jx] != ZERO) { temp = alpha*X[jx]; iy = ky; for (i = 1; i <= m; i++) { Y[iy] += temp*A(j - 1,i - 1); iy += incy; } } jx += incx; } } } else { /* Form y := alpha*A'*x + y. * */ jy = ky; if (incx == 1) { for (j = 1; j <= n; j++) { temp = ZERO; for (i = 1; i <= m; i++) { temp += A(j - 1,i - 1)*X[i]; } Y[jy] += alpha*temp; jy += incy; } } else { for (j = 1; j <= n; j++) { temp = ZERO; ix = kx; for (i = 1; i <= m; i++) { temp += A(j - 1,i - 1)*X[ix]; ix += incx; } Y[jy] += alpha*temp; jy += incy; } } } return; /* End of DGEMV . * */ #undef A } /* end of function */