ground/gcs/src/libs/eigen/blas/dspmv.f

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