SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) **************************************************************************** * * * DATA PARALLEL BLAS based on MPL * * * * Version 1.0 1/9-92 , * * For MasPar MP-1 computers * * * * para//ab, University of Bergen, NORWAY * * * * These programs must be called using F90 style array syntax. * * Note that the F77 style calling sequence has been retained * * in this version for compatibility reasons, be aware that * * parameters related to the array dimensions and shape therefore may * * be redundant and without any influence. * * The calling sequence may be changed in a future version. * * Please report any BUGs, ideas for improvement or other * * comments to * * adm@parallab.uib.no * * * * Future versions may then reflect your suggestions. * * The most current version of this software is available * * from netlib@nac.no , send the message `send index from maspar' * * * * REVISIONS: * * * **************************************************************************** implicit none * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, KL, KU, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION, array(:,:) :: A DOUBLE PRECISION, array(:) :: x, y intent(in) :: a, x intent(inout) :: y * .. * * Purpose * ======= * * DGBMV 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 band matrix, with kl sub-diagonals and ku super-diagonals. * * Parameters * ========== * * 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. * * KL - INTEGER. * On entry, KL specifies the number of sub-diagonals of the * matrix A. KL must satisfy 0 .le. KL. * Unchanged on exit. * * KU - INTEGER. * On entry, KU specifies the number of super-diagonals of the * matrix A. KU must satisfy 0 .le. KU. * 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 ( kl + ku + 1 ) by n part of the * array A must contain the matrix of coefficients, supplied * column by column, with the leading diagonal of the matrix in * row ( ku + 1 ) of the array, the first super-diagonal * starting at position 2 in row ku, the first sub-diagonal * starting at position 1 in row ( ku + 2 ), and so on. * Elements in the array A that do not correspond to elements * in the band matrix (such as the top left ku by ku triangle) * are not referenced. * The following program segment will transfer a band matrix * from conventional full matrix storage to band storage: * * DO 20, J = 1, N * K = KU + 1 - J * DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL ) * A( K + 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 * ( kl + ku + 1 ). * 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, 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 .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER m1, p1 INTEGER I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY, $ LENX, LENY, ndiags * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX, MIN INTRINSIC eoshift INTRINSIC spread INTRINSIC sum * .. Local Arrays .. DOUBLE PRECISION, array(m) :: yloc DOUBLE PRECISION, array(n) :: xloc DOUBLE PRECISION, array(ku+kl+1, max(m,n)) :: xb * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( KL.LT.0 )THEN INFO = 4 ELSE IF( KU.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN INFO = 8 ELSE IF( INCX.EQ.0 )THEN INFO = 10 ELSE IF( INCY.EQ.0 )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGBMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN m1 = -1 p1 = 1 * * 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' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the band part of A. * ndiags = ku + kl + 1 * IF ( LSAME(TRANS, 'N') ) THEN if (beta .ne. one) then if( incy.eq.1 )then if( beta.eq.zero )then yloc = zero else yloc = beta * y(1:leny) end if else if( beta.eq.zero )then yloc = zero else yloc = beta * y(ky : ky+incy*(leny-1) : incy) end if end if else if( incy.eq.1 )then yloc = y(1:leny) else yloc = y(ky : ky+incy*(leny-1) : incy) end if endif * if( incx.eq.1 )then xloc = x(1:lenx) else xloc = x(kx : kx+incx*(lenx-1) : incx) end if xloc = xloc * alpha * xb = zero xb(:, 1:n) = spread(xloc, 1, ndiags) xb(:, 1:n) = xb(:, 1:n) * a(1:ndiags, :) do i = ku, 1, -1 CPB xb(1:i,:) = eoshift( xb(1:i,:), 2, 1 ) xb(1:i,:) = eoshift( xb(1:i,:), p1, zero, 2 ) enddo do i = ku+2, ndiags CPB xb(i:ndiags,:) = eoshift( xb(i:ndiags,:), 2, -1 ) xb(i:ndiags,:) = eoshift( xb(i:ndiags,:), m1, zero, 2 ) enddo * yloc = yloc + sum(xb(:, 1:m), 1) * IF( INCY.EQ.1 )THEN Y(1:LENY) = YLOC ELSE Y(KY : KY+INCY*(LENY-1) : INCY) = YLOC END IF ELSE * * compute A'x + beta y * if (beta .ne. one) then if( incy.eq.1 )then if( beta.eq.zero )then xloc = zero else xloc = beta * y(1:leny) end if else if( beta.eq.zero )then xloc = zero else xloc = beta * y(ky:ky+incy*(leny-1):incy) end if end if else if( incy.eq.1 )then xloc = y(1:leny) else xloc = y(ky : ky+incy*(leny-1) : incy) end if endif * if( incx.eq.1 )then yloc = x(1:lenx) else yloc = x(kx : kx+incx*(lenx-1) : incx) end if yloc = yloc * alpha * xb = zero xb(:, 1:m) = spread(yloc, 1, ndiags) do i = ku, 1, -1 CPB xb(1:i,:) = eoshift( xb(1:i,:), 2, -1 ) xb(1:i,:) = eoshift( xb(1:i,:), m1, zero, 2 ) enddo do i = ku+2, ndiags CPB xb(i:ndiags,:) = eoshift( xb(i:ndiags,:), 2, 1 ) xb(i:ndiags,:) = eoshift( xb(i:ndiags,:), p1, zero, 2 ) enddo * xb(:, 1:n) = xb(:, 1:n) * a(1:ndiags, :) xloc = xloc + sum( xb(:, 1:n), 1) * if( incy.eq.1 )then y(1:leny) = xloc else y(ky : ky+incy*(leny-1) : incy) = xloc endif ENDIF * RETURN * * End of DGBMV . * END