1 SUBROUTINE DSPMVF( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) 2* .. Scalar Arguments .. 3 DOUBLE PRECISION ALPHA, BETA 4 INTEGER INCX, INCY, N 5 CHARACTER*1 UPLO 6* .. Array Arguments .. 7 DOUBLE PRECISION AP( * ), X( * ), Y( * ) 8* .. 9* 10* Purpose 11* ======= 12* 13* DSPMV performs the matrix-vector operation 14* 15* y := alpha*A*x + beta*y, 16* 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. 19* 20* Parameters 21* ========== 22* 23* UPLO - CHARACTER*1. 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: 27* 28* UPLO = 'U' or 'u' The upper triangular part of A is 29* supplied in AP. 30* 31* UPLO = 'L' or 'l' The lower triangular part of A is 32* supplied in AP. 33* 34* Unchanged on exit. 35* 36* N - INTEGER. 37* On entry, N specifies the order of the matrix A. 38* N must be at least zero. 39* Unchanged on exit. 40* 41* ALPHA - DOUBLE PRECISION. 42* On entry, ALPHA specifies the scalar alpha. 43* Unchanged on exit. 44* 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. 57* Unchanged on exit. 58* 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 62* element vector x. 63* Unchanged on exit. 64* 65* INCX - INTEGER. 66* On entry, INCX specifies the increment for the elements of 67* X. INCX must not be zero. 68* Unchanged on exit. 69* 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. 73* Unchanged on exit. 74* 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 79* vector y. 80* 81* INCY - INTEGER. 82* On entry, INCY specifies the increment for the elements of 83* Y. INCY must not be zero. 84* Unchanged on exit. 85* 86* 87* Level 2 Blas routine. 88* 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. 94* 95* 96* .. Parameters .. 97 DOUBLE PRECISION ONE , ZERO 98 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 99* .. Local Scalars .. 100 DOUBLE PRECISION TEMP1, TEMP2 101 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY 102* .. External Functions .. 103 LOGICAL LSAME 104 EXTERNAL LSAME 105* .. External Subroutines .. 106 EXTERNAL XERBLA 107* .. 108* .. Executable Statements .. 109* 110* Test the input parameters. 111* 112 INFO = 0 113 IF ( .NOT.LSAME( UPLO, 'U' ).AND. 114 $ .NOT.LSAME( UPLO, 'L' ) )THEN 115 INFO = 1 116 ELSE IF( N.LT.0 )THEN 117 INFO = 2 118 ELSE IF( INCX.EQ.0 )THEN 119 INFO = 6 120 ELSE IF( INCY.EQ.0 )THEN 121 INFO = 9 122 END IF 123 IF( INFO.NE.0 )THEN 124 CALL XERBLA( 'DSPMV ', INFO ) 125 RETURN 126 END IF 127* 128* Quick return if possible. 129* 130 IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) 131 $ RETURN 132* 133* Set up the start points in X and Y. 134* 135 IF( INCX.GT.0 )THEN 136 KX = 1 137 ELSE 138 KX = 1 - ( N - 1 )*INCX 139 END IF 140 IF( INCY.GT.0 )THEN 141 KY = 1 142 ELSE 143 KY = 1 - ( N - 1 )*INCY 144 END IF 145* 146* Start the operations. In this version the elements of the array AP 147* are accessed sequentially with one pass through AP. 148* 149* First form y := beta*y. 150* 151 IF( BETA.NE.ONE )THEN 152 IF( INCY.EQ.1 )THEN 153 IF( BETA.EQ.ZERO )THEN 154 DO 10, I = 1, N 155 Y( I ) = ZERO 156 10 CONTINUE 157 ELSE 158 DO 20, I = 1, N 159 Y( I ) = BETA*Y( I ) 160 20 CONTINUE 161 END IF 162 ELSE 163 IY = KY 164 IF( BETA.EQ.ZERO )THEN 165 DO 30, I = 1, N 166 Y( IY ) = ZERO 167 IY = IY + INCY 168 30 CONTINUE 169 ELSE 170 DO 40, I = 1, N 171 Y( IY ) = BETA*Y( IY ) 172 IY = IY + INCY 173 40 CONTINUE 174 END IF 175 END IF 176 END IF 177 IF( ALPHA.EQ.ZERO ) 178 $ RETURN 179 KK = 1 180 IF( LSAME( UPLO, 'U' ) )THEN 181* 182* Form y when AP contains the upper triangle. 183* 184 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN 185 DO 60, J = 1, N 186 TEMP1 = ALPHA*X( J ) 187 TEMP2 = ZERO 188 K = KK 189 DO 50, I = 1, J - 1 190 Y( I ) = Y( I ) + TEMP1*AP( K ) 191 TEMP2 = TEMP2 + AP( K )*X( I ) 192 K = K + 1 193 50 CONTINUE 194 Y( J ) = Y( J ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 195 KK = KK + J 196 60 CONTINUE 197 ELSE 198 JX = KX 199 JY = KY 200 DO 80, J = 1, N 201 TEMP1 = ALPHA*X( JX ) 202 TEMP2 = ZERO 203 IX = KX 204 IY = KY 205 DO 70, K = KK, KK + J - 2 206 Y( IY ) = Y( IY ) + TEMP1*AP( K ) 207 TEMP2 = TEMP2 + AP( K )*X( IX ) 208 IX = IX + INCX 209 IY = IY + INCY 210 70 CONTINUE 211 Y( JY ) = Y( JY ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 212 JX = JX + INCX 213 JY = JY + INCY 214 KK = KK + J 215 80 CONTINUE 216 END IF 217 ELSE 218* 219* Form y when AP contains the lower triangle. 220* 221 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN 222 DO 100, J = 1, N 223 TEMP1 = ALPHA*X( J ) 224 TEMP2 = ZERO 225 Y( J ) = Y( J ) + TEMP1*AP( KK ) 226 K = KK + 1 227 DO 90, I = J + 1, N 228 Y( I ) = Y( I ) + TEMP1*AP( K ) 229 TEMP2 = TEMP2 + AP( K )*X( I ) 230 K = K + 1 231 90 CONTINUE 232 Y( J ) = Y( J ) + ALPHA*TEMP2 233 KK = KK + ( N - J + 1 ) 234 100 CONTINUE 235 ELSE 236 JX = KX 237 JY = KY 238 DO 120, J = 1, N 239 TEMP1 = ALPHA*X( JX ) 240 TEMP2 = ZERO 241 Y( JY ) = Y( JY ) + TEMP1*AP( KK ) 242 IX = JX 243 IY = JY 244 DO 110, K = KK + 1, KK + N - J 245 IX = IX + INCX 246 IY = IY + INCY 247 Y( IY ) = Y( IY ) + TEMP1*AP( K ) 248 TEMP2 = TEMP2 + AP( K )*X( IX ) 249 110 CONTINUE 250 Y( JY ) = Y( JY ) + ALPHA*TEMP2 251 JX = JX + INCX 252 JY = JY + INCY 253 KK = KK + ( N - J + 1 ) 254 120 CONTINUE 255 END IF 256 END IF 257* 258 RETURN 259* 260* End of DSPMV . 261* 262 END 263