share/lapack/blas/fortran/dsymv.f

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