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

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