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