123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487 |
- #include "datatypes.h"
- int chbmv_(char *uplo, integer *n, integer *k, complex *
- alpha, complex *a, integer *lda, complex *x, integer *incx, complex *
- beta, complex *y, integer *incy, ftnlen uplo_len)
- {
-
- integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
- real r__1;
- complex q__1, q__2, q__3, q__4;
-
- void r_cnjg(complex *, complex *);
-
- integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
- complex 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_("CHBMV ", &info, (ftnlen)6);
- return 0;
- }
- if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f && (beta->r == 1.f &&
- beta->i == 0.f))) {
- 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.f || beta->i != 0.f) {
- if (*incy == 1) {
- if (beta->r == 0.f && beta->i == 0.f) {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = i__;
- y[i__2].r = 0.f, y[i__2].i = 0.f;
- }
- } else {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = i__;
- i__3 = i__;
- q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
- q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
- .r;
- y[i__2].r = q__1.r, y[i__2].i = q__1.i;
- }
- }
- } else {
- iy = ky;
- if (beta->r == 0.f && beta->i == 0.f) {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = iy;
- y[i__2].r = 0.f, y[i__2].i = 0.f;
- iy += *incy;
- }
- } else {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = iy;
- i__3 = iy;
- q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
- q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
- .r;
- y[i__2].r = q__1.r, y[i__2].i = q__1.i;
- iy += *incy;
- }
- }
- }
- }
- if (alpha->r == 0.f && alpha->i == 0.f) {
- 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;
- q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
- alpha->r * x[i__2].i + alpha->i * x[i__2].r;
- temp1.r = q__1.r, temp1.i = q__1.i;
- temp2.r = 0.f, temp2.i = 0.f;
- 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;
- q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
- y[i__2].r = q__1.r, y[i__2].i = q__1.i;
- r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
- i__2 = i__;
- q__2.r = q__3.r * x[i__2].r - q__3.i * x[i__2].i, q__2.i =
- q__3.r * x[i__2].i + q__3.i * x[i__2].r;
- q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
- temp2.r = q__1.r, temp2.i = q__1.i;
- }
- i__4 = j;
- i__2 = j;
- i__3 = kplus1 + j * a_dim1;
- r__1 = a[i__3].r;
- q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
- q__2.r = y[i__2].r + q__3.r, q__2.i = y[i__2].i + q__3.i;
- q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
- y[i__4].r = q__1.r, y[i__4].i = q__1.i;
- }
- } else {
- jx = kx;
- jy = ky;
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__4 = jx;
- q__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, q__1.i =
- alpha->r * x[i__4].i + alpha->i * x[i__4].r;
- temp1.r = q__1.r, temp1.i = q__1.i;
- temp2.r = 0.f, temp2.i = 0.f;
- 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;
- q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i;
- y[i__4].r = q__1.r, y[i__4].i = q__1.i;
- r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
- i__4 = ix;
- q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i =
- q__3.r * x[i__4].i + q__3.i * x[i__4].r;
- q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
- temp2.r = q__1.r, temp2.i = q__1.i;
- ix += *incx;
- iy += *incy;
- }
- i__3 = jy;
- i__4 = jy;
- i__2 = kplus1 + j * a_dim1;
- r__1 = a[i__2].r;
- q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
- q__2.r = y[i__4].r + q__3.r, q__2.i = y[i__4].i + q__3.i;
- q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
- y[i__3].r = q__1.r, y[i__3].i = q__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;
- q__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, q__1.i =
- alpha->r * x[i__3].i + alpha->i * x[i__3].r;
- temp1.r = q__1.r, temp1.i = q__1.i;
- temp2.r = 0.f, temp2.i = 0.f;
- i__3 = j;
- i__4 = j;
- i__2 = j * a_dim1 + 1;
- r__1 = a[i__2].r;
- q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
- q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
- y[i__3].r = q__1.r, y[i__3].i = q__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;
- q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i;
- y[i__4].r = q__1.r, y[i__4].i = q__1.i;
- r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
- i__4 = i__;
- q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i =
- q__3.r * x[i__4].i + q__3.i * x[i__4].r;
- q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
- temp2.r = q__1.r, temp2.i = q__1.i;
- }
- i__3 = j;
- i__4 = j;
- q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
- y[i__3].r = q__1.r, y[i__3].i = q__1.i;
- }
- } else {
- jx = kx;
- jy = ky;
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- i__3 = jx;
- q__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, q__1.i =
- alpha->r * x[i__3].i + alpha->i * x[i__3].r;
- temp1.r = q__1.r, temp1.i = q__1.i;
- temp2.r = 0.f, temp2.i = 0.f;
- i__3 = jy;
- i__4 = jy;
- i__2 = j * a_dim1 + 1;
- r__1 = a[i__2].r;
- q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
- q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
- y[i__3].r = q__1.r, y[i__3].i = q__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;
- q__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
- q__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
- .r;
- q__1.r = y[i__2].r + q__2.r, q__1.i = y[i__2].i + q__2.i;
- y[i__4].r = q__1.r, y[i__4].i = q__1.i;
- r_cnjg(&q__3, &a[l + i__ + j * a_dim1]);
- i__4 = ix;
- q__2.r = q__3.r * x[i__4].r - q__3.i * x[i__4].i, q__2.i =
- q__3.r * x[i__4].i + q__3.i * x[i__4].r;
- q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
- temp2.r = q__1.r, temp2.i = q__1.i;
- }
- i__3 = jy;
- i__4 = jy;
- q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i =
- alpha->r * temp2.i + alpha->i * temp2.r;
- q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
- y[i__3].r = q__1.r, y[i__3].i = q__1.i;
- jx += *incx;
- jy += *incy;
- }
- }
- }
- return 0;
- }
|