1 /* dspmv.f -- translated by f2c (version 20100827). 2 You must link the resulting object file with libf2c: 3 on Microsoft Windows system, link with libf2c.lib; 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 5 or, if you install libf2c.a in a standard place, with -lf2c -lm 6 -- in that order, at the end of the command line, as in 7 cc *.o -lf2c -lm 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 9 10 http://www.netlib.org/f2c/libf2c.zip 11 */ 12 13 #include "datatypes.h" 14 15 /* Subroutine */ int dspmv_(char *uplo, integer *n, doublereal *alpha, 16 doublereal *ap, doublereal *x, integer *incx, doublereal *beta, 17 doublereal *y, integer *incy, ftnlen uplo_len) 18 { 19 /* System generated locals */ 20 integer i__1, i__2; 21 22 /* Local variables */ 23 integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info; 24 doublereal temp1, temp2; 25 extern logical lsame_(char *, char *, ftnlen, ftnlen); 26 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen); 27 28 /* .. Scalar Arguments .. */ 29 /* .. */ 30 /* .. Array Arguments .. */ 31 /* .. */ 32 33 /* Purpose */ 34 /* ======= */ 35 36 /* DSPMV performs the matrix-vector operation */ 37 38 /* y := alpha*A*x + beta*y, */ 39 40 /* where alpha and beta are scalars, x and y are n element vectors and */ 41 /* A is an n by n symmetric matrix, supplied in packed form. */ 42 43 /* Arguments */ 44 /* ========== */ 45 46 /* UPLO - CHARACTER*1. */ 47 /* On entry, UPLO specifies whether the upper or lower */ 48 /* triangular part of the matrix A is supplied in the packed */ 49 /* array AP as follows: */ 50 51 /* UPLO = 'U' or 'u' The upper triangular part of A is */ 52 /* supplied in AP. */ 53 54 /* UPLO = 'L' or 'l' The lower triangular part of A is */ 55 /* supplied in AP. */ 56 57 /* Unchanged on exit. */ 58 59 /* N - INTEGER. */ 60 /* On entry, N specifies the order of the matrix A. */ 61 /* N must be at least zero. */ 62 /* Unchanged on exit. */ 63 64 /* ALPHA - DOUBLE PRECISION. */ 65 /* On entry, ALPHA specifies the scalar alpha. */ 66 /* Unchanged on exit. */ 67 68 /* AP - DOUBLE PRECISION array of DIMENSION at least */ 69 /* ( ( n*( n + 1 ) )/2 ). */ 70 /* Before entry with UPLO = 'U' or 'u', the array AP must */ 71 /* contain the upper triangular part of the symmetric matrix */ 72 /* packed sequentially, column by column, so that AP( 1 ) */ 73 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */ 74 /* and a( 2, 2 ) respectively, and so on. */ 75 /* Before entry with UPLO = 'L' or 'l', the array AP must */ 76 /* contain the lower triangular part of the symmetric matrix */ 77 /* packed sequentially, column by column, so that AP( 1 ) */ 78 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */ 79 /* and a( 3, 1 ) respectively, and so on. */ 80 /* Unchanged on exit. */ 81 82 /* X - DOUBLE PRECISION array of dimension at least */ 83 /* ( 1 + ( n - 1 )*abs( INCX ) ). */ 84 /* Before entry, the incremented array X must contain the n */ 85 /* element vector x. */ 86 /* Unchanged on exit. */ 87 88 /* INCX - INTEGER. */ 89 /* On entry, INCX specifies the increment for the elements of */ 90 /* X. INCX must not be zero. */ 91 /* Unchanged on exit. */ 92 93 /* BETA - DOUBLE PRECISION. */ 94 /* On entry, BETA specifies the scalar beta. When BETA is */ 95 /* supplied as zero then Y need not be set on input. */ 96 /* Unchanged on exit. */ 97 98 /* Y - DOUBLE PRECISION array of dimension at least */ 99 /* ( 1 + ( n - 1 )*abs( INCY ) ). */ 100 /* Before entry, the incremented array Y must contain the n */ 101 /* element vector y. On exit, Y is overwritten by the updated */ 102 /* vector y. */ 103 104 /* INCY - INTEGER. */ 105 /* On entry, INCY specifies the increment for the elements of */ 106 /* Y. INCY must not be zero. */ 107 /* Unchanged on exit. */ 108 109 /* Further Details */ 110 /* =============== */ 111 112 /* Level 2 Blas routine. */ 113 114 /* -- Written on 22-October-1986. */ 115 /* Jack Dongarra, Argonne National Lab. */ 116 /* Jeremy Du Croz, Nag Central Office. */ 117 /* Sven Hammarling, Nag Central Office. */ 118 /* Richard Hanson, Sandia National Labs. */ 119 120 /* ===================================================================== */ 121 122 /* .. Parameters .. */ 123 /* .. */ 124 /* .. Local Scalars .. */ 125 /* .. */ 126 /* .. External Functions .. */ 127 /* .. */ 128 /* .. External Subroutines .. */ 129 /* .. */ 130 131 /* Test the input parameters. */ 132 133 /* Parameter adjustments */ 134 --y; 135 --x; 136 --ap; 137 138 /* Function Body */ 139 info = 0; 140 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", ( 141 ftnlen)1, (ftnlen)1)) { 142 info = 1; 143 } else if (*n < 0) { 144 info = 2; 145 } else if (*incx == 0) { 146 info = 6; 147 } else if (*incy == 0) { 148 info = 9; 149 } 150 if (info != 0) { 151 xerbla_("DSPMV ", &info, (ftnlen)6); 152 return 0; 153 } 154 155 /* Quick return if possible. */ 156 157 if (*n == 0 || (*alpha == 0. && *beta == 1.)) { 158 return 0; 159 } 160 161 /* Set up the start points in X and Y. */ 162 163 if (*incx > 0) { 164 kx = 1; 165 } else { 166 kx = 1 - (*n - 1) * *incx; 167 } 168 if (*incy > 0) { 169 ky = 1; 170 } else { 171 ky = 1 - (*n - 1) * *incy; 172 } 173 174 /* Start the operations. In this version the elements of the array AP */ 175 /* are accessed sequentially with one pass through AP. */ 176 177 /* First form y := beta*y. */ 178 179 if (*beta != 1.) { 180 if (*incy == 1) { 181 if (*beta == 0.) { 182 i__1 = *n; 183 for (i__ = 1; i__ <= i__1; ++i__) { 184 y[i__] = 0.; 185 /* L10: */ 186 } 187 } else { 188 i__1 = *n; 189 for (i__ = 1; i__ <= i__1; ++i__) { 190 y[i__] = *beta * y[i__]; 191 /* L20: */ 192 } 193 } 194 } else { 195 iy = ky; 196 if (*beta == 0.) { 197 i__1 = *n; 198 for (i__ = 1; i__ <= i__1; ++i__) { 199 y[iy] = 0.; 200 iy += *incy; 201 /* L30: */ 202 } 203 } else { 204 i__1 = *n; 205 for (i__ = 1; i__ <= i__1; ++i__) { 206 y[iy] = *beta * y[iy]; 207 iy += *incy; 208 /* L40: */ 209 } 210 } 211 } 212 } 213 if (*alpha == 0.) { 214 return 0; 215 } 216 kk = 1; 217 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { 218 219 /* Form y when AP contains the upper triangle. */ 220 221 if (*incx == 1 && *incy == 1) { 222 i__1 = *n; 223 for (j = 1; j <= i__1; ++j) { 224 temp1 = *alpha * x[j]; 225 temp2 = 0.; 226 k = kk; 227 i__2 = j - 1; 228 for (i__ = 1; i__ <= i__2; ++i__) { 229 y[i__] += temp1 * ap[k]; 230 temp2 += ap[k] * x[i__]; 231 ++k; 232 /* L50: */ 233 } 234 y[j] = y[j] + temp1 * ap[kk + j - 1] + *alpha * temp2; 235 kk += j; 236 /* L60: */ 237 } 238 } else { 239 jx = kx; 240 jy = ky; 241 i__1 = *n; 242 for (j = 1; j <= i__1; ++j) { 243 temp1 = *alpha * x[jx]; 244 temp2 = 0.; 245 ix = kx; 246 iy = ky; 247 i__2 = kk + j - 2; 248 for (k = kk; k <= i__2; ++k) { 249 y[iy] += temp1 * ap[k]; 250 temp2 += ap[k] * x[ix]; 251 ix += *incx; 252 iy += *incy; 253 /* L70: */ 254 } 255 y[jy] = y[jy] + temp1 * ap[kk + j - 1] + *alpha * temp2; 256 jx += *incx; 257 jy += *incy; 258 kk += j; 259 /* L80: */ 260 } 261 } 262 } else { 263 264 /* Form y when AP contains the lower triangle. */ 265 266 if (*incx == 1 && *incy == 1) { 267 i__1 = *n; 268 for (j = 1; j <= i__1; ++j) { 269 temp1 = *alpha * x[j]; 270 temp2 = 0.; 271 y[j] += temp1 * ap[kk]; 272 k = kk + 1; 273 i__2 = *n; 274 for (i__ = j + 1; i__ <= i__2; ++i__) { 275 y[i__] += temp1 * ap[k]; 276 temp2 += ap[k] * x[i__]; 277 ++k; 278 /* L90: */ 279 } 280 y[j] += *alpha * temp2; 281 kk += *n - j + 1; 282 /* L100: */ 283 } 284 } else { 285 jx = kx; 286 jy = ky; 287 i__1 = *n; 288 for (j = 1; j <= i__1; ++j) { 289 temp1 = *alpha * x[jx]; 290 temp2 = 0.; 291 y[jy] += temp1 * ap[kk]; 292 ix = jx; 293 iy = jy; 294 i__2 = kk + *n - j; 295 for (k = kk + 1; k <= i__2; ++k) { 296 ix += *incx; 297 iy += *incy; 298 y[iy] += temp1 * ap[k]; 299 temp2 += ap[k] * x[ix]; 300 /* L110: */ 301 } 302 y[jy] += *alpha * temp2; 303 jx += *incx; 304 jy += *incy; 305 kk += *n - j + 1; 306 /* L120: */ 307 } 308 } 309 } 310 311 return 0; 312 313 /* End of DSPMV . */ 314 315 } /* dspmv_ */ 316 317
