1 /* dsbmv.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 dsbmv_(char *uplo, integer *n, integer *k, doublereal * 16 alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, 17 doublereal *beta, doublereal *y, integer *incy, ftnlen uplo_len) 18 { 19 /* System generated locals */ 20 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 21 22 /* Local variables */ 23 integer i__, j, l, ix, iy, jx, jy, kx, ky, info; 24 doublereal temp1, temp2; 25 extern logical lsame_(char *, char *, ftnlen, ftnlen); 26 integer kplus1; 27 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen); 28 29 /* .. Scalar Arguments .. */ 30 /* .. */ 31 /* .. Array Arguments .. */ 32 /* .. */ 33 34 /* Purpose */ 35 /* ======= */ 36 37 /* DSBMV performs the matrix-vector operation */ 38 39 /* y := alpha*A*x + beta*y, */ 40 41 /* where alpha and beta are scalars, x and y are n element vectors and */ 42 /* A is an n by n symmetric band matrix, with k super-diagonals. */ 43 44 /* Arguments */ 45 /* ========== */ 46 47 /* UPLO - CHARACTER*1. */ 48 /* On entry, UPLO specifies whether the upper or lower */ 49 /* triangular part of the band matrix A is being supplied as */ 50 /* follows: */ 51 52 /* UPLO = 'U' or 'u' The upper triangular part of A is */ 53 /* being supplied. */ 54 55 /* UPLO = 'L' or 'l' The lower triangular part of A is */ 56 /* being supplied. */ 57 58 /* Unchanged on exit. */ 59 60 /* N - INTEGER. */ 61 /* On entry, N specifies the order of the matrix A. */ 62 /* N must be at least zero. */ 63 /* Unchanged on exit. */ 64 65 /* K - INTEGER. */ 66 /* On entry, K specifies the number of super-diagonals of the */ 67 /* matrix A. K must satisfy 0 .le. K. */ 68 /* Unchanged on exit. */ 69 70 /* ALPHA - DOUBLE PRECISION. */ 71 /* On entry, ALPHA specifies the scalar alpha. */ 72 /* Unchanged on exit. */ 73 74 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */ 75 /* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */ 76 /* by n part of the array A must contain the upper triangular */ 77 /* band part of the symmetric matrix, supplied column by */ 78 /* column, with the leading diagonal of the matrix in row */ 79 /* ( k + 1 ) of the array, the first super-diagonal starting at */ 80 /* position 2 in row k, and so on. The top left k by k triangle */ 81 /* of the array A is not referenced. */ 82 /* The following program segment will transfer the upper */ 83 /* triangular part of a symmetric band matrix from conventional */ 84 /* full matrix storage to band storage: */ 85 86 /* DO 20, J = 1, N */ 87 /* M = K + 1 - J */ 88 /* DO 10, I = MAX( 1, J - K ), J */ 89 /* A( M + I, J ) = matrix( I, J ) */ 90 /* 10 CONTINUE */ 91 /* 20 CONTINUE */ 92 93 /* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */ 94 /* by n part of the array A must contain the lower triangular */ 95 /* band part of the symmetric matrix, supplied column by */ 96 /* column, with the leading diagonal of the matrix in row 1 of */ 97 /* the array, the first sub-diagonal starting at position 1 in */ 98 /* row 2, and so on. The bottom right k by k triangle of the */ 99 /* array A is not referenced. */ 100 /* The following program segment will transfer the lower */ 101 /* triangular part of a symmetric band matrix from conventional */ 102 /* full matrix storage to band storage: */ 103 104 /* DO 20, J = 1, N */ 105 /* M = 1 - J */ 106 /* DO 10, I = J, MIN( N, J + K ) */ 107 /* A( M + I, J ) = matrix( I, J ) */ 108 /* 10 CONTINUE */ 109 /* 20 CONTINUE */ 110 111 /* Unchanged on exit. */ 112 113 /* LDA - INTEGER. */ 114 /* On entry, LDA specifies the first dimension of A as declared */ 115 /* in the calling (sub) program. LDA must be at least */ 116 /* ( k + 1 ). */ 117 /* Unchanged on exit. */ 118 119 /* X - DOUBLE PRECISION array of DIMENSION at least */ 120 /* ( 1 + ( n - 1 )*abs( INCX ) ). */ 121 /* Before entry, the incremented array X must contain the */ 122 /* vector x. */ 123 /* Unchanged on exit. */ 124 125 /* INCX - INTEGER. */ 126 /* On entry, INCX specifies the increment for the elements of */ 127 /* X. INCX must not be zero. */ 128 /* Unchanged on exit. */ 129 130 /* BETA - DOUBLE PRECISION. */ 131 /* On entry, BETA specifies the scalar beta. */ 132 /* Unchanged on exit. */ 133 134 /* Y - DOUBLE PRECISION array of DIMENSION at least */ 135 /* ( 1 + ( n - 1 )*abs( INCY ) ). */ 136 /* Before entry, the incremented array Y must contain the */ 137 /* vector y. On exit, Y is overwritten by the updated vector y. */ 138 139 /* INCY - INTEGER. */ 140 /* On entry, INCY specifies the increment for the elements of */ 141 /* Y. INCY must not be zero. */ 142 /* Unchanged on exit. */ 143 144 145 /* Level 2 Blas routine. */ 146 147 /* -- Written on 22-October-1986. */ 148 /* Jack Dongarra, Argonne National Lab. */ 149 /* Jeremy Du Croz, Nag Central Office. */ 150 /* Sven Hammarling, Nag Central Office. */ 151 /* Richard Hanson, Sandia National Labs. */ 152 153 /* ===================================================================== */ 154 155 /* .. Parameters .. */ 156 /* .. */ 157 /* .. Local Scalars .. */ 158 /* .. */ 159 /* .. External Functions .. */ 160 /* .. */ 161 /* .. External Subroutines .. */ 162 /* .. */ 163 /* .. Intrinsic Functions .. */ 164 /* .. */ 165 166 /* Test the input parameters. */ 167 168 /* Parameter adjustments */ 169 a_dim1 = *lda; 170 a_offset = 1 + a_dim1; 171 a -= a_offset; 172 --x; 173 --y; 174 175 /* Function Body */ 176 info = 0; 177 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", ( 178 ftnlen)1, (ftnlen)1)) { 179 info = 1; 180 } else if (*n < 0) { 181 info = 2; 182 } else if (*k < 0) { 183 info = 3; 184 } else if (*lda < *k + 1) { 185 info = 6; 186 } else if (*incx == 0) { 187 info = 8; 188 } else if (*incy == 0) { 189 info = 11; 190 } 191 if (info != 0) { 192 xerbla_("DSBMV ", &info, (ftnlen)6); 193 return 0; 194 } 195 196 /* Quick return if possible. */ 197 198 if (*n == 0 || (*alpha == 0. && *beta == 1.)) { 199 return 0; 200 } 201 202 /* Set up the start points in X and Y. */ 203 204 if (*incx > 0) { 205 kx = 1; 206 } else { 207 kx = 1 - (*n - 1) * *incx; 208 } 209 if (*incy > 0) { 210 ky = 1; 211 } else { 212 ky = 1 - (*n - 1) * *incy; 213 } 214 215 /* Start the operations. In this version the elements of the array A */ 216 /* are accessed sequentially with one pass through A. */ 217 218 /* First form y := beta*y. */ 219 220 if (*beta != 1.) { 221 if (*incy == 1) { 222 if (*beta == 0.) { 223 i__1 = *n; 224 for (i__ = 1; i__ <= i__1; ++i__) { 225 y[i__] = 0.; 226 /* L10: */ 227 } 228 } else { 229 i__1 = *n; 230 for (i__ = 1; i__ <= i__1; ++i__) { 231 y[i__] = *beta * y[i__]; 232 /* L20: */ 233 } 234 } 235 } else { 236 iy = ky; 237 if (*beta == 0.) { 238 i__1 = *n; 239 for (i__ = 1; i__ <= i__1; ++i__) { 240 y[iy] = 0.; 241 iy += *incy; 242 /* L30: */ 243 } 244 } else { 245 i__1 = *n; 246 for (i__ = 1; i__ <= i__1; ++i__) { 247 y[iy] = *beta * y[iy]; 248 iy += *incy; 249 /* L40: */ 250 } 251 } 252 } 253 } 254 if (*alpha == 0.) { 255 return 0; 256 } 257 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { 258 259 /* Form y when upper triangle of A is stored. */ 260 261 kplus1 = *k + 1; 262 if (*incx == 1 && *incy == 1) { 263 i__1 = *n; 264 for (j = 1; j <= i__1; ++j) { 265 temp1 = *alpha * x[j]; 266 temp2 = 0.; 267 l = kplus1 - j; 268 /* Computing MAX */ 269 i__2 = 1, i__3 = j - *k; 270 i__4 = j - 1; 271 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { 272 y[i__] += temp1 * a[l + i__ + j * a_dim1]; 273 temp2 += a[l + i__ + j * a_dim1] * x[i__]; 274 /* L50: */ 275 } 276 y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2; 277 /* L60: */ 278 } 279 } else { 280 jx = kx; 281 jy = ky; 282 i__1 = *n; 283 for (j = 1; j <= i__1; ++j) { 284 temp1 = *alpha * x[jx]; 285 temp2 = 0.; 286 ix = kx; 287 iy = ky; 288 l = kplus1 - j; 289 /* Computing MAX */ 290 i__4 = 1, i__2 = j - *k; 291 i__3 = j - 1; 292 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) { 293 y[iy] += temp1 * a[l + i__ + j * a_dim1]; 294 temp2 += a[l + i__ + j * a_dim1] * x[ix]; 295 ix += *incx; 296 iy += *incy; 297 /* L70: */ 298 } 299 y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha * 300 temp2; 301 jx += *incx; 302 jy += *incy; 303 if (j > *k) { 304 kx += *incx; 305 ky += *incy; 306 } 307 /* L80: */ 308 } 309 } 310 } else { 311 312 /* Form y when lower triangle of A is stored. */ 313 314 if (*incx == 1 && *incy == 1) { 315 i__1 = *n; 316 for (j = 1; j <= i__1; ++j) { 317 temp1 = *alpha * x[j]; 318 temp2 = 0.; 319 y[j] += temp1 * a[j * a_dim1 + 1]; 320 l = 1 - j; 321 /* Computing MIN */ 322 i__4 = *n, i__2 = j + *k; 323 i__3 = min(i__4,i__2); 324 for (i__ = j + 1; i__ <= i__3; ++i__) { 325 y[i__] += temp1 * a[l + i__ + j * a_dim1]; 326 temp2 += a[l + i__ + j * a_dim1] * x[i__]; 327 /* L90: */ 328 } 329 y[j] += *alpha * temp2; 330 /* L100: */ 331 } 332 } else { 333 jx = kx; 334 jy = ky; 335 i__1 = *n; 336 for (j = 1; j <= i__1; ++j) { 337 temp1 = *alpha * x[jx]; 338 temp2 = 0.; 339 y[jy] += temp1 * a[j * a_dim1 + 1]; 340 l = 1 - j; 341 ix = jx; 342 iy = jy; 343 /* Computing MIN */ 344 i__4 = *n, i__2 = j + *k; 345 i__3 = min(i__4,i__2); 346 for (i__ = j + 1; i__ <= i__3; ++i__) { 347 ix += *incx; 348 iy += *incy; 349 y[iy] += temp1 * a[l + i__ + j * a_dim1]; 350 temp2 += a[l + i__ + j * a_dim1] * x[ix]; 351 /* L110: */ 352 } 353 y[jy] += *alpha * temp2; 354 jx += *incx; 355 jy += *incy; 356 /* L120: */ 357 } 358 } 359 } 360 361 return 0; 362 363 /* End of DSBMV . */ 364 365 } /* dsbmv_ */ 366 367
