ctpmv.c
Go to the documentation of this file.
00001 /* ctpmv.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 ctpmv_(char *uplo, char *trans, char *diag, integer *n, 
00017         complex *ap, complex *x, integer *incx)
00018 {
00019     /* System generated locals */
00020     integer i__1, i__2, i__3, i__4, i__5;
00021     complex q__1, q__2, q__3;
00022 
00023     /* Builtin functions */
00024     void r_cnjg(complex *, complex *);
00025 
00026     /* Local variables */
00027     integer i__, j, k, kk, ix, jx, kx, info;
00028     complex temp;
00029     extern logical lsame_(char *, char *);
00030     extern /* Subroutine */ int xerbla_(char *, integer *);
00031     logical noconj, nounit;
00032 
00033 /*     .. Scalar Arguments .. */
00034 /*     .. */
00035 /*     .. Array Arguments .. */
00036 /*     .. */
00037 
00038 /*  Purpose */
00039 /*  ======= */
00040 
00041 /*  CTPMV  performs one of the matrix-vector operations */
00042 
00043 /*     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*x, */
00044 
00045 /*  where x is an n element vector and  A is an n by n unit, or non-unit, */
00046 /*  upper or lower triangular matrix, supplied in packed form. */
00047 
00048 /*  Arguments */
00049 /*  ========== */
00050 
00051 /*  UPLO   - CHARACTER*1. */
00052 /*           On entry, UPLO specifies whether the matrix is an upper or */
00053 /*           lower triangular matrix as follows: */
00054 
00055 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
00056 
00057 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
00058 
00059 /*           Unchanged on exit. */
00060 
00061 /*  TRANS  - CHARACTER*1. */
00062 /*           On entry, TRANS specifies the operation to be performed as */
00063 /*           follows: */
00064 
00065 /*              TRANS = 'N' or 'n'   x := A*x. */
00066 
00067 /*              TRANS = 'T' or 't'   x := A'*x. */
00068 
00069 /*              TRANS = 'C' or 'c'   x := conjg( A' )*x. */
00070 
00071 /*           Unchanged on exit. */
00072 
00073 /*  DIAG   - CHARACTER*1. */
00074 /*           On entry, DIAG specifies whether or not A is unit */
00075 /*           triangular as follows: */
00076 
00077 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00078 
00079 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00080 /*                                  triangular. */
00081 
00082 /*           Unchanged on exit. */
00083 
00084 /*  N      - INTEGER. */
00085 /*           On entry, N specifies the order of the matrix A. */
00086 /*           N must be at least zero. */
00087 /*           Unchanged on exit. */
00088 
00089 /*  AP     - COMPLEX          array of DIMENSION at least */
00090 /*           ( ( n*( n + 1 ) )/2 ). */
00091 /*           Before entry with  UPLO = 'U' or 'u', the array AP must */
00092 /*           contain the upper triangular matrix packed sequentially, */
00093 /*           column by column, so that AP( 1 ) contains a( 1, 1 ), */
00094 /*           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) */
00095 /*           respectively, and so on. */
00096 /*           Before entry with UPLO = 'L' or 'l', the array AP must */
00097 /*           contain the lower triangular matrix packed sequentially, */
00098 /*           column by column, so that AP( 1 ) contains a( 1, 1 ), */
00099 /*           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) */
00100 /*           respectively, and so on. */
00101 /*           Note that when  DIAG = 'U' or 'u', the diagonal elements of */
00102 /*           A are not referenced, but are assumed to be unity. */
00103 /*           Unchanged on exit. */
00104 
00105 /*  X      - COMPLEX          array of dimension at least */
00106 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00107 /*           Before entry, the incremented array X must contain the n */
00108 /*           element vector x. On exit, X is overwritten with the */
00109 /*           tranformed vector x. */
00110 
00111 /*  INCX   - INTEGER. */
00112 /*           On entry, INCX specifies the increment for the elements of */
00113 /*           X. INCX must not be zero. */
00114 /*           Unchanged on exit. */
00115 
00116 
00117 /*  Level 2 Blas routine. */
00118 
00119 /*  -- Written on 22-October-1986. */
00120 /*     Jack Dongarra, Argonne National Lab. */
00121 /*     Jeremy Du Croz, Nag Central Office. */
00122 /*     Sven Hammarling, Nag Central Office. */
00123 /*     Richard Hanson, Sandia National Labs. */
00124 
00125 
00126 /*     .. Parameters .. */
00127 /*     .. */
00128 /*     .. Local Scalars .. */
00129 /*     .. */
00130 /*     .. External Functions .. */
00131 /*     .. */
00132 /*     .. External Subroutines .. */
00133 /*     .. */
00134 /*     .. Intrinsic Functions .. */
00135 /*     .. */
00136 
00137 /*     Test the input parameters. */
00138 
00139     /* Parameter adjustments */
00140     --x;
00141     --ap;
00142 
00143     /* Function Body */
00144     info = 0;
00145     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00146         info = 1;
00147     } else if (! lsame_(trans, "N") && ! lsame_(trans, 
00148             "T") && ! lsame_(trans, "C")) {
00149         info = 2;
00150     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00151             "N")) {
00152         info = 3;
00153     } else if (*n < 0) {
00154         info = 4;
00155     } else if (*incx == 0) {
00156         info = 7;
00157     }
00158     if (info != 0) {
00159         xerbla_("CTPMV ", &info);
00160         return 0;
00161     }
00162 
00163 /*     Quick return if possible. */
00164 
00165     if (*n == 0) {
00166         return 0;
00167     }
00168 
00169     noconj = lsame_(trans, "T");
00170     nounit = lsame_(diag, "N");
00171 
00172 /*     Set up the start point in X if the increment is not unity. This */
00173 /*     will be  ( N - 1 )*INCX  too small for descending loops. */
00174 
00175     if (*incx <= 0) {
00176         kx = 1 - (*n - 1) * *incx;
00177     } else if (*incx != 1) {
00178         kx = 1;
00179     }
00180 
00181 /*     Start the operations. In this version the elements of AP are */
00182 /*     accessed sequentially with one pass through AP. */
00183 
00184     if (lsame_(trans, "N")) {
00185 
00186 /*        Form  x:= A*x. */
00187 
00188         if (lsame_(uplo, "U")) {
00189             kk = 1;
00190             if (*incx == 1) {
00191                 i__1 = *n;
00192                 for (j = 1; j <= i__1; ++j) {
00193                     i__2 = j;
00194                     if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
00195                         i__2 = j;
00196                         temp.r = x[i__2].r, temp.i = x[i__2].i;
00197                         k = kk;
00198                         i__2 = j - 1;
00199                         for (i__ = 1; i__ <= i__2; ++i__) {
00200                             i__3 = i__;
00201                             i__4 = i__;
00202                             i__5 = k;
00203                             q__2.r = temp.r * ap[i__5].r - temp.i * ap[i__5]
00204                                     .i, q__2.i = temp.r * ap[i__5].i + temp.i 
00205                                     * ap[i__5].r;
00206                             q__1.r = x[i__4].r + q__2.r, q__1.i = x[i__4].i + 
00207                                     q__2.i;
00208                             x[i__3].r = q__1.r, x[i__3].i = q__1.i;
00209                             ++k;
00210 /* L10: */
00211                         }
00212                         if (nounit) {
00213                             i__2 = j;
00214                             i__3 = j;
00215                             i__4 = kk + j - 1;
00216                             q__1.r = x[i__3].r * ap[i__4].r - x[i__3].i * ap[
00217                                     i__4].i, q__1.i = x[i__3].r * ap[i__4].i 
00218                                     + x[i__3].i * ap[i__4].r;
00219                             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00220                         }
00221                     }
00222                     kk += j;
00223 /* L20: */
00224                 }
00225             } else {
00226                 jx = kx;
00227                 i__1 = *n;
00228                 for (j = 1; j <= i__1; ++j) {
00229                     i__2 = jx;
00230                     if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
00231                         i__2 = jx;
00232                         temp.r = x[i__2].r, temp.i = x[i__2].i;
00233                         ix = kx;
00234                         i__2 = kk + j - 2;
00235                         for (k = kk; k <= i__2; ++k) {
00236                             i__3 = ix;
00237                             i__4 = ix;
00238                             i__5 = k;
00239                             q__2.r = temp.r * ap[i__5].r - temp.i * ap[i__5]
00240                                     .i, q__2.i = temp.r * ap[i__5].i + temp.i 
00241                                     * ap[i__5].r;
00242                             q__1.r = x[i__4].r + q__2.r, q__1.i = x[i__4].i + 
00243                                     q__2.i;
00244                             x[i__3].r = q__1.r, x[i__3].i = q__1.i;
00245                             ix += *incx;
00246 /* L30: */
00247                         }
00248                         if (nounit) {
00249                             i__2 = jx;
00250                             i__3 = jx;
00251                             i__4 = kk + j - 1;
00252                             q__1.r = x[i__3].r * ap[i__4].r - x[i__3].i * ap[
00253                                     i__4].i, q__1.i = x[i__3].r * ap[i__4].i 
00254                                     + x[i__3].i * ap[i__4].r;
00255                             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00256                         }
00257                     }
00258                     jx += *incx;
00259                     kk += j;
00260 /* L40: */
00261                 }
00262             }
00263         } else {
00264             kk = *n * (*n + 1) / 2;
00265             if (*incx == 1) {
00266                 for (j = *n; j >= 1; --j) {
00267                     i__1 = j;
00268                     if (x[i__1].r != 0.f || x[i__1].i != 0.f) {
00269                         i__1 = j;
00270                         temp.r = x[i__1].r, temp.i = x[i__1].i;
00271                         k = kk;
00272                         i__1 = j + 1;
00273                         for (i__ = *n; i__ >= i__1; --i__) {
00274                             i__2 = i__;
00275                             i__3 = i__;
00276                             i__4 = k;
00277                             q__2.r = temp.r * ap[i__4].r - temp.i * ap[i__4]
00278                                     .i, q__2.i = temp.r * ap[i__4].i + temp.i 
00279                                     * ap[i__4].r;
00280                             q__1.r = x[i__3].r + q__2.r, q__1.i = x[i__3].i + 
00281                                     q__2.i;
00282                             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00283                             --k;
00284 /* L50: */
00285                         }
00286                         if (nounit) {
00287                             i__1 = j;
00288                             i__2 = j;
00289                             i__3 = kk - *n + j;
00290                             q__1.r = x[i__2].r * ap[i__3].r - x[i__2].i * ap[
00291                                     i__3].i, q__1.i = x[i__2].r * ap[i__3].i 
00292                                     + x[i__2].i * ap[i__3].r;
00293                             x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00294                         }
00295                     }
00296                     kk -= *n - j + 1;
00297 /* L60: */
00298                 }
00299             } else {
00300                 kx += (*n - 1) * *incx;
00301                 jx = kx;
00302                 for (j = *n; j >= 1; --j) {
00303                     i__1 = jx;
00304                     if (x[i__1].r != 0.f || x[i__1].i != 0.f) {
00305                         i__1 = jx;
00306                         temp.r = x[i__1].r, temp.i = x[i__1].i;
00307                         ix = kx;
00308                         i__1 = kk - (*n - (j + 1));
00309                         for (k = kk; k >= i__1; --k) {
00310                             i__2 = ix;
00311                             i__3 = ix;
00312                             i__4 = k;
00313                             q__2.r = temp.r * ap[i__4].r - temp.i * ap[i__4]
00314                                     .i, q__2.i = temp.r * ap[i__4].i + temp.i 
00315                                     * ap[i__4].r;
00316                             q__1.r = x[i__3].r + q__2.r, q__1.i = x[i__3].i + 
00317                                     q__2.i;
00318                             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00319                             ix -= *incx;
00320 /* L70: */
00321                         }
00322                         if (nounit) {
00323                             i__1 = jx;
00324                             i__2 = jx;
00325                             i__3 = kk - *n + j;
00326                             q__1.r = x[i__2].r * ap[i__3].r - x[i__2].i * ap[
00327                                     i__3].i, q__1.i = x[i__2].r * ap[i__3].i 
00328                                     + x[i__2].i * ap[i__3].r;
00329                             x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00330                         }
00331                     }
00332                     jx -= *incx;
00333                     kk -= *n - j + 1;
00334 /* L80: */
00335                 }
00336             }
00337         }
00338     } else {
00339 
00340 /*        Form  x := A'*x  or  x := conjg( A' )*x. */
00341 
00342         if (lsame_(uplo, "U")) {
00343             kk = *n * (*n + 1) / 2;
00344             if (*incx == 1) {
00345                 for (j = *n; j >= 1; --j) {
00346                     i__1 = j;
00347                     temp.r = x[i__1].r, temp.i = x[i__1].i;
00348                     k = kk - 1;
00349                     if (noconj) {
00350                         if (nounit) {
00351                             i__1 = kk;
00352                             q__1.r = temp.r * ap[i__1].r - temp.i * ap[i__1]
00353                                     .i, q__1.i = temp.r * ap[i__1].i + temp.i 
00354                                     * ap[i__1].r;
00355                             temp.r = q__1.r, temp.i = q__1.i;
00356                         }
00357                         for (i__ = j - 1; i__ >= 1; --i__) {
00358                             i__1 = k;
00359                             i__2 = i__;
00360                             q__2.r = ap[i__1].r * x[i__2].r - ap[i__1].i * x[
00361                                     i__2].i, q__2.i = ap[i__1].r * x[i__2].i 
00362                                     + ap[i__1].i * x[i__2].r;
00363                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00364                                     q__2.i;
00365                             temp.r = q__1.r, temp.i = q__1.i;
00366                             --k;
00367 /* L90: */
00368                         }
00369                     } else {
00370                         if (nounit) {
00371                             r_cnjg(&q__2, &ap[kk]);
00372                             q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
00373                                     q__1.i = temp.r * q__2.i + temp.i * 
00374                                     q__2.r;
00375                             temp.r = q__1.r, temp.i = q__1.i;
00376                         }
00377                         for (i__ = j - 1; i__ >= 1; --i__) {
00378                             r_cnjg(&q__3, &ap[k]);
00379                             i__1 = i__;
00380                             q__2.r = q__3.r * x[i__1].r - q__3.i * x[i__1].i, 
00381                                     q__2.i = q__3.r * x[i__1].i + q__3.i * x[
00382                                     i__1].r;
00383                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00384                                     q__2.i;
00385                             temp.r = q__1.r, temp.i = q__1.i;
00386                             --k;
00387 /* L100: */
00388                         }
00389                     }
00390                     i__1 = j;
00391                     x[i__1].r = temp.r, x[i__1].i = temp.i;
00392                     kk -= j;
00393 /* L110: */
00394                 }
00395             } else {
00396                 jx = kx + (*n - 1) * *incx;
00397                 for (j = *n; j >= 1; --j) {
00398                     i__1 = jx;
00399                     temp.r = x[i__1].r, temp.i = x[i__1].i;
00400                     ix = jx;
00401                     if (noconj) {
00402                         if (nounit) {
00403                             i__1 = kk;
00404                             q__1.r = temp.r * ap[i__1].r - temp.i * ap[i__1]
00405                                     .i, q__1.i = temp.r * ap[i__1].i + temp.i 
00406                                     * ap[i__1].r;
00407                             temp.r = q__1.r, temp.i = q__1.i;
00408                         }
00409                         i__1 = kk - j + 1;
00410                         for (k = kk - 1; k >= i__1; --k) {
00411                             ix -= *incx;
00412                             i__2 = k;
00413                             i__3 = ix;
00414                             q__2.r = ap[i__2].r * x[i__3].r - ap[i__2].i * x[
00415                                     i__3].i, q__2.i = ap[i__2].r * x[i__3].i 
00416                                     + ap[i__2].i * x[i__3].r;
00417                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00418                                     q__2.i;
00419                             temp.r = q__1.r, temp.i = q__1.i;
00420 /* L120: */
00421                         }
00422                     } else {
00423                         if (nounit) {
00424                             r_cnjg(&q__2, &ap[kk]);
00425                             q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
00426                                     q__1.i = temp.r * q__2.i + temp.i * 
00427                                     q__2.r;
00428                             temp.r = q__1.r, temp.i = q__1.i;
00429                         }
00430                         i__1 = kk - j + 1;
00431                         for (k = kk - 1; k >= i__1; --k) {
00432                             ix -= *incx;
00433                             r_cnjg(&q__3, &ap[k]);
00434                             i__2 = ix;
00435                             q__2.r = q__3.r * x[i__2].r - q__3.i * x[i__2].i, 
00436                                     q__2.i = q__3.r * x[i__2].i + q__3.i * x[
00437                                     i__2].r;
00438                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00439                                     q__2.i;
00440                             temp.r = q__1.r, temp.i = q__1.i;
00441 /* L130: */
00442                         }
00443                     }
00444                     i__1 = jx;
00445                     x[i__1].r = temp.r, x[i__1].i = temp.i;
00446                     jx -= *incx;
00447                     kk -= j;
00448 /* L140: */
00449                 }
00450             }
00451         } else {
00452             kk = 1;
00453             if (*incx == 1) {
00454                 i__1 = *n;
00455                 for (j = 1; j <= i__1; ++j) {
00456                     i__2 = j;
00457                     temp.r = x[i__2].r, temp.i = x[i__2].i;
00458                     k = kk + 1;
00459                     if (noconj) {
00460                         if (nounit) {
00461                             i__2 = kk;
00462                             q__1.r = temp.r * ap[i__2].r - temp.i * ap[i__2]
00463                                     .i, q__1.i = temp.r * ap[i__2].i + temp.i 
00464                                     * ap[i__2].r;
00465                             temp.r = q__1.r, temp.i = q__1.i;
00466                         }
00467                         i__2 = *n;
00468                         for (i__ = j + 1; i__ <= i__2; ++i__) {
00469                             i__3 = k;
00470                             i__4 = i__;
00471                             q__2.r = ap[i__3].r * x[i__4].r - ap[i__3].i * x[
00472                                     i__4].i, q__2.i = ap[i__3].r * x[i__4].i 
00473                                     + ap[i__3].i * x[i__4].r;
00474                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00475                                     q__2.i;
00476                             temp.r = q__1.r, temp.i = q__1.i;
00477                             ++k;
00478 /* L150: */
00479                         }
00480                     } else {
00481                         if (nounit) {
00482                             r_cnjg(&q__2, &ap[kk]);
00483                             q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
00484                                     q__1.i = temp.r * q__2.i + temp.i * 
00485                                     q__2.r;
00486                             temp.r = q__1.r, temp.i = q__1.i;
00487                         }
00488                         i__2 = *n;
00489                         for (i__ = j + 1; i__ <= i__2; ++i__) {
00490                             r_cnjg(&q__3, &ap[k]);
00491                             i__3 = i__;
00492                             q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, 
00493                                     q__2.i = q__3.r * x[i__3].i + q__3.i * x[
00494                                     i__3].r;
00495                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00496                                     q__2.i;
00497                             temp.r = q__1.r, temp.i = q__1.i;
00498                             ++k;
00499 /* L160: */
00500                         }
00501                     }
00502                     i__2 = j;
00503                     x[i__2].r = temp.r, x[i__2].i = temp.i;
00504                     kk += *n - j + 1;
00505 /* L170: */
00506                 }
00507             } else {
00508                 jx = kx;
00509                 i__1 = *n;
00510                 for (j = 1; j <= i__1; ++j) {
00511                     i__2 = jx;
00512                     temp.r = x[i__2].r, temp.i = x[i__2].i;
00513                     ix = jx;
00514                     if (noconj) {
00515                         if (nounit) {
00516                             i__2 = kk;
00517                             q__1.r = temp.r * ap[i__2].r - temp.i * ap[i__2]
00518                                     .i, q__1.i = temp.r * ap[i__2].i + temp.i 
00519                                     * ap[i__2].r;
00520                             temp.r = q__1.r, temp.i = q__1.i;
00521                         }
00522                         i__2 = kk + *n - j;
00523                         for (k = kk + 1; k <= i__2; ++k) {
00524                             ix += *incx;
00525                             i__3 = k;
00526                             i__4 = ix;
00527                             q__2.r = ap[i__3].r * x[i__4].r - ap[i__3].i * x[
00528                                     i__4].i, q__2.i = ap[i__3].r * x[i__4].i 
00529                                     + ap[i__3].i * x[i__4].r;
00530                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00531                                     q__2.i;
00532                             temp.r = q__1.r, temp.i = q__1.i;
00533 /* L180: */
00534                         }
00535                     } else {
00536                         if (nounit) {
00537                             r_cnjg(&q__2, &ap[kk]);
00538                             q__1.r = temp.r * q__2.r - temp.i * q__2.i, 
00539                                     q__1.i = temp.r * q__2.i + temp.i * 
00540                                     q__2.r;
00541                             temp.r = q__1.r, temp.i = q__1.i;
00542                         }
00543                         i__2 = kk + *n - j;
00544                         for (k = kk + 1; k <= i__2; ++k) {
00545                             ix += *incx;
00546                             r_cnjg(&q__3, &ap[k]);
00547                             i__3 = ix;
00548                             q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, 
00549                                     q__2.i = q__3.r * x[i__3].i + q__3.i * x[
00550                                     i__3].r;
00551                             q__1.r = temp.r + q__2.r, q__1.i = temp.i + 
00552                                     q__2.i;
00553                             temp.r = q__1.r, temp.i = q__1.i;
00554 /* L190: */
00555                         }
00556                     }
00557                     i__2 = jx;
00558                     x[i__2].r = temp.r, x[i__2].i = temp.i;
00559                     jx += *incx;
00560                     kk += *n - j + 1;
00561 /* L200: */
00562                 }
00563             }
00564         }
00565     }
00566 
00567     return 0;
00568 
00569 /*     End of CTPMV . */
00570 
00571 } /* ctpmv_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:34