1 /* ssbmv.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 ssbmv_(char *uplo, integer *n, integer *k, real *alpha, 16 real *a, integer *lda, real *x, integer *incx, real *beta, real *y, 17 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 real 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 /* SSBMV 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 - REAL . */ 71 /* On entry, ALPHA specifies the scalar alpha. */ 72 /* Unchanged on exit. */ 73 74 /* A - REAL 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 - REAL 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 - REAL . */ 131 /* On entry, BETA specifies the scalar beta. */ 132 /* Unchanged on exit. */ 133 134 /* Y - REAL 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 /* Further Details */ 145 /* =============== */ 146 147 /* Level 2 Blas routine. */ 148 149 /* -- Written on 22-October-1986. */ 150 /* Jack Dongarra, Argonne National Lab. */ 151 /* Jeremy Du Croz, Nag Central Office. */ 152 /* Sven Hammarling, Nag Central Office. */ 153 /* Richard Hanson, Sandia National Labs. */ 154 155 /* ===================================================================== */ 156 157 /* .. Parameters .. */ 158 /* .. */ 159 /* .. Local Scalars .. */ 160 /* .. */ 161 /* .. External Functions .. */ 162 /* .. */ 163 /* .. External Subroutines .. */ 164 /* .. */ 165 /* .. Intrinsic Functions .. */ 166 /* .. */ 167 168 /* Test the input parameters. */ 169 170 /* Parameter adjustments */ 171 a_dim1 = *lda; 172 a_offset = 1 + a_dim1; 173 a -= a_offset; 174 --x; 175 --y; 176 177 /* Function Body */ 178 info = 0; 179 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", ( 180 ftnlen)1, (ftnlen)1)) { 181 info = 1; 182 } else if (*n < 0) { 183 info = 2; 184 } else if (*k < 0) { 185 info = 3; 186 } else if (*lda < *k + 1) { 187 info = 6; 188 } else if (*incx == 0) { 189 info = 8; 190 } else if (*incy == 0) { 191 info = 11; 192 } 193 if (info != 0) { 194 xerbla_("SSBMV ", &info, (ftnlen)6); 195 return 0; 196 } 197 198 /* Quick return if possible. */ 199 200 if (*n == 0 || (*alpha == 0.f && *beta == 1.f)) { 201 return 0; 202 } 203 204 /* Set up the start points in X and Y. */ 205 206 if (*incx > 0) { 207 kx = 1; 208 } else { 209 kx = 1 - (*n - 1) * *incx; 210 } 211 if (*incy > 0) { 212 ky = 1; 213 } else { 214 ky = 1 - (*n - 1) * *incy; 215 } 216 217 /* Start the operations. In this version the elements of the array A */ 218 /* are accessed sequentially with one pass through A. */ 219 220 /* First form y := beta*y. */ 221 222 if (*beta != 1.f) { 223 if (*incy == 1) { 224 if (*beta == 0.f) { 225 i__1 = *n; 226 for (i__ = 1; i__ <= i__1; ++i__) { 227 y[i__] = 0.f; 228 /* L10: */ 229 } 230 } else { 231 i__1 = *n; 232 for (i__ = 1; i__ <= i__1; ++i__) { 233 y[i__] = *beta * y[i__]; 234 /* L20: */ 235 } 236 } 237 } else { 238 iy = ky; 239 if (*beta == 0.f) { 240 i__1 = *n; 241 for (i__ = 1; i__ <= i__1; ++i__) { 242 y[iy] = 0.f; 243 iy += *incy; 244 /* L30: */ 245 } 246 } else { 247 i__1 = *n; 248 for (i__ = 1; i__ <= i__1; ++i__) { 249 y[iy] = *beta * y[iy]; 250 iy += *incy; 251 /* L40: */ 252 } 253 } 254 } 255 } 256 if (*alpha == 0.f) { 257 return 0; 258 } 259 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { 260 261 /* Form y when upper triangle of A is stored. */ 262 263 kplus1 = *k + 1; 264 if (*incx == 1 && *incy == 1) { 265 i__1 = *n; 266 for (j = 1; j <= i__1; ++j) { 267 temp1 = *alpha * x[j]; 268 temp2 = 0.f; 269 l = kplus1 - j; 270 /* Computing MAX */ 271 i__2 = 1, i__3 = j - *k; 272 i__4 = j - 1; 273 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { 274 y[i__] += temp1 * a[l + i__ + j * a_dim1]; 275 temp2 += a[l + i__ + j * a_dim1] * x[i__]; 276 /* L50: */ 277 } 278 y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2; 279 /* L60: */ 280 } 281 } else { 282 jx = kx; 283 jy = ky; 284 i__1 = *n; 285 for (j = 1; j <= i__1; ++j) { 286 temp1 = *alpha * x[jx]; 287 temp2 = 0.f; 288 ix = kx; 289 iy = ky; 290 l = kplus1 - j; 291 /* Computing MAX */ 292 i__4 = 1, i__2 = j - *k; 293 i__3 = j - 1; 294 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) { 295 y[iy] += temp1 * a[l + i__ + j * a_dim1]; 296 temp2 += a[l + i__ + j * a_dim1] * x[ix]; 297 ix += *incx; 298 iy += *incy; 299 /* L70: */ 300 } 301 y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha * 302 temp2; 303 jx += *incx; 304 jy += *incy; 305 if (j > *k) { 306 kx += *incx; 307 ky += *incy; 308 } 309 /* L80: */ 310 } 311 } 312 } else { 313 314 /* Form y when lower triangle of A is stored. */ 315 316 if (*incx == 1 && *incy == 1) { 317 i__1 = *n; 318 for (j = 1; j <= i__1; ++j) { 319 temp1 = *alpha * x[j]; 320 temp2 = 0.f; 321 y[j] += temp1 * a[j * a_dim1 + 1]; 322 l = 1 - j; 323 /* Computing MIN */ 324 i__4 = *n, i__2 = j + *k; 325 i__3 = min(i__4,i__2); 326 for (i__ = j + 1; i__ <= i__3; ++i__) { 327 y[i__] += temp1 * a[l + i__ + j * a_dim1]; 328 temp2 += a[l + i__ + j * a_dim1] * x[i__]; 329 /* L90: */ 330 } 331 y[j] += *alpha * temp2; 332 /* L100: */ 333 } 334 } else { 335 jx = kx; 336 jy = ky; 337 i__1 = *n; 338 for (j = 1; j <= i__1; ++j) { 339 temp1 = *alpha * x[jx]; 340 temp2 = 0.f; 341 y[jy] += temp1 * a[j * a_dim1 + 1]; 342 l = 1 - j; 343 ix = jx; 344 iy = jy; 345 /* Computing MIN */ 346 i__4 = *n, i__2 = j + *k; 347 i__3 = min(i__4,i__2); 348 for (i__ = j + 1; i__ <= i__3; ++i__) { 349 ix += *incx; 350 iy += *incy; 351 y[iy] += temp1 * a[l + i__ + j * a_dim1]; 352 temp2 += a[l + i__ + j * a_dim1] * x[ix]; 353 /* L110: */ 354 } 355 y[jy] += *alpha * temp2; 356 jx += *incx; 357 jy += *incy; 358 /* L120: */ 359 } 360 } 361 } 362 363 return 0; 364 365 /* End of SSBMV . */ 366 367 } /* ssbmv_ */ 368 369
