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