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