1*> \brief \b DSYMV
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 DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
12*
13*       .. Scalar Arguments ..
14*       DOUBLE PRECISION ALPHA,BETA
15*       INTEGER INCX,INCY,LDA,N
16*       CHARACTER UPLO
17*       ..
18*       .. Array Arguments ..
19*       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
20*       ..
21*
22*
23*> \par Purpose:
24*  =============
25*>
26*> \verbatim
27*>
28*> DSYMV  performs the matrix-vector  operation
29*>
30*>    y := alpha*A*x + beta*y,
31*>
32*> where alpha and beta are scalars, x and y are n element vectors and
33*> A is an n by n symmetric matrix.
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 array A is to be referenced as
44*>           follows:
45*>
46*>              UPLO = 'U' or 'u'   Only the upper triangular part of A
47*>                                  is to be referenced.
48*>
49*>              UPLO = 'L' or 'l'   Only the lower triangular part of A
50*>                                  is to be referenced.
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] A
67*> \verbatim
68*>          A is DOUBLE PRECISION array, dimension ( LDA, N )
69*>           Before entry with  UPLO = 'U' or 'u', the leading n by n
70*>           upper triangular part of the array A must contain the upper
71*>           triangular part of the symmetric matrix and the strictly
72*>           lower triangular part of A is not referenced.
73*>           Before entry with UPLO = 'L' or 'l', the leading n by n
74*>           lower triangular part of the array A must contain the lower
75*>           triangular part of the symmetric matrix and the strictly
76*>           upper triangular part of A is not referenced.
77*> \endverbatim
78*>
79*> \param[in] LDA
80*> \verbatim
81*>          LDA is INTEGER
82*>           On entry, LDA specifies the first dimension of A as declared
83*>           in the calling (sub) program. LDA must be at least
84*>           max( 1, n ).
85*> \endverbatim
86*>
87*> \param[in] X
88*> \verbatim
89*>          X is DOUBLE PRECISION array, dimension at least
90*>           ( 1 + ( n - 1 )*abs( INCX ) ).
91*>           Before entry, the incremented array X must contain the n
92*>           element vector x.
93*> \endverbatim
94*>
95*> \param[in] INCX
96*> \verbatim
97*>          INCX is INTEGER
98*>           On entry, INCX specifies the increment for the elements of
99*>           X. INCX must not be zero.
100*> \endverbatim
101*>
102*> \param[in] BETA
103*> \verbatim
104*>          BETA is DOUBLE PRECISION.
105*>           On entry, BETA specifies the scalar beta. When BETA is
106*>           supplied as zero then Y need not be set on input.
107*> \endverbatim
108*>
109*> \param[in,out] Y
110*> \verbatim
111*>          Y is DOUBLE PRECISION array, dimension at least
112*>           ( 1 + ( n - 1 )*abs( INCY ) ).
113*>           Before entry, the incremented array Y must contain the n
114*>           element vector y. On exit, Y is overwritten by the updated
115*>           vector y.
116*> \endverbatim
117*>
118*> \param[in] INCY
119*> \verbatim
120*>          INCY is INTEGER
121*>           On entry, INCY specifies the increment for the elements of
122*>           Y. INCY must not be zero.
123*> \endverbatim
124*
125*  Authors:
126*  ========
127*
128*> \author Univ. of Tennessee
129*> \author Univ. of California Berkeley
130*> \author Univ. of Colorado Denver
131*> \author NAG Ltd.
132*
133*> \ingroup double_blas_level2
134*
135*> \par Further Details:
136*  =====================
137*>
138*> \verbatim
139*>
140*>  Level 2 Blas routine.
141*>  The vector and matrix arguments are not referenced when N = 0, or M = 0
142*>
143*>  -- Written on 22-October-1986.
144*>     Jack Dongarra, Argonne National Lab.
145*>     Jeremy Du Croz, Nag Central Office.
146*>     Sven Hammarling, Nag Central Office.
147*>     Richard Hanson, Sandia National Labs.
148*> \endverbatim
149*>
150*  =====================================================================
151      SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
152*
153*  -- Reference BLAS level2 routine --
154*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
155*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156*
157*     .. Scalar Arguments ..
158      DOUBLE PRECISION ALPHA,BETA
159      INTEGER INCX,INCY,LDA,N
160      CHARACTER UPLO
161*     ..
162*     .. Array Arguments ..
163      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
164*     ..
165*
166*  =====================================================================
167*
168*     .. Parameters ..
169      DOUBLE PRECISION ONE,ZERO
170      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
171*     ..
172*     .. Local Scalars ..
173      DOUBLE PRECISION TEMP1,TEMP2
174      INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY
175*     ..
176*     .. External Functions ..
177      LOGICAL LSAME
178      EXTERNAL LSAME
179*     ..
180*     .. External Subroutines ..
181      EXTERNAL XERBLA
182*     ..
183*     .. Intrinsic Functions ..
184      INTRINSIC MAX
185*     ..
186*
187*     Test the input parameters.
188*
189      INFO = 0
190      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
191          INFO = 1
192      ELSE IF (N.LT.0) THEN
193          INFO = 2
194      ELSE IF (LDA.LT.MAX(1,N)) THEN
195          INFO = 5
196      ELSE IF (INCX.EQ.0) THEN
197          INFO = 7
198      ELSE IF (INCY.EQ.0) THEN
199          INFO = 10
200      END IF
201      IF (INFO.NE.0) THEN
202          CALL XERBLA('DSYMV ',INFO)
203          RETURN
204      END IF
205*
206*     Quick return if possible.
207*
208      IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
209*
210*     Set up the start points in  X  and  Y.
211*
212      IF (INCX.GT.0) THEN
213          KX = 1
214      ELSE
215          KX = 1 - (N-1)*INCX
216      END IF
217      IF (INCY.GT.0) THEN
218          KY = 1
219      ELSE
220          KY = 1 - (N-1)*INCY
221      END IF
222*
223*     Start the operations. In this version the elements of A are
224*     accessed sequentially with one pass through the triangular part
225*     of A.
226*
227*     First form  y := beta*y.
228*
229      IF (BETA.NE.ONE) THEN
230          IF (INCY.EQ.1) THEN
231              IF (BETA.EQ.ZERO) THEN
232                  DO 10 I = 1,N
233                      Y(I) = ZERO
234   10             CONTINUE
235              ELSE
236                  DO 20 I = 1,N
237                      Y(I) = BETA*Y(I)
238   20             CONTINUE
239              END IF
240          ELSE
241              IY = KY
242              IF (BETA.EQ.ZERO) THEN
243                  DO 30 I = 1,N
244                      Y(IY) = ZERO
245                      IY = IY + INCY
246   30             CONTINUE
247              ELSE
248                  DO 40 I = 1,N
249                      Y(IY) = BETA*Y(IY)
250                      IY = IY + INCY
251   40             CONTINUE
252              END IF
253          END IF
254      END IF
255      IF (ALPHA.EQ.ZERO) RETURN
256      IF (LSAME(UPLO,'U')) THEN
257*
258*        Form  y  when A is stored in upper triangle.
259*
260          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
261              DO 60 J = 1,N
262                  TEMP1 = ALPHA*X(J)
263                  TEMP2 = ZERO
264                  DO 50 I = 1,J - 1
265                      Y(I) = Y(I) + TEMP1*A(I,J)
266                      TEMP2 = TEMP2 + A(I,J)*X(I)
267   50             CONTINUE
268                  Y(J) = Y(J) + TEMP1*A(J,J) + ALPHA*TEMP2
269   60         CONTINUE
270          ELSE
271              JX = KX
272              JY = KY
273              DO 80 J = 1,N
274                  TEMP1 = ALPHA*X(JX)
275                  TEMP2 = ZERO
276                  IX = KX
277                  IY = KY
278                  DO 70 I = 1,J - 1
279                      Y(IY) = Y(IY) + TEMP1*A(I,J)
280                      TEMP2 = TEMP2 + A(I,J)*X(IX)
281                      IX = IX + INCX
282                      IY = IY + INCY
283   70             CONTINUE
284                  Y(JY) = Y(JY) + TEMP1*A(J,J) + ALPHA*TEMP2
285                  JX = JX + INCX
286                  JY = JY + INCY
287   80         CONTINUE
288          END IF
289      ELSE
290*
291*        Form  y  when A is stored in lower triangle.
292*
293          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
294              DO 100 J = 1,N
295                  TEMP1 = ALPHA*X(J)
296                  TEMP2 = ZERO
297                  Y(J) = Y(J) + TEMP1*A(J,J)
298                  DO 90 I = J + 1,N
299                      Y(I) = Y(I) + TEMP1*A(I,J)
300                      TEMP2 = TEMP2 + A(I,J)*X(I)
301   90             CONTINUE
302                  Y(J) = Y(J) + ALPHA*TEMP2
303  100         CONTINUE
304          ELSE
305              JX = KX
306              JY = KY
307              DO 120 J = 1,N
308                  TEMP1 = ALPHA*X(JX)
309                  TEMP2 = ZERO
310                  Y(JY) = Y(JY) + TEMP1*A(J,J)
311                  IX = JX
312                  IY = JY
313                  DO 110 I = J + 1,N
314                      IX = IX + INCX
315                      IY = IY + INCY
316                      Y(IY) = Y(IY) + TEMP1*A(I,J)
317                      TEMP2 = TEMP2 + A(I,J)*X(IX)
318  110             CONTINUE
319                  Y(JY) = Y(JY) + ALPHA*TEMP2
320                  JX = JX + INCX
321                  JY = JY + INCY
322  120         CONTINUE
323          END IF
324      END IF
325*
326      RETURN
327*
328*     End of DSYMV
329*
330      END
331