zhbmv.c
Go to the documentation of this file.
00001 /* zhbmv.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int zhbmv_(char *uplo, integer *n, integer *k, doublecomplex 
00017         *alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *
00018         incx, doublecomplex *beta, doublecomplex *y, integer *incy)
00019 {
00020     /* System generated locals */
00021     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
00022     doublereal d__1;
00023     doublecomplex z__1, z__2, z__3, z__4;
00024 
00025     /* Builtin functions */
00026     void d_cnjg(doublecomplex *, doublecomplex *);
00027 
00028     /* Local variables */
00029     integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
00030     doublecomplex temp1, temp2;
00031     extern logical lsame_(char *, char *);
00032     integer kplus1;
00033     extern /* Subroutine */ int xerbla_(char *, integer *);
00034 
00035 /*     .. Scalar Arguments .. */
00036 /*     .. */
00037 /*     .. Array Arguments .. */
00038 /*     .. */
00039 
00040 /*  Purpose */
00041 /*  ======= */
00042 
00043 /*  ZHBMV  performs the matrix-vector  operation */
00044 
00045 /*     y := alpha*A*x + beta*y, */
00046 
00047 /*  where alpha and beta are scalars, x and y are n element vectors and */
00048 /*  A is an n by n hermitian band matrix, with k super-diagonals. */
00049 
00050 /*  Arguments */
00051 /*  ========== */
00052 
00053 /*  UPLO   - CHARACTER*1. */
00054 /*           On entry, UPLO specifies whether the upper or lower */
00055 /*           triangular part of the band matrix A is being supplied as */
00056 /*           follows: */
00057 
00058 /*              UPLO = 'U' or 'u'   The upper triangular part of A is */
00059 /*                                  being supplied. */
00060 
00061 /*              UPLO = 'L' or 'l'   The lower triangular part of A is */
00062 /*                                  being supplied. */
00063 
00064 /*           Unchanged on exit. */
00065 
00066 /*  N      - INTEGER. */
00067 /*           On entry, N specifies the order of the matrix A. */
00068 /*           N must be at least zero. */
00069 /*           Unchanged on exit. */
00070 
00071 /*  K      - INTEGER. */
00072 /*           On entry, K specifies the number of super-diagonals of the */
00073 /*           matrix A. K must satisfy  0 .le. K. */
00074 /*           Unchanged on exit. */
00075 
00076 /*  ALPHA  - COMPLEX*16      . */
00077 /*           On entry, ALPHA specifies the scalar alpha. */
00078 /*           Unchanged on exit. */
00079 
00080 /*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ). */
00081 /*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
00082 /*           by n part of the array A must contain the upper triangular */
00083 /*           band part of the hermitian matrix, supplied column by */
00084 /*           column, with the leading diagonal of the matrix in row */
00085 /*           ( k + 1 ) of the array, the first super-diagonal starting at */
00086 /*           position 2 in row k, and so on. The top left k by k triangle */
00087 /*           of the array A is not referenced. */
00088 /*           The following program segment will transfer the upper */
00089 /*           triangular part of a hermitian band matrix from conventional */
00090 /*           full matrix storage to band storage: */
00091 
00092 /*                 DO 20, J = 1, N */
00093 /*                    M = K + 1 - J */
00094 /*                    DO 10, I = MAX( 1, J - K ), J */
00095 /*                       A( M + I, J ) = matrix( I, J ) */
00096 /*              10    CONTINUE */
00097 /*              20 CONTINUE */
00098 
00099 /*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
00100 /*           by n part of the array A must contain the lower triangular */
00101 /*           band part of the hermitian matrix, supplied column by */
00102 /*           column, with the leading diagonal of the matrix in row 1 of */
00103 /*           the array, the first sub-diagonal starting at position 1 in */
00104 /*           row 2, and so on. The bottom right k by k triangle of the */
00105 /*           array A is not referenced. */
00106 /*           The following program segment will transfer the lower */
00107 /*           triangular part of a hermitian band matrix from conventional */
00108 /*           full matrix storage to band storage: */
00109 
00110 /*                 DO 20, J = 1, N */
00111 /*                    M = 1 - J */
00112 /*                    DO 10, I = J, MIN( N, J + K ) */
00113 /*                       A( M + I, J ) = matrix( I, J ) */
00114 /*              10    CONTINUE */
00115 /*              20 CONTINUE */
00116 
00117 /*           Note that the imaginary parts of the diagonal elements need */
00118 /*           not be set and are assumed to be zero. */
00119 /*           Unchanged on exit. */
00120 
00121 /*  LDA    - INTEGER. */
00122 /*           On entry, LDA specifies the first dimension of A as declared */
00123 /*           in the calling (sub) program. LDA must be at least */
00124 /*           ( k + 1 ). */
00125 /*           Unchanged on exit. */
00126 
00127 /*  X      - COMPLEX*16       array of DIMENSION at least */
00128 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00129 /*           Before entry, the incremented array X must contain the */
00130 /*           vector x. */
00131 /*           Unchanged on exit. */
00132 
00133 /*  INCX   - INTEGER. */
00134 /*           On entry, INCX specifies the increment for the elements of */
00135 /*           X. INCX must not be zero. */
00136 /*           Unchanged on exit. */
00137 
00138 /*  BETA   - COMPLEX*16      . */
00139 /*           On entry, BETA specifies the scalar beta. */
00140 /*           Unchanged on exit. */
00141 
00142 /*  Y      - COMPLEX*16       array of DIMENSION at least */
00143 /*           ( 1 + ( n - 1 )*abs( INCY ) ). */
00144 /*           Before entry, the incremented array Y must contain the */
00145 /*           vector y. On exit, Y is overwritten by the updated vector y. */
00146 
00147 /*  INCY   - INTEGER. */
00148 /*           On entry, INCY specifies the increment for the elements of */
00149 /*           Y. INCY must not be zero. */
00150 /*           Unchanged on exit. */
00151 
00152 
00153 /*  Level 2 Blas routine. */
00154 
00155 /*  -- Written on 22-October-1986. */
00156 /*     Jack Dongarra, Argonne National Lab. */
00157 /*     Jeremy Du Croz, Nag Central Office. */
00158 /*     Sven Hammarling, Nag Central Office. */
00159 /*     Richard Hanson, Sandia National Labs. */
00160 
00161 
00162 /*     .. Parameters .. */
00163 /*     .. */
00164 /*     .. Local Scalars .. */
00165 /*     .. */
00166 /*     .. External Functions .. */
00167 /*     .. */
00168 /*     .. External Subroutines .. */
00169 /*     .. */
00170 /*     .. Intrinsic Functions .. */
00171 /*     .. */
00172 
00173 /*     Test the input parameters. */
00174 
00175     /* Parameter adjustments */
00176     a_dim1 = *lda;
00177     a_offset = 1 + a_dim1;
00178     a -= a_offset;
00179     --x;
00180     --y;
00181 
00182     /* Function Body */
00183     info = 0;
00184     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00185         info = 1;
00186     } else if (*n < 0) {
00187         info = 2;
00188     } else if (*k < 0) {
00189         info = 3;
00190     } else if (*lda < *k + 1) {
00191         info = 6;
00192     } else if (*incx == 0) {
00193         info = 8;
00194     } else if (*incy == 0) {
00195         info = 11;
00196     }
00197     if (info != 0) {
00198         xerbla_("ZHBMV ", &info);
00199         return 0;
00200     }
00201 
00202 /*     Quick return if possible. */
00203 
00204     if (*n == 0 || alpha->r == 0. && alpha->i == 0. && (beta->r == 1. && 
00205             beta->i == 0.)) {
00206         return 0;
00207     }
00208 
00209 /*     Set up the start points in  X  and  Y. */
00210 
00211     if (*incx > 0) {
00212         kx = 1;
00213     } else {
00214         kx = 1 - (*n - 1) * *incx;
00215     }
00216     if (*incy > 0) {
00217         ky = 1;
00218     } else {
00219         ky = 1 - (*n - 1) * *incy;
00220     }
00221 
00222 /*     Start the operations. In this version the elements of the array A */
00223 /*     are accessed sequentially with one pass through A. */
00224 
00225 /*     First form  y := beta*y. */
00226 
00227     if (beta->r != 1. || beta->i != 0.) {
00228         if (*incy == 1) {
00229             if (beta->r == 0. && beta->i == 0.) {
00230                 i__1 = *n;
00231                 for (i__ = 1; i__ <= i__1; ++i__) {
00232                     i__2 = i__;
00233                     y[i__2].r = 0., y[i__2].i = 0.;
00234 /* L10: */
00235                 }
00236             } else {
00237                 i__1 = *n;
00238                 for (i__ = 1; i__ <= i__1; ++i__) {
00239                     i__2 = i__;
00240                     i__3 = i__;
00241                     z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
00242                             z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
00243                             .r;
00244                     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00245 /* L20: */
00246                 }
00247             }
00248         } else {
00249             iy = ky;
00250             if (beta->r == 0. && beta->i == 0.) {
00251                 i__1 = *n;
00252                 for (i__ = 1; i__ <= i__1; ++i__) {
00253                     i__2 = iy;
00254                     y[i__2].r = 0., y[i__2].i = 0.;
00255                     iy += *incy;
00256 /* L30: */
00257                 }
00258             } else {
00259                 i__1 = *n;
00260                 for (i__ = 1; i__ <= i__1; ++i__) {
00261                     i__2 = iy;
00262                     i__3 = iy;
00263                     z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
00264                             z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
00265                             .r;
00266                     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00267                     iy += *incy;
00268 /* L40: */
00269                 }
00270             }
00271         }
00272     }
00273     if (alpha->r == 0. && alpha->i == 0.) {
00274         return 0;
00275     }
00276     if (lsame_(uplo, "U")) {
00277 
00278 /*        Form  y  when upper triangle of A is stored. */
00279 
00280         kplus1 = *k + 1;
00281         if (*incx == 1 && *incy == 1) {
00282             i__1 = *n;
00283             for (j = 1; j <= i__1; ++j) {
00284                 i__2 = j;
00285                 z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
00286                          alpha->r * x[i__2].i + alpha->i * x[i__2].r;
00287                 temp1.r = z__1.r, temp1.i = z__1.i;
00288                 temp2.r = 0., temp2.i = 0.;
00289                 l = kplus1 - j;
00290 /* Computing MAX */
00291                 i__2 = 1, i__3 = j - *k;
00292                 i__4 = j - 1;
00293                 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
00294                     i__2 = i__;
00295                     i__3 = i__;
00296                     i__5 = l + i__ + j * a_dim1;
00297                     z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
00298                             z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
00299                             .r;
00300                     z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
00301                     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00302                     d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
00303                     i__2 = i__;
00304                     z__2.r = z__3.r * x[i__2].r - z__3.i * x[i__2].i, z__2.i =
00305                              z__3.r * x[i__2].i + z__3.i * x[i__2].r;
00306                     z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
00307                     temp2.r = z__1.r, temp2.i = z__1.i;
00308 /* L50: */
00309                 }
00310                 i__4 = j;
00311                 i__2 = j;
00312                 i__3 = kplus1 + j * a_dim1;
00313                 d__1 = a[i__3].r;
00314                 z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
00315                 z__2.r = y[i__2].r + z__3.r, z__2.i = y[i__2].i + z__3.i;
00316                 z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 
00317                         alpha->r * temp2.i + alpha->i * temp2.r;
00318                 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
00319                 y[i__4].r = z__1.r, y[i__4].i = z__1.i;
00320 /* L60: */
00321             }
00322         } else {
00323             jx = kx;
00324             jy = ky;
00325             i__1 = *n;
00326             for (j = 1; j <= i__1; ++j) {
00327                 i__4 = jx;
00328                 z__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, z__1.i =
00329                          alpha->r * x[i__4].i + alpha->i * x[i__4].r;
00330                 temp1.r = z__1.r, temp1.i = z__1.i;
00331                 temp2.r = 0., temp2.i = 0.;
00332                 ix = kx;
00333                 iy = ky;
00334                 l = kplus1 - j;
00335 /* Computing MAX */
00336                 i__4 = 1, i__2 = j - *k;
00337                 i__3 = j - 1;
00338                 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
00339                     i__4 = iy;
00340                     i__2 = iy;
00341                     i__5 = l + i__ + j * a_dim1;
00342                     z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
00343                             z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
00344                             .r;
00345                     z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
00346                     y[i__4].r = z__1.r, y[i__4].i = z__1.i;
00347                     d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
00348                     i__4 = ix;
00349                     z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
00350                              z__3.r * x[i__4].i + z__3.i * x[i__4].r;
00351                     z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
00352                     temp2.r = z__1.r, temp2.i = z__1.i;
00353                     ix += *incx;
00354                     iy += *incy;
00355 /* L70: */
00356                 }
00357                 i__3 = jy;
00358                 i__4 = jy;
00359                 i__2 = kplus1 + j * a_dim1;
00360                 d__1 = a[i__2].r;
00361                 z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
00362                 z__2.r = y[i__4].r + z__3.r, z__2.i = y[i__4].i + z__3.i;
00363                 z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 
00364                         alpha->r * temp2.i + alpha->i * temp2.r;
00365                 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
00366                 y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00367                 jx += *incx;
00368                 jy += *incy;
00369                 if (j > *k) {
00370                     kx += *incx;
00371                     ky += *incy;
00372                 }
00373 /* L80: */
00374             }
00375         }
00376     } else {
00377 
00378 /*        Form  y  when lower triangle of A is stored. */
00379 
00380         if (*incx == 1 && *incy == 1) {
00381             i__1 = *n;
00382             for (j = 1; j <= i__1; ++j) {
00383                 i__3 = j;
00384                 z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
00385                          alpha->r * x[i__3].i + alpha->i * x[i__3].r;
00386                 temp1.r = z__1.r, temp1.i = z__1.i;
00387                 temp2.r = 0., temp2.i = 0.;
00388                 i__3 = j;
00389                 i__4 = j;
00390                 i__2 = j * a_dim1 + 1;
00391                 d__1 = a[i__2].r;
00392                 z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
00393                 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
00394                 y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00395                 l = 1 - j;
00396 /* Computing MIN */
00397                 i__4 = *n, i__2 = j + *k;
00398                 i__3 = min(i__4,i__2);
00399                 for (i__ = j + 1; i__ <= i__3; ++i__) {
00400                     i__4 = i__;
00401                     i__2 = i__;
00402                     i__5 = l + i__ + j * a_dim1;
00403                     z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
00404                             z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
00405                             .r;
00406                     z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
00407                     y[i__4].r = z__1.r, y[i__4].i = z__1.i;
00408                     d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
00409                     i__4 = i__;
00410                     z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
00411                              z__3.r * x[i__4].i + z__3.i * x[i__4].r;
00412                     z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
00413                     temp2.r = z__1.r, temp2.i = z__1.i;
00414 /* L90: */
00415                 }
00416                 i__3 = j;
00417                 i__4 = j;
00418                 z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 
00419                         alpha->r * temp2.i + alpha->i * temp2.r;
00420                 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
00421                 y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00422 /* L100: */
00423             }
00424         } else {
00425             jx = kx;
00426             jy = ky;
00427             i__1 = *n;
00428             for (j = 1; j <= i__1; ++j) {
00429                 i__3 = jx;
00430                 z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
00431                          alpha->r * x[i__3].i + alpha->i * x[i__3].r;
00432                 temp1.r = z__1.r, temp1.i = z__1.i;
00433                 temp2.r = 0., temp2.i = 0.;
00434                 i__3 = jy;
00435                 i__4 = jy;
00436                 i__2 = j * a_dim1 + 1;
00437                 d__1 = a[i__2].r;
00438                 z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
00439                 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
00440                 y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00441                 l = 1 - j;
00442                 ix = jx;
00443                 iy = jy;
00444 /* Computing MIN */
00445                 i__4 = *n, i__2 = j + *k;
00446                 i__3 = min(i__4,i__2);
00447                 for (i__ = j + 1; i__ <= i__3; ++i__) {
00448                     ix += *incx;
00449                     iy += *incy;
00450                     i__4 = iy;
00451                     i__2 = iy;
00452                     i__5 = l + i__ + j * a_dim1;
00453                     z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i, 
00454                             z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
00455                             .r;
00456                     z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
00457                     y[i__4].r = z__1.r, y[i__4].i = z__1.i;
00458                     d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
00459                     i__4 = ix;
00460                     z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
00461                              z__3.r * x[i__4].i + z__3.i * x[i__4].r;
00462                     z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
00463                     temp2.r = z__1.r, temp2.i = z__1.i;
00464 /* L110: */
00465                 }
00466                 i__3 = jy;
00467                 i__4 = jy;
00468                 z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 
00469                         alpha->r * temp2.i + alpha->i * temp2.r;
00470                 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
00471                 y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00472                 jx += *incx;
00473                 jy += *incy;
00474 /* L120: */
00475             }
00476         }
00477     }
00478 
00479     return 0;
00480 
00481 /*     End of ZHBMV . */
00482 
00483 } /* zhbmv_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:36