swiftnav: zgemv.c Source File

Go to the documentation of this file.
00001 /* zgemv.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 zgemv_(char *trans, integer *m, integer *n, 
00017         doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *
00018         x, integer *incx, doublecomplex *beta, doublecomplex *y, integer *
00019         incy)
00020 {
00021     /* System generated locals */
00022     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
00023     doublecomplex z__1, z__2, z__3;
00024 
00025     /* Builtin functions */
00026     void d_cnjg(doublecomplex *, doublecomplex *);
00027 
00028     /* Local variables */
00029     integer i__, j, ix, iy, jx, jy, kx, ky, info;
00030     doublecomplex temp;
00031     integer lenx, leny;
00032     extern logical lsame_(char *, char *);
00033     extern /* Subroutine */ int xerbla_(char *, integer *);
00034     logical noconj;
00035 
00036 /*     .. Scalar Arguments .. */
00037 /*     .. */
00038 /*     .. Array Arguments .. */
00039 /*     .. */
00040 
00041 /*  Purpose */
00042 /*  ======= */
00043 
00044 /*  ZGEMV  performs one of the matrix-vector operations */
00045 
00046 /*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or */
00047 
00048 /*     y := alpha*conjg( A' )*x + beta*y, */
00049 
00050 /*  where alpha and beta are scalars, x and y are vectors and A is an */
00051 /*  m by n matrix. */
00052 
00053 /*  Arguments */
00054 /*  ========== */
00055 
00056 /*  TRANS  - CHARACTER*1. */
00057 /*           On entry, TRANS specifies the operation to be performed as */
00058 /*           follows: */
00059 
00060 /*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y. */
00061 
00062 /*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y. */
00063 
00064 /*              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y. */
00065 
00066 /*           Unchanged on exit. */
00067 
00068 /*  M      - INTEGER. */
00069 /*           On entry, M specifies the number of rows of the matrix A. */
00070 /*           M must be at least zero. */
00071 /*           Unchanged on exit. */
00072 
00073 /*  N      - INTEGER. */
00074 /*           On entry, N specifies the number of columns of the matrix A. */
00075 /*           N must be at least zero. */
00076 /*           Unchanged on exit. */
00077 
00078 /*  ALPHA  - COMPLEX*16      . */
00079 /*           On entry, ALPHA specifies the scalar alpha. */
00080 /*           Unchanged on exit. */
00081 
00082 /*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ). */
00083 /*           Before entry, the leading m by n part of the array A must */
00084 /*           contain the matrix of coefficients. */
00085 /*           Unchanged on exit. */
00086 
00087 /*  LDA    - INTEGER. */
00088 /*           On entry, LDA specifies the first dimension of A as declared */
00089 /*           in the calling (sub) program. LDA must be at least */
00090 /*           max( 1, m ). */
00091 /*           Unchanged on exit. */
00092 
00093 /*  X      - COMPLEX*16       array of DIMENSION at least */
00094 /*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
00095 /*           and at least */
00096 /*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
00097 /*           Before entry, the incremented array X must contain the */
00098 /*           vector x. */
00099 /*           Unchanged on exit. */
00100 
00101 /*  INCX   - INTEGER. */
00102 /*           On entry, INCX specifies the increment for the elements of */
00103 /*           X. INCX must not be zero. */
00104 /*           Unchanged on exit. */
00105 
00106 /*  BETA   - COMPLEX*16      . */
00107 /*           On entry, BETA specifies the scalar beta. When BETA is */
00108 /*           supplied as zero then Y need not be set on input. */
00109 /*           Unchanged on exit. */
00110 
00111 /*  Y      - COMPLEX*16       array of DIMENSION at least */
00112 /*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
00113 /*           and at least */
00114 /*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
00115 /*           Before entry with BETA non-zero, the incremented array Y */
00116 /*           must contain the vector y. On exit, Y is overwritten by the */
00117 /*           updated vector y. */
00118 
00119 /*  INCY   - INTEGER. */
00120 /*           On entry, INCY specifies the increment for the elements of */
00121 /*           Y. INCY must not be zero. */
00122 /*           Unchanged on exit. */
00123 
00124 
00125 /*  Level 2 Blas routine. */
00126 
00127 /*  -- Written on 22-October-1986. */
00128 /*     Jack Dongarra, Argonne National Lab. */
00129 /*     Jeremy Du Croz, Nag Central Office. */
00130 /*     Sven Hammarling, Nag Central Office. */
00131 /*     Richard Hanson, Sandia National Labs. */
00132 
00133 
00134 /*     .. Parameters .. */
00135 /*     .. */
00136 /*     .. Local Scalars .. */
00137 /*     .. */
00138 /*     .. External Functions .. */
00139 /*     .. */
00140 /*     .. External Subroutines .. */
00141 /*     .. */
00142 /*     .. Intrinsic Functions .. */
00143 /*     .. */
00144 
00145 /*     Test the input parameters. */
00146 
00147     /* Parameter adjustments */
00148     a_dim1 = *lda;
00149     a_offset = 1 + a_dim1;
00150     a -= a_offset;
00151     --x;
00152     --y;
00153 
00154     /* Function Body */
00155     info = 0;
00156     if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
00157             ) {
00158         info = 1;
00159     } else if (*m < 0) {
00160         info = 2;
00161     } else if (*n < 0) {
00162         info = 3;
00163     } else if (*lda < max(1,*m)) {
00164         info = 6;
00165     } else if (*incx == 0) {
00166         info = 8;
00167     } else if (*incy == 0) {
00168         info = 11;
00169     }
00170     if (info != 0) {
00171         xerbla_("ZGEMV ", &info);
00172         return 0;
00173     }
00174 
00175 /*     Quick return if possible. */
00176 
00177     if (*m == 0 || *n == 0 || alpha->r == 0. && alpha->i == 0. && (beta->r == 
00178             1. && beta->i == 0.)) {
00179         return 0;
00180     }
00181 
00182     noconj = lsame_(trans, "T");
00183 
00184 /*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set */
00185 /*     up the start points in  X  and  Y. */
00186 
00187     if (lsame_(trans, "N")) {
00188         lenx = *n;
00189         leny = *m;
00190     } else {
00191         lenx = *m;
00192         leny = *n;
00193     }
00194     if (*incx > 0) {
00195         kx = 1;
00196     } else {
00197         kx = 1 - (lenx - 1) * *incx;
00198     }
00199     if (*incy > 0) {
00200         ky = 1;
00201     } else {
00202         ky = 1 - (leny - 1) * *incy;
00203     }
00204 
00205 /*     Start the operations. In this version the elements of A are */
00206 /*     accessed sequentially with one pass through A. */
00207 
00208 /*     First form  y := beta*y. */
00209 
00210     if (beta->r != 1. || beta->i != 0.) {
00211         if (*incy == 1) {
00212             if (beta->r == 0. && beta->i == 0.) {
00213                 i__1 = leny;
00214                 for (i__ = 1; i__ <= i__1; ++i__) {
00215                     i__2 = i__;
00216                     y[i__2].r = 0., y[i__2].i = 0.;
00217 /* L10: */
00218                 }
00219             } else {
00220                 i__1 = leny;
00221                 for (i__ = 1; i__ <= i__1; ++i__) {
00222                     i__2 = i__;
00223                     i__3 = i__;
00224                     z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
00225                             z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
00226                             .r;
00227                     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00228 /* L20: */
00229                 }
00230             }
00231         } else {
00232             iy = ky;
00233             if (beta->r == 0. && beta->i == 0.) {
00234                 i__1 = leny;
00235                 for (i__ = 1; i__ <= i__1; ++i__) {
00236                     i__2 = iy;
00237                     y[i__2].r = 0., y[i__2].i = 0.;
00238                     iy += *incy;
00239 /* L30: */
00240                 }
00241             } else {
00242                 i__1 = leny;
00243                 for (i__ = 1; i__ <= i__1; ++i__) {
00244                     i__2 = iy;
00245                     i__3 = iy;
00246                     z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 
00247                             z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
00248                             .r;
00249                     y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00250                     iy += *incy;
00251 /* L40: */
00252                 }
00253             }
00254         }
00255     }
00256     if (alpha->r == 0. && alpha->i == 0.) {
00257         return 0;
00258     }
00259     if (lsame_(trans, "N")) {
00260 
00261 /*        Form  y := alpha*A*x + y. */
00262 
00263         jx = kx;
00264         if (*incy == 1) {
00265             i__1 = *n;
00266             for (j = 1; j <= i__1; ++j) {
00267                 i__2 = jx;
00268                 if (x[i__2].r != 0. || x[i__2].i != 0.) {
00269                     i__2 = jx;
00270                     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00271                             z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00272                             .r;
00273                     temp.r = z__1.r, temp.i = z__1.i;
00274                     i__2 = *m;
00275                     for (i__ = 1; i__ <= i__2; ++i__) {
00276                         i__3 = i__;
00277                         i__4 = i__;
00278                         i__5 = i__ + j * a_dim1;
00279                         z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i, 
00280                                 z__2.i = temp.r * a[i__5].i + temp.i * a[i__5]
00281                                 .r;
00282                         z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + 
00283                                 z__2.i;
00284                         y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00285 /* L50: */
00286                     }
00287                 }
00288                 jx += *incx;
00289 /* L60: */
00290             }
00291         } else {
00292             i__1 = *n;
00293             for (j = 1; j <= i__1; ++j) {
00294                 i__2 = jx;
00295                 if (x[i__2].r != 0. || x[i__2].i != 0.) {
00296                     i__2 = jx;
00297                     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00298                             z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00299                             .r;
00300                     temp.r = z__1.r, temp.i = z__1.i;
00301                     iy = ky;
00302                     i__2 = *m;
00303                     for (i__ = 1; i__ <= i__2; ++i__) {
00304                         i__3 = iy;
00305                         i__4 = iy;
00306                         i__5 = i__ + j * a_dim1;
00307                         z__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i, 
00308                                 z__2.i = temp.r * a[i__5].i + temp.i * a[i__5]
00309                                 .r;
00310                         z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + 
00311                                 z__2.i;
00312                         y[i__3].r = z__1.r, y[i__3].i = z__1.i;
00313                         iy += *incy;
00314 /* L70: */
00315                     }
00316                 }
00317                 jx += *incx;
00318 /* L80: */
00319             }
00320         }
00321     } else {
00322 
00323 /*        Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y. */
00324 
00325         jy = ky;
00326         if (*incx == 1) {
00327             i__1 = *n;
00328             for (j = 1; j <= i__1; ++j) {
00329                 temp.r = 0., temp.i = 0.;
00330                 if (noconj) {
00331                     i__2 = *m;
00332                     for (i__ = 1; i__ <= i__2; ++i__) {
00333                         i__3 = i__ + j * a_dim1;
00334                         i__4 = i__;
00335                         z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4]
00336                                 .i, z__2.i = a[i__3].r * x[i__4].i + a[i__3]
00337                                 .i * x[i__4].r;
00338                         z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
00339                         temp.r = z__1.r, temp.i = z__1.i;
00340 /* L90: */
00341                     }
00342                 } else {
00343                     i__2 = *m;
00344                     for (i__ = 1; i__ <= i__2; ++i__) {
00345                         d_cnjg(&z__3, &a[i__ + j * a_dim1]);
00346                         i__3 = i__;
00347                         z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, 
00348                                 z__2.i = z__3.r * x[i__3].i + z__3.i * x[i__3]
00349                                 .r;
00350                         z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
00351                         temp.r = z__1.r, temp.i = z__1.i;
00352 /* L100: */
00353                     }
00354                 }
00355                 i__2 = jy;
00356                 i__3 = jy;
00357                 z__2.r = alpha->r * temp.r - alpha->i * temp.i, z__2.i = 
00358                         alpha->r * temp.i + alpha->i * temp.r;
00359                 z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
00360                 y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00361                 jy += *incy;
00362 /* L110: */
00363             }
00364         } else {
00365             i__1 = *n;
00366             for (j = 1; j <= i__1; ++j) {
00367                 temp.r = 0., temp.i = 0.;
00368                 ix = kx;
00369                 if (noconj) {
00370                     i__2 = *m;
00371                     for (i__ = 1; i__ <= i__2; ++i__) {
00372                         i__3 = i__ + j * a_dim1;
00373                         i__4 = ix;
00374                         z__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4]
00375                                 .i, z__2.i = a[i__3].r * x[i__4].i + a[i__3]
00376                                 .i * x[i__4].r;
00377                         z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
00378                         temp.r = z__1.r, temp.i = z__1.i;
00379                         ix += *incx;
00380 /* L120: */
00381                     }
00382                 } else {
00383                     i__2 = *m;
00384                     for (i__ = 1; i__ <= i__2; ++i__) {
00385                         d_cnjg(&z__3, &a[i__ + j * a_dim1]);
00386                         i__3 = ix;
00387                         z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, 
00388                                 z__2.i = z__3.r * x[i__3].i + z__3.i * x[i__3]
00389                                 .r;
00390                         z__1.r = temp.r + z__2.r, z__1.i = temp.i + z__2.i;
00391                         temp.r = z__1.r, temp.i = z__1.i;
00392                         ix += *incx;
00393 /* L130: */
00394                     }
00395                 }
00396                 i__2 = jy;
00397                 i__3 = jy;
00398                 z__2.r = alpha->r * temp.r - alpha->i * temp.i, z__2.i = 
00399                         alpha->r * temp.i + alpha->i * temp.r;
00400                 z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
00401                 y[i__2].r = z__1.r, y[i__2].i = z__1.i;
00402                 jy += *incy;
00403 /* L140: */
00404             }
00405         }
00406     }
00407 
00408     return 0;
00409 
00410 /*     End of ZGEMV . */
00411 
00412 } /* zgemv_ */