var/external/blas/dsymv.inc

   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 SUBROUTINE DSYMV