2
************************************************************************
4
SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
6
* .. Scalar Arguments ..
7
DOUBLE PRECISION ALPHA, BETA
8
INTEGER INCX, INCY, LDA, N
10
* .. Array Arguments ..
11
DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
17
* DSYMV performs the matrix-vector operation
19
* y := alpha*A*x + beta*y,
21
* where alpha and beta are scalars, x and y are n element vectors and
22
* A is an n by n symmetric matrix.
28
* On entry, UPLO specifies whether the upper or lower
29
* triangular part of the array A is to be referenced as
32
* UPLO = 'U' or 'u' Only the upper triangular part of A
33
* is to be referenced.
35
* UPLO = 'L' or 'l' Only the lower triangular part of A
36
* is to be referenced.
41
* On entry, N specifies the order of the matrix A.
42
* N must be at least zero.
45
* ALPHA - DOUBLE PRECISION.
46
* On entry, ALPHA specifies the scalar alpha.
49
* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
50
* Before entry with UPLO = 'U' or 'u', the leading n by n
51
* upper triangular part of the array A must contain the upper
52
* triangular part of the symmetric matrix and the strictly
53
* lower triangular part of A is not referenced.
54
* Before entry with UPLO = 'L' or 'l', the leading n by n
55
* lower triangular part of the array A must contain the lower
56
* triangular part of the symmetric matrix and the strictly
57
* upper triangular part of A is not referenced.
61
* On entry, LDA specifies the first dimension of A as declared
62
* in the calling (sub) program. LDA must be at least
66
* X - DOUBLE PRECISION array of dimension at least
67
* ( 1 + ( n - 1 )*abs( INCX ) ).
68
* Before entry, the incremented array X must contain the n
73
* On entry, INCX specifies the increment for the elements of
74
* X. INCX must not be zero.
77
* BETA - DOUBLE PRECISION.
78
* On entry, BETA specifies the scalar beta. When BETA is
79
* supplied as zero then Y need not be set on input.
82
* Y - DOUBLE PRECISION array of dimension at least
83
* ( 1 + ( n - 1 )*abs( INCY ) ).
84
* Before entry, the incremented array Y must contain the n
85
* element vector y. On exit, Y is overwritten by the updated
89
* On entry, INCY specifies the increment for the elements of
90
* Y. INCY must not be zero.
94
* Level 2 Blas routine.
96
* -- Written on 22-October-1986.
97
* Jack Dongarra, Argonne National Lab.
98
* Jeremy Du Croz, Nag Central Office.
99
* Sven Hammarling, Nag Central Office.
100
* Richard Hanson, Sandia National Labs.
104
DOUBLE PRECISION ONE , ZERO
105
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
106
* .. Local Scalars ..
107
DOUBLE PRECISION TEMP1, TEMP2
108
INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY
109
* .. External Functions ..
112
* .. External Subroutines ..
114
* .. Intrinsic Functions ..
117
* .. Executable Statements ..
119
* Test the input parameters.
122
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
123
$ .NOT.LSAME( UPLO, 'L' ) )THEN
125
ELSE IF( N.LT.0 )THEN
127
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
129
ELSE IF( INCX.EQ.0 )THEN
131
ELSE IF( INCY.EQ.0 )THEN
135
CALL XERBLA( 'DSYMV ', INFO )
139
* Quick return if possible.
141
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
144
* Set up the start points in X and Y.
149
KX = 1 - ( N - 1 )*INCX
154
KY = 1 - ( N - 1 )*INCY
157
* Start the operations. In this version the elements of A are
158
* accessed sequentially with one pass through the triangular part
161
* First form y := beta*y.
163
IF( BETA.NE.ONE )THEN
165
IF( BETA.EQ.ZERO )THEN
176
IF( BETA.EQ.ZERO )THEN
183
Y( IY ) = BETA*Y( IY )
191
IF( LSAME( UPLO, 'U' ) )THEN
193
* Form y when A is stored in upper triangle.
195
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
200
Y( I ) = Y( I ) + TEMP1*A( I, J )
201
TEMP2 = TEMP2 + A( I, J )*X( I )
203
Y( J ) = Y( J ) + TEMP1*A( J, J ) + ALPHA*TEMP2
209
TEMP1 = ALPHA*X( JX )
214
Y( IY ) = Y( IY ) + TEMP1*A( I, J )
215
TEMP2 = TEMP2 + A( I, J )*X( IX )
219
Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2
226
* Form y when A is stored in lower triangle.
228
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
232
Y( J ) = Y( J ) + TEMP1*A( J, J )
234
Y( I ) = Y( I ) + TEMP1*A( I, J )
235
TEMP2 = TEMP2 + A( I, J )*X( I )
237
Y( J ) = Y( J ) + ALPHA*TEMP2
243
TEMP1 = ALPHA*X( JX )
245
Y( JY ) = Y( JY ) + TEMP1*A( J, J )
251
Y( IY ) = Y( IY ) + TEMP1*A( I, J )
252
TEMP2 = TEMP2 + A( I, J )*X( IX )
254
Y( JY ) = Y( JY ) + ALPHA*TEMP2