1 /* chbmv.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 chbmv_(char *uplo, integer *n, integer *k, complex * 16 alpha, complex *a, integer *lda, complex *x, integer *incx, complex * 17 beta, complex *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, i__5; 21 real r__1; 22 complex q__1, q__2, q__3, q__4; 23 24 /* Builtin functions */ 25 void r_cnjg(complex *, complex *); 26 27 /* Local variables */ 28 integer i__, j, l, ix, iy, jx, jy, kx, ky, info; 29 complex temp1, temp2; 30 extern logical lsame_(char *, char *, ftnlen, ftnlen); 31 integer kplus1; 32 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen); 33 34 /* .. Scalar Arguments .. */ 35 /* .. */ 36 /* .. Array Arguments .. */ 37 /* .. */ 38 39 /* Purpose */ 40 /* ======= */ 41 42 /* CHBMV performs the matrix-vector operation */ 43 44 /* y := alpha*A*x + beta*y, */ 45 46 /* where alpha and beta are scalars, x and y are n element vectors and */ 47 /* A is an n by n hermitian band matrix, with k super-diagonals. */ 48 49 /* Arguments */ 50 /* ========== */ 51 52 /* UPLO - CHARACTER*1. */ 53 /* On entry, UPLO specifies whether the upper or lower */ 54 /* triangular part of the band matrix A is being supplied as */ 55 /* follows: */ 56 57 /* UPLO = 'U' or 'u' The upper triangular part of A is */ 58 /* being supplied. */ 59 60 /* UPLO = 'L' or 'l' The lower triangular part of A is */ 61 /* being supplied. */ 62 63 /* Unchanged on exit. */ 64 65 /* N - INTEGER. */ 66 /* On entry, N specifies the order of the matrix A. */ 67 /* N must be at least zero. */ 68 /* Unchanged on exit. */ 69 70 /* K - INTEGER. */ 71 /* On entry, K specifies the number of super-diagonals of the */ 72 /* matrix A. K must satisfy 0 .le. K. */ 73 /* Unchanged on exit. */ 74 75 /* ALPHA - COMPLEX . */ 76 /* On entry, ALPHA specifies the scalar alpha. */ 77 /* Unchanged on exit. */ 78 79 /* A - COMPLEX array of DIMENSION ( LDA, n ). */ 80 /* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */ 81 /* by n part of the array A must contain the upper triangular */ 82 /* band part of the hermitian matrix, supplied column by */ 83 /* column, with the leading diagonal of the matrix in row */ 84 /* ( k + 1 ) of the array, the first super-diagonal starting at */ 85 /* position 2 in row k, and so on. The top left k by k triangle */ 86 /* of the array A is not referenced. */ 87 /* The following program segment will transfer the upper */ 88 /* triangular part of a hermitian band matrix from conventional */ 89 /* full matrix storage to band storage: */ 90 91 /* DO 20, J = 1, N */ 92 /* M = K + 1 - J */ 93 /* DO 10, I = MAX( 1, J - K ), J */ 94 /* A( M + I, J ) = matrix( I, J ) */ 95 /* 10 CONTINUE */ 96 /* 20 CONTINUE */ 97 98 /* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */ 99 /* by n part of the array A must contain the lower triangular */ 100 /* band part of the hermitian matrix, supplied column by */ 101 /* column, with the leading diagonal of the matrix in row 1 of */ 102 /* the array, the first sub-diagonal starting at position 1 in */ 103 /* row 2, and so on. The bottom right k by k triangle of the */ 104 /* array A is not referenced. */ 105 /* The following program segment will transfer the lower */ 106 /* triangular part of a hermitian band matrix from conventional */ 107 /* full matrix storage to band storage: */ 108 109 /* DO 20, J = 1, N */ 110 /* M = 1 - J */ 111 /* DO 10, I = J, MIN( N, J + K ) */ 112 /* A( M + I, J ) = matrix( I, J ) */ 113 /* 10 CONTINUE */ 114 /* 20 CONTINUE */ 115 116 /* Note that the imaginary parts of the diagonal elements need */ 117 /* not be set and are assumed to be zero. */ 118 /* Unchanged on exit. */ 119 120 /* LDA - INTEGER. */ 121 /* On entry, LDA specifies the first dimension of A as declared */ 122 /* in the calling (sub) program. LDA must be at least */ 123 /* ( k + 1 ). */ 124 /* Unchanged on exit. */ 125 126 /* X - COMPLEX array of DIMENSION at least */ 127 /* ( 1 + ( n - 1 )*abs( INCX ) ). */ 128 /* Before entry, the incremented array X must contain the */ 129 /* vector x. */ 130 /* Unchanged on exit. */ 131 132 /* INCX - INTEGER. */ 133 /* On entry, INCX specifies the increment for the elements of */ 134 /* X. INCX must not be zero. */ 135 /* Unchanged on exit. */ 136 137 /* BETA - COMPLEX . */ 138 /* On entry, BETA specifies the scalar beta. */ 139 /* Unchanged on exit. */ 140 141 /* Y - COMPLEX array of DIMENSION at least */ 142 /* ( 1 + ( n - 1 )*abs( INCY ) ). */ 143 /* Before entry, the incremented array Y must contain the */ 144 /* vector y. On exit, Y is overwritten by the updated vector y. */ 145 146 /* INCY - INTEGER. */ 147 /* On entry, INCY specifies the increment for the elements of */ 148 /* Y. INCY must not be zero. */ 149 /* Unchanged on exit. */ 150 151 /* Further Details */ 152 /* =============== */ 153 154 /* Level 2 Blas routine. */ 155 156 /* -- Written on 22-October-1986. */ 157 /* Jack Dongarra, Argonne National Lab. */ 158 /* Jeremy Du Croz, Nag Central Office. */ 159 /* Sven Hammarling, Nag Central Office. */ 160 /* Richard Hanson, Sandia National Labs. */ 161 162 /* ===================================================================== */ 163 164 /* .. Parameters .. */ 165 /* .. */ 166 /* .. Local Scalars .. */ 167 /* .. */ 168 /* .. External Functions .. */ 169 /* .. */ 170 /* .. External Subroutines .. */ 171 /* .. */ 172 /* .. Intrinsic Functions .. */ 173 /* .. */ 174 175 /* Test the input parameters. */ 176 177 /* Parameter adjustments */ 178 a_dim1 = *lda; 179 a_offset = 1 + a_dim1; 180 a -= a_offset; 181 --x; 182 --y; 183 184 /* Function Body */ 185 info = 0; 186 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", ( 187 ftnlen)1, (ftnlen)1)) { 188 info = 1; 189 } else if (*n < 0) { 190 info = 2; 191 } else if (*k < 0) { 192 info = 3; 193 } else if (*lda < *k + 1) { 194 info = 6; 195 } else if (*incx == 0) { 196 info = 8; 197 } else if (*incy == 0) { 198 info = 11; 199 } 200 if (info != 0) { 201 xerbla_("CHBMV ", &info, (ftnlen)6); 202 return 0; 203 } 204 205 /* Quick return if possible. */ 206 207 if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f && (beta->r == 1.f && 208 beta->i == 0.f))) { 209 return 0; 210 } 211 212 /* Set up the start points in X and Y. */ 213 214 if (*incx > 0) { 215 kx = 1; 216 } else { 217 kx = 1 - (*n - 1) * *incx; 218 } 219 if (*incy > 0) { 220 ky = 1; 221 } else { 222 ky = 1 - (*n - 1) * *incy; 223 } 224 225 /* Start the operations. In this version the elements of the array A */ 226 /* are accessed sequentially with one pass through A. */ 227 228 /* First form y := beta*y. */ 229 230 if (beta->r != 1.f || beta->i != 0.f) { 231 if (*incy == 1) { 232 if (beta->r == 0.f && beta->i == 0.f) { 233 i__1 = *n; 234 for (i__ = 1; i__ <= i__1; ++i__) { 235 i__2 = i__; 236 y[i__2].r = 0.f, y[i__2].i = 0.f; 237 /* L10: */ 238 } 239 } else { 240 i__1 = *n; 241 for (i__ = 1; i__ <= i__1; ++i__) { 242 i__2 = i__; 243 i__3 = i__; 244 q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 245 q__1.i = beta->r * y[i__3].i + beta->i * y[i__3] 246 .r; 247 y[i__2].r = q__1.r, y[i__2].i = q__1.i; 248 /* L20: */ 249 } 250 } 251 } else { 252 iy = ky; 253 if (beta->r == 0.f && beta->i == 0.f) { 254 i__1 = *n; 255 for (i__ = 1; i__ <= i__1; ++i__) { 256 i__2 = iy; 257 y[i__2].r = 0.f, y[i__2].i = 0.f; 258 iy += *incy; 259 /* L30: */ 260 } 261 } else { 262 i__1 = *n; 263 for (i__ = 1; i__ <= i__1; ++i__) { 264 i__2 = iy; 265 i__3 = iy; 266 q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 267 q__1.i = beta->r * y[i__3].i + beta->i * y[i__3] 268 .r; 269 y[i__2].r = q__1.r, y[i__2].i = q__1.i; 270 iy += *incy; 271 /* L40: */ 272 } 273 } 274 } 275 } 276 if (alpha->r == 0.f && alpha->i == 0.f) { 277 return 0; 278 } 279 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { 280 281 /* Form y when upper triangle of A is stored. */ 282 283 kplus1 = *k + 1; 284 if (*incx == 1 && *incy == 1) { 285 i__1 = *n; 286 for (j = 1; j <= i__1; ++j) { 287 i__2 = j; 288 q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i = 289 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 290 temp1.r = q__1.r, temp1.i = q__1.i; 291 temp2.r = 0.f, temp2.i = 0.f; 292 l = kplus1 - j; 293 /* Computing MAX */ 294 i__2 = 1, i__3 = j - *k; 295 i__4 = j - 1; 296 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { 297 i__2 = i__; 298 i__3 = i__; 299 i__5 = l + i__ + j * a_dim1; 300 q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 301 q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5] 302 .r; 303 q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i; 304 y[i__2].r = q__1.r, y[i__2].i = q__1.i; 305 r_cnjg(&q__3, &a[l + i__ + j * a_dim1]); 306 i__2 = i__; 307 q__2.r = q__3.r * x[i__2].r - q__3.i * x[i__2].i, q__2.i = 308 q__3.r * x[i__2].i + q__3.i * x[i__2].r; 309 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 310 temp2.r = q__1.r, temp2.i = q__1.i; 311 /* L50: */ 312 } 313 i__4 = j; 314 i__2 = j; 315 i__3 = kplus1 + j * a_dim1; 316 r__1 = a[i__3].r; 317 q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i; 318 q__2.r = y[i__2].r + q__3.r, q__2.i = y[i__2].i + q__3.i; 319 q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 320 alpha->r * temp2.i + alpha->i * temp2.r; 321 q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i; 322 y[i__4].r = q__1.r, y[i__4].i = q__1.i; 323 /* L60: */ 324 } 325 } else { 326 jx = kx; 327 jy = ky; 328 i__1 = *n; 329 for (j = 1; j <= i__1; ++j) { 330 i__4 = jx; 331 q__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, q__1.i = 332 alpha->r * x[i__4].i + alpha->i * x[i__4].r; 333 temp1.r = q__1.r, temp1.i = q__1.i; 334 temp2.r = 0.f, temp2.i = 0.f; 335 ix = kx; 336 iy = ky; 337 l = kplus1 - j; 338 /* Computing MAX */ 339 i__4 = 1, i__2 = j - *k; 340 i__3 = j - 1; 341 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) { 342 i__4 = iy; 343 i__2 = iy; 344 i__5 = l + i__ + j * a_dim1; 345 q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 346 q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5] 347 .r; 348 q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i; 349 y[i__4].r = q__1.r, y[i__4].i = q__1.i; 350 r_cnjg(&q__3, &a[l + i__ + j * a_dim1]); 351 i__4 = ix; 352 q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i = 353 q__3.r * x[i__4].i + q__3.i * x[i__4].r; 354 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 355 temp2.r = q__1.r, temp2.i = q__1.i; 356 ix += *incx; 357 iy += *incy; 358 /* L70: */ 359 } 360 i__3 = jy; 361 i__4 = jy; 362 i__2 = kplus1 + j * a_dim1; 363 r__1 = a[i__2].r; 364 q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i; 365 q__2.r = y[i__4].r + q__3.r, q__2.i = y[i__4].i + q__3.i; 366 q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 367 alpha->r * temp2.i + alpha->i * temp2.r; 368 q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i; 369 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 370 jx += *incx; 371 jy += *incy; 372 if (j > *k) { 373 kx += *incx; 374 ky += *incy; 375 } 376 /* L80: */ 377 } 378 } 379 } else { 380 381 /* Form y when lower triangle of A is stored. */ 382 383 if (*incx == 1 && *incy == 1) { 384 i__1 = *n; 385 for (j = 1; j <= i__1; ++j) { 386 i__3 = j; 387 q__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, q__1.i = 388 alpha->r * x[i__3].i + alpha->i * x[i__3].r; 389 temp1.r = q__1.r, temp1.i = q__1.i; 390 temp2.r = 0.f, temp2.i = 0.f; 391 i__3 = j; 392 i__4 = j; 393 i__2 = j * a_dim1 + 1; 394 r__1 = a[i__2].r; 395 q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i; 396 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 397 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 398 l = 1 - j; 399 /* Computing MIN */ 400 i__4 = *n, i__2 = j + *k; 401 i__3 = min(i__4,i__2); 402 for (i__ = j + 1; i__ <= i__3; ++i__) { 403 i__4 = i__; 404 i__2 = i__; 405 i__5 = l + i__ + j * a_dim1; 406 q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 407 q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5] 408 .r; 409 q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i; 410 y[i__4].r = q__1.r, y[i__4].i = q__1.i; 411 r_cnjg(&q__3, &a[l + i__ + j * a_dim1]); 412 i__4 = i__; 413 q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i = 414 q__3.r * x[i__4].i + q__3.i * x[i__4].r; 415 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 416 temp2.r = q__1.r, temp2.i = q__1.i; 417 /* L90: */ 418 } 419 i__3 = j; 420 i__4 = j; 421 q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 422 alpha->r * temp2.i + alpha->i * temp2.r; 423 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 424 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 425 /* L100: */ 426 } 427 } else { 428 jx = kx; 429 jy = ky; 430 i__1 = *n; 431 for (j = 1; j <= i__1; ++j) { 432 i__3 = jx; 433 q__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, q__1.i = 434 alpha->r * x[i__3].i + alpha->i * x[i__3].r; 435 temp1.r = q__1.r, temp1.i = q__1.i; 436 temp2.r = 0.f, temp2.i = 0.f; 437 i__3 = jy; 438 i__4 = jy; 439 i__2 = j * a_dim1 + 1; 440 r__1 = a[i__2].r; 441 q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i; 442 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 443 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 444 l = 1 - j; 445 ix = jx; 446 iy = jy; 447 /* Computing MIN */ 448 i__4 = *n, i__2 = j + *k; 449 i__3 = min(i__4,i__2); 450 for (i__ = j + 1; i__ <= i__3; ++i__) { 451 ix += *incx; 452 iy += *incy; 453 i__4 = iy; 454 i__2 = iy; 455 i__5 = l + i__ + j * a_dim1; 456 q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 457 q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5] 458 .r; 459 q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i; 460 y[i__4].r = q__1.r, y[i__4].i = q__1.i; 461 r_cnjg(&q__3, &a[l + i__ + j * a_dim1]); 462 i__4 = ix; 463 q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i = 464 q__3.r * x[i__4].i + q__3.i * x[i__4].r; 465 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 466 temp2.r = q__1.r, temp2.i = q__1.i; 467 /* L110: */ 468 } 469 i__3 = jy; 470 i__4 = jy; 471 q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 472 alpha->r * temp2.i + alpha->i * temp2.r; 473 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 474 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 475 jx += *incx; 476 jy += *incy; 477 /* L120: */ 478 } 479 } 480 } 481 482 return 0; 483 484 /* End of CHBMV . */ 485 486 } /* chbmv_ */ 487 488
