1
SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
2
* .. Scalar Arguments ..
3
DOUBLE PRECISION ALPHA, BETA
6
* .. Array Arguments ..
7
DOUBLE PRECISION AP( * ), X( * ), Y( * )
13
* DSPMV performs the matrix-vector operation
15
* y := alpha*A*x + beta*y,
17
* where alpha and beta are scalars, x and y are n element vectors and
18
* A is an n by n symmetric matrix, supplied in packed form.
24
* On entry, UPLO specifies whether the upper or lower
25
* triangular part of the matrix A is supplied in the packed
26
* array AP as follows:
28
* UPLO = 'U' or 'u' The upper triangular part of A is
31
* UPLO = 'L' or 'l' The lower triangular part of A is
37
* On entry, N specifies the order of the matrix A.
38
* N must be at least zero.
41
* ALPHA - DOUBLE PRECISION.
42
* On entry, ALPHA specifies the scalar alpha.
45
* AP - DOUBLE PRECISION array of DIMENSION at least
46
* ( ( n*( n + 1 ) )/2 ).
47
* Before entry with UPLO = 'U' or 'u', the array AP must
48
* contain the upper triangular part of the symmetric matrix
49
* packed sequentially, column by column, so that AP( 1 )
50
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
51
* and a( 2, 2 ) respectively, and so on.
52
* Before entry with UPLO = 'L' or 'l', the array AP must
53
* contain the lower triangular part of the symmetric matrix
54
* packed sequentially, column by column, so that AP( 1 )
55
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
56
* and a( 3, 1 ) respectively, and so on.
59
* X - DOUBLE PRECISION array of dimension at least
60
* ( 1 + ( n - 1 )*abs( INCX ) ).
61
* Before entry, the incremented array X must contain the n
66
* On entry, INCX specifies the increment for the elements of
67
* X. INCX must not be zero.
70
* BETA - DOUBLE PRECISION.
71
* On entry, BETA specifies the scalar beta. When BETA is
72
* supplied as zero then Y need not be set on input.
75
* Y - DOUBLE PRECISION array of dimension at least
76
* ( 1 + ( n - 1 )*abs( INCY ) ).
77
* Before entry, the incremented array Y must contain the n
78
* element vector y. On exit, Y is overwritten by the updated
82
* On entry, INCY specifies the increment for the elements of
83
* Y. INCY must not be zero.
87
* Level 2 Blas routine.
89
* -- Written on 22-October-1986.
90
* Jack Dongarra, Argonne National Lab.
91
* Jeremy Du Croz, Nag Central Office.
92
* Sven Hammarling, Nag Central Office.
93
* Richard Hanson, Sandia National Labs.
97
DOUBLE PRECISION ONE , ZERO
98
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
100
DOUBLE PRECISION TEMP1, TEMP2
101
INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
102
* .. External Functions ..
105
* .. External Subroutines ..
108
* .. Executable Statements ..
110
* Test the input parameters.
113
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
114
$ .NOT.LSAME( UPLO, 'L' ) )THEN
116
ELSE IF( N.LT.0 )THEN
118
ELSE IF( INCX.EQ.0 )THEN
120
ELSE IF( INCY.EQ.0 )THEN
124
CALL XERBLA( 'DSPMV ', INFO )
128
* Quick return if possible.
130
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
133
* Set up the start points in X and Y.
138
KX = 1 - ( N - 1 )*INCX
143
KY = 1 - ( N - 1 )*INCY
146
* Start the operations. In this version the elements of the array AP
147
* are accessed sequentially with one pass through AP.
149
* First form y := beta*y.
151
IF( BETA.NE.ONE )THEN
153
IF( BETA.EQ.ZERO )THEN
164
IF( BETA.EQ.ZERO )THEN
171
Y( IY ) = BETA*Y( IY )
180
IF( LSAME( UPLO, 'U' ) )THEN
182
* Form y when AP contains the upper triangle.
184
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
190
Y( I ) = Y( I ) + TEMP1*AP( K )
191
TEMP2 = TEMP2 + AP( K )*X( I )
194
Y( J ) = Y( J ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2
201
TEMP1 = ALPHA*X( JX )
205
DO 70, K = KK, KK + J - 2
206
Y( IY ) = Y( IY ) + TEMP1*AP( K )
207
TEMP2 = TEMP2 + AP( K )*X( IX )
211
Y( JY ) = Y( JY ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2
219
* Form y when AP contains the lower triangle.
221
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
225
Y( J ) = Y( J ) + TEMP1*AP( KK )
228
Y( I ) = Y( I ) + TEMP1*AP( K )
229
TEMP2 = TEMP2 + AP( K )*X( I )
232
Y( J ) = Y( J ) + ALPHA*TEMP2
233
KK = KK + ( N - J + 1 )
239
TEMP1 = ALPHA*X( JX )
241
Y( JY ) = Y( JY ) + TEMP1*AP( KK )
244
DO 110, K = KK + 1, KK + N - J
247
Y( IY ) = Y( IY ) + TEMP1*AP( K )
248
TEMP2 = TEMP2 + AP( K )*X( IX )
250
Y( JY ) = Y( JY ) + ALPHA*TEMP2
253
KK = KK + ( N - J + 1 )