123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368 |
- #include "datatypes.h"
- int ssbmv_(char *uplo, integer *n, integer *k, real *alpha,
- real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
- integer *incy, ftnlen uplo_len)
- {
-
- integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
-
- integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
- real 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_("SSBMV ", &info, (ftnlen)6);
- return 0;
- }
- if (*n == 0 || (*alpha == 0.f && *beta == 1.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 != 1.f) {
- if (*incy == 1) {
- if (*beta == 0.f) {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[i__] = 0.f;
- }
- } else {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[i__] = *beta * y[i__];
- }
- }
- } else {
- iy = ky;
- if (*beta == 0.f) {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[iy] = 0.f;
- iy += *incy;
- }
- } else {
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- y[iy] = *beta * y[iy];
- iy += *incy;
- }
- }
- }
- }
- if (*alpha == 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) {
- temp1 = *alpha * x[j];
- temp2 = 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__) {
- y[i__] += temp1 * a[l + i__ + j * a_dim1];
- temp2 += a[l + i__ + j * a_dim1] * x[i__];
- }
- y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
- }
- } else {
- jx = kx;
- jy = ky;
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- temp1 = *alpha * x[jx];
- temp2 = 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__) {
- y[iy] += temp1 * a[l + i__ + j * a_dim1];
- temp2 += a[l + i__ + j * a_dim1] * x[ix];
- ix += *incx;
- iy += *incy;
- }
- y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha *
- temp2;
- 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) {
- temp1 = *alpha * x[j];
- temp2 = 0.f;
- y[j] += temp1 * a[j * a_dim1 + 1];
- 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__) {
- y[i__] += temp1 * a[l + i__ + j * a_dim1];
- temp2 += a[l + i__ + j * a_dim1] * x[i__];
- }
- y[j] += *alpha * temp2;
- }
- } else {
- jx = kx;
- jy = ky;
- i__1 = *n;
- for (j = 1; j <= i__1; ++j) {
- temp1 = *alpha * x[jx];
- temp2 = 0.f;
- y[jy] += temp1 * a[j * a_dim1 + 1];
- 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;
- y[iy] += temp1 * a[l + i__ + j * a_dim1];
- temp2 += a[l + i__ + j * a_dim1] * x[ix];
- }
- y[jy] += *alpha * temp2;
- jx += *incx;
- jy += *incy;
- }
- }
- }
- return 0;
- }
|