src/BLAS/dspr.f

   1       SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
   2 *     .. Scalar Arguments ..
   3       DOUBLE PRECISION ALPHA
   4       INTEGER INCX,N
   5       CHARACTER UPLO
   6 *     ..
   7 *     .. Array Arguments ..
   8       DOUBLE PRECISION AP(*),X(*)
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  DSPR    performs the symmetric rank 1 operation
  15 *
  16 *     A := alpha*x*x' + A,
  17 *
  18 *  where alpha is a real scalar, x is an n element vector and A is an
  19 *  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 *  X      - DOUBLE PRECISION array of dimension at least
  47 *           ( 1 + ( n - 1 )*abs( INCX ) ).
  48 *           Before entry, the incremented array X must contain the n
  49 *           element vector x.
  50 *           Unchanged on exit.
  51 *
  52 *  INCX   - INTEGER.
  53 *           On entry, INCX specifies the increment for the elements of
  54 *           X. INCX must not be zero.
  55 *           Unchanged on exit.
  56 *
  57 *  AP     - DOUBLE PRECISION array of DIMENSION at least
  58 *           ( ( n*( n + 1 ) )/2 ).
  59 *           Before entry with  UPLO = 'U' or 'u', the array AP must
  60 *           contain the upper triangular part of the symmetric matrix
  61 *           packed sequentially, column by column, so that AP( 1 )
  62 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  63 *           and a( 2, 2 ) respectively, and so on. On exit, the array
  64 *           AP is overwritten by the upper triangular part of the
  65 *           updated matrix.
  66 *           Before entry with UPLO = 'L' or 'l', the array AP must
  67 *           contain the lower triangular part of the symmetric matrix
  68 *           packed sequentially, column by column, so that AP( 1 )
  69 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  70 *           and a( 3, 1 ) respectively, and so on. On exit, the array
  71 *           AP is overwritten by the lower triangular part of the
  72 *           updated matrix.
  73 *
  74 *
  75 *  Level 2 Blas routine.
  76 *
  77 *  -- Written on 22-October-1986.
  78 *     Jack Dongarra, Argonne National Lab.
  79 *     Jeremy Du Croz, Nag Central Office.
  80 *     Sven Hammarling, Nag Central Office.
  81 *     Richard Hanson, Sandia National Labs.
  82 *
  83 *
  84 *     .. Parameters ..
  85       DOUBLE PRECISION ZERO
  86       PARAMETER (ZERO=0.0D+0)
  87 *     ..
  88 *     .. Local Scalars ..
  89       DOUBLE PRECISION TEMP
  90       INTEGER I,INFO,IX,J,JX,K,KK,KX
  91 *     ..
  92 *     .. External Functions ..
  93       LOGICAL LSAME
  94       EXTERNAL LSAME
  95 *     ..
  96 *     .. External Subroutines ..
  97       EXTERNAL XERBLA
  98 *     ..
  99 *
 100 *     Test the input parameters.
 101 *
 102       INFO = 0
 103       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
 104           INFO = 1
 105       ELSE IF (N.LT.0) THEN
 106           INFO = 2
 107       ELSE IF (INCX.EQ.0) THEN
 108           INFO = 5
 109       END IF
 110       IF (INFO.NE.0) THEN
 111           CALL XERBLA('DSPR  ',INFO)
 112           RETURN
 113       END IF
 114 *
 115 *     Quick return if possible.
 116 *
 117       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
 118 *
 119 *     Set the start point in X if the increment is not unity.
 120 *
 121       IF (INCX.LE.0) THEN
 122           KX = 1 - (N-1)*INCX
 123       ELSE IF (INCX.NE.1) THEN
 124           KX = 1
 125       END IF
 126 *
 127 *     Start the operations. In this version the elements of the array AP
 128 *     are accessed sequentially with one pass through AP.
 129 *
 130       KK = 1
 131       IF (LSAME(UPLO,'U')) THEN
 132 *
 133 *        Form  A  when upper triangle is stored in AP.
 134 *
 135           IF (INCX.EQ.1) THEN
 136               DO 20 J = 1,N
 137                   IF (X(J).NE.ZERO) THEN
 138                       TEMP = ALPHA*X(J)
 139                       K = KK
 140                       DO 10 I = 1,J
 141                           AP(K) = AP(K) + X(I)*TEMP
 142                           K = K + 1
 143    10                 CONTINUE
 144                   END IF
 145                   KK = KK + J
 146    20         CONTINUE
 147           ELSE
 148               JX = KX
 149               DO 40 J = 1,N
 150                   IF (X(JX).NE.ZERO) THEN
 151                       TEMP = ALPHA*X(JX)
 152                       IX = KX
 153                       DO 30 K = KK,KK + J - 1
 154                           AP(K) = AP(K) + X(IX)*TEMP
 155                           IX = IX + INCX
 156    30                 CONTINUE
 157                   END IF
 158                   JX = JX + INCX
 159                   KK = KK + J
 160    40         CONTINUE
 161           END IF
 162       ELSE
 163 *
 164 *        Form  A  when lower triangle is stored in AP.
 165 *
 166           IF (INCX.EQ.1) THEN
 167               DO 60 J = 1,N
 168                   IF (X(J).NE.ZERO) THEN
 169                       TEMP = ALPHA*X(J)
 170                       K = KK
 171                       DO 50 I = J,N
 172                           AP(K) = AP(K) + X(I)*TEMP
 173                           K = K + 1
 174    50                 CONTINUE
 175                   END IF
 176                   KK = KK + N - J + 1
 177    60         CONTINUE
 178           ELSE
 179               JX = KX
 180               DO 80 J = 1,N
 181                   IF (X(JX).NE.ZERO) THEN
 182                       TEMP = ALPHA*X(JX)
 183                       IX = JX
 184                       DO 70 K = KK,KK + N - J
 185                           AP(K) = AP(K) + X(IX)*TEMP
 186                           IX = IX + INCX
 187    70                 CONTINUE
 188                   END IF
 189                   JX = JX + INCX
 190                   KK = KK + N - J + 1
 191    80         CONTINUE
 192           END IF
 193       END IF
 194 *
 195       RETURN
 196 *
 197 *     End of DSPR  .
 198 *
 199       END