share/lapack/blas/fortran/dgemv.f

   1       SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
   2      $                   BETA, Y, INCY )
   3 *     .. Scalar Arguments ..
   4       DOUBLE PRECISION   ALPHA, BETA
   5       INTEGER            INCX, INCY, LDA, M, N
   6       CHARACTER*1        TRANS
   7 *     .. Array Arguments ..
   8       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  DGEMV  performs one of the matrix-vector operations
  15 *
  16 *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
  17 *
  18 *  where alpha and beta are scalars, x and y are vectors and A is an
  19 *  m by n matrix.
  20 *
  21 *  Parameters
  22 *  ==========
  23 *
  24 *  TRANS  - CHARACTER*1.
  25 *           On entry, TRANS specifies the operation to be performed as
  26 *           follows:
  27 *
  28 *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  29 *
  30 *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  31 *
  32 *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
  33 *
  34 *           Unchanged on exit.
  35 *
  36 *  M      - INTEGER.
  37 *           On entry, M specifies the number of rows of the matrix A.
  38 *           M must be at least zero.
  39 *           Unchanged on exit.
  40 *
  41 *  N      - INTEGER.
  42 *           On entry, N specifies the number of columns of the matrix A.
  43 *           N must be at least zero.
  44 *           Unchanged on exit.
  45 *
  46 *  ALPHA  - DOUBLE PRECISION.
  47 *           On entry, ALPHA specifies the scalar alpha.
  48 *           Unchanged on exit.
  49 *
  50 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  51 *           Before entry, the leading m by n part of the array A must
  52 *           contain the matrix of coefficients.
  53 *           Unchanged on exit.
  54 *
  55 *  LDA    - INTEGER.
  56 *           On entry, LDA specifies the first dimension of A as declared
  57 *           in the calling (sub) program. LDA must be at least
  58 *           max( 1, m ).
  59 *           Unchanged on exit.
  60 *
  61 *  X      - DOUBLE PRECISION array of DIMENSION at least
  62 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  63 *           and at least
  64 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  65 *           Before entry, the incremented array X must contain the
  66 *           vector x.
  67 *           Unchanged on exit.
  68 *
  69 *  INCX   - INTEGER.
  70 *           On entry, INCX specifies the increment for the elements of
  71 *           X. INCX must not be zero.
  72 *           Unchanged on exit.
  73 *
  74 *  BETA   - DOUBLE PRECISION.
  75 *           On entry, BETA specifies the scalar beta. When BETA is
  76 *           supplied as zero then Y need not be set on input.
  77 *           Unchanged on exit.
  78 *
  79 *  Y      - DOUBLE PRECISION array of DIMENSION at least
  80 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  81 *           and at least
  82 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  83 *           Before entry with BETA non-zero, the incremented array Y
  84 *           must contain the vector y. On exit, Y is overwritten by the
  85 *           updated vector y.
  86 *
  87 *  INCY   - INTEGER.
  88 *           On entry, INCY specifies the increment for the elements of
  89 *           Y. INCY must not be zero.
  90 *           Unchanged on exit.
  91 *
  92 *
  93 *  Level 2 Blas routine.
  94 *
  95 *  -- Written on 22-October-1986.
  96 *     Jack Dongarra, Argonne National Lab.
  97 *     Jeremy Du Croz, Nag Central Office.
  98 *     Sven Hammarling, Nag Central Office.
  99 *     Richard Hanson, Sandia National Labs.
 100 *
 101 *
 102 *     .. Parameters ..
 103       DOUBLE PRECISION   ONE         , ZERO
 104       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 105 *     .. Local Scalars ..
 106       DOUBLE PRECISION   TEMP
 107       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
 108 *     .. External Functions ..
 109       LOGICAL            LSAME
 110       EXTERNAL           LSAME
 111 *     .. External Subroutines ..
 112       EXTERNAL           XERBLA
 113 *     .. Intrinsic Functions ..
 114       INTRINSIC          MAX
 115 *     ..
 116 *     .. Executable Statements ..
 117 *
 118 *     Test the input parameters.
 119 *
 120       INFO = 0
 121       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
 122      $         .NOT.LSAME( TRANS, 'T' ).AND.
 123      $         .NOT.LSAME( TRANS, 'C' )      )THEN
 124          INFO = 1
 125       ELSE IF( M.LT.0 )THEN
 126          INFO = 2
 127       ELSE IF( N.LT.0 )THEN
 128          INFO = 3
 129       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
 130          INFO = 6
 131       ELSE IF( INCX.EQ.0 )THEN
 132          INFO = 8
 133       ELSE IF( INCY.EQ.0 )THEN
 134          INFO = 11
 135       END IF
 136       IF( INFO.NE.0 )THEN
 137          CALL XERBLA( 'DGEMV ', INFO )
 138          RETURN
 139       END IF
 140 *
 141 *     Quick return if possible.
 142 *
 143       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
 144      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 145      $   RETURN
 146 *
 147 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
 148 *     up the start points in  X  and  Y.
 149 *
 150       IF( LSAME( TRANS, 'N' ) )THEN
 151          LENX = N
 152          LENY = M
 153       ELSE
 154          LENX = M
 155          LENY = N
 156       END IF
 157       IF( INCX.GT.0 )THEN
 158          KX = 1
 159       ELSE
 160          KX = 1 - ( LENX - 1 )*INCX
 161       END IF
 162       IF( INCY.GT.0 )THEN
 163          KY = 1
 164       ELSE
 165          KY = 1 - ( LENY - 1 )*INCY
 166       END IF
 167 *
 168 *     Start the operations. In this version the elements of A are
 169 *     accessed sequentially with one pass through A.
 170 *
 171 *     First form  y := beta*y.
 172 *
 173       IF( BETA.NE.ONE )THEN
 174          IF( INCY.EQ.1 )THEN
 175             IF( BETA.EQ.ZERO )THEN
 176                DO 10, I = 1, LENY
 177                   Y( I ) = ZERO
 178    10          CONTINUE
 179             ELSE
 180                DO 20, I = 1, LENY
 181                   Y( I ) = BETA*Y( I )
 182    20          CONTINUE
 183             END IF
 184          ELSE
 185             IY = KY
 186             IF( BETA.EQ.ZERO )THEN
 187                DO 30, I = 1, LENY
 188                   Y( IY ) = ZERO
 189                   IY      = IY   + INCY
 190    30          CONTINUE
 191             ELSE
 192                DO 40, I = 1, LENY
 193                   Y( IY ) = BETA*Y( IY )
 194                   IY      = IY           + INCY
 195    40          CONTINUE
 196             END IF
 197          END IF
 198       END IF
 199       IF( ALPHA.EQ.ZERO )
 200      $   RETURN
 201       IF( LSAME( TRANS, 'N' ) )THEN
 202 *
 203 *        Form  y := alpha*A*x + y.
 204 *
 205          JX = KX
 206          IF( INCY.EQ.1 )THEN
 207             DO 60, J = 1, N
 208                IF( X( JX ).NE.ZERO )THEN
 209                   TEMP = ALPHA*X( JX )
 210                   DO 50, I = 1, M
 211                      Y( I ) = Y( I ) + TEMP*A( I, J )
 212    50             CONTINUE
 213                END IF
 214                JX = JX + INCX
 215    60       CONTINUE
 216          ELSE
 217             DO 80, J = 1, N
 218                IF( X( JX ).NE.ZERO )THEN
 219                   TEMP = ALPHA*X( JX )
 220                   IY   = KY
 221                   DO 70, I = 1, M
 222                      Y( IY ) = Y( IY ) + TEMP*A( I, J )
 223                      IY      = IY      + INCY
 224    70             CONTINUE
 225                END IF
 226                JX = JX + INCX
 227    80       CONTINUE
 228          END IF
 229       ELSE
 230 *
 231 *        Form  y := alpha*A'*x + y.
 232 *
 233          JY = KY
 234          IF( INCX.EQ.1 )THEN
 235             DO 100, J = 1, N
 236                TEMP = ZERO
 237                DO 90, I = 1, M
 238                   TEMP = TEMP + A( I, J )*X( I )
 239    90          CONTINUE
 240                Y( JY ) = Y( JY ) + ALPHA*TEMP
 241                JY      = JY      + INCY
 242   100       CONTINUE
 243          ELSE
 244             DO 120, J = 1, N
 245                TEMP = ZERO
 246                IX   = KX
 247                DO 110, I = 1, M
 248                   TEMP = TEMP + A( I, J )*X( IX )
 249                   IX   = IX   + INCX
 250   110          CONTINUE
 251                Y( JY ) = Y( JY ) + ALPHA*TEMP
 252                JY      = JY      + INCY
 253   120       CONTINUE
 254          END IF
 255       END IF
 256 *
 257       RETURN
 258 *
 259 *     End of DGEMV .
 260 *
 261       END