1*> \brief \b DSPR
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
12*
13*       .. Scalar Arguments ..
14*       DOUBLE PRECISION ALPHA
15*       INTEGER INCX,N
16*       CHARACTER UPLO
17*       ..
18*       .. Array Arguments ..
19*       DOUBLE PRECISION AP(*),X(*)
20*       ..
21*
22*
23*> \par Purpose:
24*  =============
25*>
26*> \verbatim
27*>
28*> DSPR    performs the symmetric rank 1 operation
29*>
30*>    A := alpha*x*x**T + A,
31*>
32*> where alpha is a real scalar, x is an n element vector and A is an
33*> n by n symmetric matrix, supplied in packed form.
34*> \endverbatim
35*
36*  Arguments:
37*  ==========
38*
39*> \param[in] UPLO
40*> \verbatim
41*>          UPLO is CHARACTER*1
42*>           On entry, UPLO specifies whether the upper or lower
43*>           triangular part of the matrix A is supplied in the packed
44*>           array AP as follows:
45*>
46*>              UPLO = 'U' or 'u'   The upper triangular part of A is
47*>                                  supplied in AP.
48*>
49*>              UPLO = 'L' or 'l'   The lower triangular part of A is
50*>                                  supplied in AP.
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*>          N is INTEGER
56*>           On entry, N specifies the order of the matrix A.
57*>           N must be at least zero.
58*> \endverbatim
59*>
60*> \param[in] ALPHA
61*> \verbatim
62*>          ALPHA is DOUBLE PRECISION.
63*>           On entry, ALPHA specifies the scalar alpha.
64*> \endverbatim
65*>
66*> \param[in] X
67*> \verbatim
68*>          X is DOUBLE PRECISION array, dimension at least
69*>           ( 1 + ( n - 1 )*abs( INCX ) ).
70*>           Before entry, the incremented array X must contain the n
71*>           element vector x.
72*> \endverbatim
73*>
74*> \param[in] INCX
75*> \verbatim
76*>          INCX is INTEGER
77*>           On entry, INCX specifies the increment for the elements of
78*>           X. INCX must not be zero.
79*> \endverbatim
80*>
81*> \param[in,out] AP
82*> \verbatim
83*>          AP is DOUBLE PRECISION array, dimension at least
84*>           ( ( n*( n + 1 ) )/2 ).
85*>           Before entry with  UPLO = 'U' or 'u', the array AP must
86*>           contain the upper triangular part of the symmetric matrix
87*>           packed sequentially, column by column, so that AP( 1 )
88*>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
89*>           and a( 2, 2 ) respectively, and so on. On exit, the array
90*>           AP is overwritten by the upper triangular part of the
91*>           updated matrix.
92*>           Before entry with UPLO = 'L' or 'l', the array AP must
93*>           contain the lower triangular part of the symmetric matrix
94*>           packed sequentially, column by column, so that AP( 1 )
95*>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
96*>           and a( 3, 1 ) respectively, and so on. On exit, the array
97*>           AP is overwritten by the lower triangular part of the
98*>           updated matrix.
99*> \endverbatim
100*
101*  Authors:
102*  ========
103*
104*> \author Univ. of Tennessee
105*> \author Univ. of California Berkeley
106*> \author Univ. of Colorado Denver
107*> \author NAG Ltd.
108*
109*> \ingroup double_blas_level2
110*
111*> \par Further Details:
112*  =====================
113*>
114*> \verbatim
115*>
116*>  Level 2 Blas routine.
117*>
118*>  -- Written on 22-October-1986.
119*>     Jack Dongarra, Argonne National Lab.
120*>     Jeremy Du Croz, Nag Central Office.
121*>     Sven Hammarling, Nag Central Office.
122*>     Richard Hanson, Sandia National Labs.
123*> \endverbatim
124*>
125*  =====================================================================
126      SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
127*
128*  -- Reference BLAS level2 routine --
129*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
130*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
131*
132*     .. Scalar Arguments ..
133      DOUBLE PRECISION ALPHA
134      INTEGER INCX,N
135      CHARACTER UPLO
136*     ..
137*     .. Array Arguments ..
138      DOUBLE PRECISION AP(*),X(*)
139*     ..
140*
141*  =====================================================================
142*
143*     .. Parameters ..
144      DOUBLE PRECISION ZERO
145      PARAMETER (ZERO=0.0D+0)
146*     ..
147*     .. Local Scalars ..
148      DOUBLE PRECISION TEMP
149      INTEGER I,INFO,IX,J,JX,K,KK,KX
150*     ..
151*     .. External Functions ..
152      LOGICAL LSAME
153      EXTERNAL LSAME
154*     ..
155*     .. External Subroutines ..
156      EXTERNAL XERBLA
157*     ..
158*
159*     Test the input parameters.
160*
161      INFO = 0
162      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
163          INFO = 1
164      ELSE IF (N.LT.0) THEN
165          INFO = 2
166      ELSE IF (INCX.EQ.0) THEN
167          INFO = 5
168      END IF
169      IF (INFO.NE.0) THEN
170          CALL XERBLA('DSPR  ',INFO)
171          RETURN
172      END IF
173*
174*     Quick return if possible.
175*
176      IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
177*
178*     Set the start point in X if the increment is not unity.
179*
180      IF (INCX.LE.0) THEN
181          KX = 1 - (N-1)*INCX
182      ELSE IF (INCX.NE.1) THEN
183          KX = 1
184      END IF
185*
186*     Start the operations. In this version the elements of the array AP
187*     are accessed sequentially with one pass through AP.
188*
189      KK = 1
190      IF (LSAME(UPLO,'U')) THEN
191*
192*        Form  A  when upper triangle is stored in AP.
193*
194          IF (INCX.EQ.1) THEN
195              DO 20 J = 1,N
196                  IF (X(J).NE.ZERO) THEN
197                      TEMP = ALPHA*X(J)
198                      K = KK
199                      DO 10 I = 1,J
200                          AP(K) = AP(K) + X(I)*TEMP
201                          K = K + 1
202   10                 CONTINUE
203                  END IF
204                  KK = KK + J
205   20         CONTINUE
206          ELSE
207              JX = KX
208              DO 40 J = 1,N
209                  IF (X(JX).NE.ZERO) THEN
210                      TEMP = ALPHA*X(JX)
211                      IX = KX
212                      DO 30 K = KK,KK + J - 1
213                          AP(K) = AP(K) + X(IX)*TEMP
214                          IX = IX + INCX
215   30                 CONTINUE
216                  END IF
217                  JX = JX + INCX
218                  KK = KK + J
219   40         CONTINUE
220          END IF
221      ELSE
222*
223*        Form  A  when lower triangle is stored in AP.
224*
225          IF (INCX.EQ.1) THEN
226              DO 60 J = 1,N
227                  IF (X(J).NE.ZERO) THEN
228                      TEMP = ALPHA*X(J)
229                      K = KK
230                      DO 50 I = J,N
231                          AP(K) = AP(K) + X(I)*TEMP
232                          K = K + 1
233   50                 CONTINUE
234                  END IF
235                  KK = KK + N - J + 1
236   60         CONTINUE
237          ELSE
238              JX = KX
239              DO 80 J = 1,N
240                  IF (X(JX).NE.ZERO) THEN
241                      TEMP = ALPHA*X(JX)
242                      IX = JX
243                      DO 70 K = KK,KK + N - J
244                          AP(K) = AP(K) + X(IX)*TEMP
245                          IX = IX + INCX
246   70                 CONTINUE
247                  END IF
248                  JX = JX + INCX
249                  KK = KK + N - J + 1
250   80         CONTINUE
251          END IF
252      END IF
253*
254      RETURN
255*
256*     End of DSPR
257*
258      END
259