123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488 |
- #include "datatypes.h"
- int zhbmv_(char *uplo, integer *n, integer *k, doublecomplex
- *alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *
- incx, doublecomplex *beta, doublecomplex *y, integer *incy, ftnlen
- uplo_len)
- {
-
- integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
- doublereal d__1;
- doublecomplex z__1, z__2, z__3, z__4;
-
- void d_cnjg(doublecomplex *, doublecomplex *);
-
- integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
- doublecomplex temp1, temp2;
- extern logical lsame_(char *, char *, ftnlen, ftnlen);
- integer kplus1;
- extern int xerbla_(char *, integer *, ftnlen);
-
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- --x;
- --y;
-
- info = 0;
- if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
- ftnlen)1, (ftnlen)1)) {
- info = 1;
- } else if (*n < 0) {
- info = 2;
- } else if (*k < 0) {
- info = 3;
- } else if (*lda < *k + 1) {
- info = 6;
- } else if (*incx == 0) {
- info = 8;
- } else if (*incy == 0) {
- info = 11;
- }
- if (info != 0) {
- xerbla_("ZHBMV ", &info, (ftnlen)6);
- return 0;
- }
- if (*n == 0 || (alpha->r == 0. && alpha->i == 0. && (beta->r == 1. &&
- beta->i == 0.))) {
- return 0;
- }
- if (*incx > 0) {
- kx = 1;
- } else {
- kx = 1 - (*n - 1) * *incx;
- }
- if (*incy > 0) {
- ky = 1;
- } else {
- ky = 1 - (*n - 1) * *incy;
- }
- if (beta->r != 1. || beta->i != 0.) {
- if (*incy == 1) {
- if (beta->r == 0. && beta->i == 0.) {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = i__;
- y[i__2].r = 0., y[i__2].i = 0.;
- }
- } else {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = i__;
- i__3 = i__;
- z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
- z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
- .r;
- y[i__2].r = z__1.r, y[i__2].i = z__1.i;
- }
- }
- } else {
- iy = ky;
- if (beta->r == 0. && beta->i == 0.) {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = iy;
- y[i__2].r = 0., y[i__2].i = 0.;
- iy += *incy;
- }
- } else {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = iy;
- i__3 = iy;
- z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
- z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
- .r;
- y[i__2].r = z__1.r, y[i__2].i = z__1.i;
- iy += *incy;
- }
- }
- }
- }
- if (alpha->r == 0. && alpha->i == 0.) {
- return 0;
- }
- if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
- kplus1 = *k + 1;
- if (*incx == 1 && *incy == 1) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__2 = j;
- z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
- alpha->r * x[i__2].i + alpha->i * x[i__2].r;
- temp1.r = z__1.r, temp1.i = z__1.i;
- temp2.r = 0., temp2.i = 0.;
- l = kplus1 - j;
- i__2 = 1, i__3 = j - *k;
- i__4 = j - 1;
- for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
- i__2 = i__;
- i__3 = i__;
- i__5 = l + i__ + j * a_dim1;
- z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
- y[i__2].r = z__1.r, y[i__2].i = z__1.i;
- d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
- i__2 = i__;
- z__2.r = z__3.r * x[i__2].r - z__3.i * x[i__2].i, z__2.i =
- z__3.r * x[i__2].i + z__3.i * x[i__2].r;
- z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
- temp2.r = z__1.r, temp2.i = z__1.i;
- }
- i__4 = j;
- i__2 = j;
- i__3 = kplus1 + j * a_dim1;
- d__1 = a[i__3].r;
- z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
- z__2.r = y[i__2].r + z__3.r, z__2.i = y[i__2].i + z__3.i;
- z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
- y[i__4].r = z__1.r, y[i__4].i = z__1.i;
- }
- } else {
- jx = kx;
- jy = ky;
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__4 = jx;
- z__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, z__1.i =
- alpha->r * x[i__4].i + alpha->i * x[i__4].r;
- temp1.r = z__1.r, temp1.i = z__1.i;
- temp2.r = 0., temp2.i = 0.;
- ix = kx;
- iy = ky;
- l = kplus1 - j;
- i__4 = 1, i__2 = j - *k;
- i__3 = j - 1;
- for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
- i__4 = iy;
- i__2 = iy;
- i__5 = l + i__ + j * a_dim1;
- z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
- y[i__4].r = z__1.r, y[i__4].i = z__1.i;
- d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
- i__4 = ix;
- z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
- z__3.r * x[i__4].i + z__3.i * x[i__4].r;
- z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
- temp2.r = z__1.r, temp2.i = z__1.i;
- ix += *incx;
- iy += *incy;
- }
- i__3 = jy;
- i__4 = jy;
- i__2 = kplus1 + j * a_dim1;
- d__1 = a[i__2].r;
- z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
- z__2.r = y[i__4].r + z__3.r, z__2.i = y[i__4].i + z__3.i;
- z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
- y[i__3].r = z__1.r, y[i__3].i = z__1.i;
- jx += *incx;
- jy += *incy;
- if (j > *k) {
- kx += *incx;
- ky += *incy;
- }
- }
- }
- } else {
- if (*incx == 1 && *incy == 1) {
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__3 = j;
- z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
- alpha->r * x[i__3].i + alpha->i * x[i__3].r;
- temp1.r = z__1.r, temp1.i = z__1.i;
- temp2.r = 0., temp2.i = 0.;
- i__3 = j;
- i__4 = j;
- i__2 = j * a_dim1 + 1;
- d__1 = a[i__2].r;
- z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
- z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
- y[i__3].r = z__1.r, y[i__3].i = z__1.i;
- l = 1 - j;
- i__4 = *n, i__2 = j + *k;
- i__3 = min(i__4,i__2);
- for (i__ = j + 1; i__ <= i__3; ++i__) {
- i__4 = i__;
- i__2 = i__;
- i__5 = l + i__ + j * a_dim1;
- z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
- y[i__4].r = z__1.r, y[i__4].i = z__1.i;
- d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
- i__4 = i__;
- z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
- z__3.r * x[i__4].i + z__3.i * x[i__4].r;
- z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
- temp2.r = z__1.r, temp2.i = z__1.i;
- }
- i__3 = j;
- i__4 = j;
- z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
- y[i__3].r = z__1.r, y[i__3].i = z__1.i;
- }
- } else {
- jx = kx;
- jy = ky;
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__3 = jx;
- z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
- alpha->r * x[i__3].i + alpha->i * x[i__3].r;
- temp1.r = z__1.r, temp1.i = z__1.i;
- temp2.r = 0., temp2.i = 0.;
- i__3 = jy;
- i__4 = jy;
- i__2 = j * a_dim1 + 1;
- d__1 = a[i__2].r;
- z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
- z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
- y[i__3].r = z__1.r, y[i__3].i = z__1.i;
- l = 1 - j;
- ix = jx;
- iy = jy;
- i__4 = *n, i__2 = j + *k;
- i__3 = min(i__4,i__2);
- for (i__ = j + 1; i__ <= i__3; ++i__) {
- ix += *incx;
- iy += *incy;
- i__4 = iy;
- i__2 = iy;
- i__5 = l + i__ + j * a_dim1;
- z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
- y[i__4].r = z__1.r, y[i__4].i = z__1.i;
- d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
- i__4 = ix;
- z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
- z__3.r * x[i__4].i + z__3.i * x[i__4].r;
- z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
- temp2.r = z__1.r, temp2.i = z__1.i;
- }
- i__3 = jy;
- i__4 = jy;
- z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
- y[i__3].r = z__1.r, y[i__3].i = z__1.i;
- jx += *incx;
- jy += *incy;
- }
- }
- }
- return 0;
- }
|