Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages | Examples

dtrsv.c

Go to the documentation of this file.
00001 /* DTRSV.F -- translated by f2c (version 19941215).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 
00009 /* *********************************************************************** */
00010 
00011 /* Subroutine */ int dtrsv_(uplo, trans, diag, n, a, lda, x, incx, uplo_len, 
00012         trans_len, diag_len)
00013 char *uplo, *trans, *diag;
00014 integer *n;
00015 doublereal *a;
00016 integer *lda;
00017 doublereal *x;
00018 integer *incx;
00019 ftnlen uplo_len;
00020 ftnlen trans_len;
00021 ftnlen diag_len;
00022 {
00023     /* System generated locals */
00024     integer a_dim1, a_offset, i__1, i__2;
00025 
00026     /* Local variables */
00027     static integer info;
00028     static doublereal temp;
00029     static integer i, j;
00030     extern logical lsame_();
00031     static integer ix, jx, kx;
00032     extern /* Subroutine */ int xerbla_();
00033     static logical nounit;
00034 
00035 /*     .. Scalar Arguments .. */
00036 /*     .. Array Arguments .. */
00037 /*     .. */
00038 
00039 /*  Purpose */
00040 /*  ======= */
00041 
00042 /*  DTRSV  solves one of the systems of equations */
00043 
00044 /*     A*x = b,   or   A'*x = b, */
00045 
00046 /*  where b and x are n element vectors and A is an n by n unit, or */
00047 /*  non-unit, upper or lower triangular matrix. */
00048 
00049 /*  No test for singularity or near-singularity is included in this */
00050 /*  routine. Such tests must be performed before calling this routine. */
00051 
00052 /*  Parameters */
00053 /*  ========== */
00054 
00055 /*  UPLO   - CHARACTER*1. */
00056 /*           On entry, UPLO specifies whether the matrix is an upper or */
00057 /*           lower triangular matrix as follows: */
00058 
00059 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
00060 
00061 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
00062 
00063 /*           Unchanged on exit. */
00064 
00065 /*  TRANS  - CHARACTER*1. */
00066 /*           On entry, TRANS specifies the equations to be solved as */
00067 /*           follows: */
00068 
00069 /*              TRANS = 'N' or 'n'   A*x = b. */
00070 
00071 /*              TRANS = 'T' or 't'   A'*x = b. */
00072 
00073 /*              TRANS = 'C' or 'c'   A'*x = b. */
00074 
00075 /*           Unchanged on exit. */
00076 
00077 /*  DIAG   - CHARACTER*1. */
00078 /*           On entry, DIAG specifies whether or not A is unit */
00079 /*           triangular as follows: */
00080 
00081 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00082 
00083 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00084 /*                                  triangular. */
00085 
00086 /*           Unchanged on exit. */
00087 
00088 /*  N      - INTEGER. */
00089 /*           On entry, N specifies the order of the matrix A. */
00090 /*           N must be at least zero. */
00091 /*           Unchanged on exit. */
00092 
00093 /*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
00094 /*           Before entry with  UPLO = 'U' or 'u', the leading n by n */
00095 /*           upper triangular part of the array A must contain the upper 
00096 */
00097 /*           triangular matrix and the strictly lower triangular part of 
00098 */
00099 /*           A is not referenced. */
00100 /*           Before entry with UPLO = 'L' or 'l', the leading n by n */
00101 /*           lower triangular part of the array A must contain the lower 
00102 */
00103 /*           triangular matrix and the strictly upper triangular part of 
00104 */
00105 /*           A is not referenced. */
00106 /*           Note that when  DIAG = 'U' or 'u', the diagonal elements of 
00107 */
00108 /*           A are not referenced either, but are assumed to be unity. */
00109 /*           Unchanged on exit. */
00110 
00111 /*  LDA    - INTEGER. */
00112 /*           On entry, LDA specifies the first dimension of A as declared 
00113 */
00114 /*           in the calling (sub) program. LDA must be at least */
00115 /*           max( 1, n ). */
00116 /*           Unchanged on exit. */
00117 
00118 /*  X      - DOUBLE PRECISION array of dimension at least */
00119 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00120 /*           Before entry, the incremented array X must contain the n */
00121 /*           element right-hand side vector b. On exit, X is overwritten 
00122 */
00123 /*           with the solution vector x. */
00124 
00125 /*  INCX   - INTEGER. */
00126 /*           On entry, INCX specifies the increment for the elements of */
00127 /*           X. INCX must not be zero. */
00128 /*           Unchanged on exit. */
00129 
00130 
00131 /*  Level 2 Blas routine. */
00132 
00133 /*  -- Written on 22-October-1986. */
00134 /*     Jack Dongarra, Argonne National Lab. */
00135 /*     Jeremy Du Croz, Nag Central Office. */
00136 /*     Sven Hammarling, Nag Central Office. */
00137 /*     Richard Hanson, Sandia National Labs. */
00138 
00139 
00140 /*     .. Parameters .. */
00141 /*     .. Local Scalars .. */
00142 /*     .. External Functions .. */
00143 /*     .. External Subroutines .. */
00144 /*     .. Intrinsic Functions .. */
00145 /*     .. */
00146 /*     .. Executable Statements .. */
00147 
00148 /*     Test the input parameters. */
00149 
00150     /* Parameter adjustments */
00151     a_dim1 = *lda;
00152     a_offset = a_dim1 + 1;
00153     a -= a_offset;
00154     --x;
00155 
00156     /* Function Body */
00157     info = 0;
00158     if (! lsame_(uplo, "U", 1L, 1L) && ! lsame_(uplo, "L", 1L, 1L)) {
00159         info = 1;
00160     } else if (! lsame_(trans, "N", 1L, 1L) && ! lsame_(trans, "T", 1L, 1L) &&
00161              ! lsame_(trans, "C", 1L, 1L)) {
00162         info = 2;
00163     } else if (! lsame_(diag, "U", 1L, 1L) && ! lsame_(diag, "N", 1L, 1L)) {
00164         info = 3;
00165     } else if (*n < 0) {
00166         info = 4;
00167     } else if (*lda < max(1,*n)) {
00168         info = 6;
00169     } else if (*incx == 0) {
00170         info = 8;
00171     }
00172     if (info != 0) {
00173         xerbla_("DTRSV ", &info, 6L);
00174         return 0;
00175     }
00176 
00177 /*     Quick return if possible. */
00178 
00179     if (*n == 0) {
00180         return 0;
00181     }
00182 
00183     nounit = lsame_(diag, "N", 1L, 1L);
00184 
00185 /*     Set up the start point in X if the increment is not unity. This */
00186 /*     will be  ( N - 1 )*INCX  too small for descending loops. */
00187 
00188     if (*incx <= 0) {
00189         kx = 1 - (*n - 1) * *incx;
00190     } else if (*incx != 1) {
00191         kx = 1;
00192     }
00193 
00194 /*     Start the operations. In this version the elements of A are */
00195 /*     accessed sequentially with one pass through A. */
00196 
00197     if (lsame_(trans, "N", 1L, 1L)) {
00198 
00199 /*        Form  x := inv( A )*x. */
00200 
00201         if (lsame_(uplo, "U", 1L, 1L)) {
00202             if (*incx == 1) {
00203                 for (j = *n; j >= 1; --j) {
00204                     if (x[j] != 0.) {
00205                         if (nounit) {
00206                             x[j] /= a[j + j * a_dim1];
00207                         }
00208                         temp = x[j];
00209                         for (i = j - 1; i >= 1; --i) {
00210                             x[i] -= temp * a[i + j * a_dim1];
00211 /* L10: */
00212                         }
00213                     }
00214 /* L20: */
00215                 }
00216             } else {
00217                 jx = kx + (*n - 1) * *incx;
00218                 for (j = *n; j >= 1; --j) {
00219                     if (x[jx] != 0.) {
00220                         if (nounit) {
00221                             x[jx] /= a[j + j * a_dim1];
00222                         }
00223                         temp = x[jx];
00224                         ix = jx;
00225                         for (i = j - 1; i >= 1; --i) {
00226                             ix -= *incx;
00227                             x[ix] -= temp * a[i + j * a_dim1];
00228 /* L30: */
00229                         }
00230                     }
00231                     jx -= *incx;
00232 /* L40: */
00233                 }
00234             }
00235         } else {
00236             if (*incx == 1) {
00237                 i__1 = *n;
00238                 for (j = 1; j <= i__1; ++j) {
00239                     if (x[j] != 0.) {
00240                         if (nounit) {
00241                             x[j] /= a[j + j * a_dim1];
00242                         }
00243                         temp = x[j];
00244                         i__2 = *n;
00245                         for (i = j + 1; i <= i__2; ++i) {
00246                             x[i] -= temp * a[i + j * a_dim1];
00247 /* L50: */
00248                         }
00249                     }
00250 /* L60: */
00251                 }
00252             } else {
00253                 jx = kx;
00254                 i__1 = *n;
00255                 for (j = 1; j <= i__1; ++j) {
00256                     if (x[jx] != 0.) {
00257                         if (nounit) {
00258                             x[jx] /= a[j + j * a_dim1];
00259                         }
00260                         temp = x[jx];
00261                         ix = jx;
00262                         i__2 = *n;
00263                         for (i = j + 1; i <= i__2; ++i) {
00264                             ix += *incx;
00265                             x[ix] -= temp * a[i + j * a_dim1];
00266 /* L70: */
00267                         }
00268                     }
00269                     jx += *incx;
00270 /* L80: */
00271                 }
00272             }
00273         }
00274     } else {
00275 
00276 /*        Form  x := inv( A' )*x. */
00277 
00278         if (lsame_(uplo, "U", 1L, 1L)) {
00279             if (*incx == 1) {
00280                 i__1 = *n;
00281                 for (j = 1; j <= i__1; ++j) {
00282                     temp = x[j];
00283                     i__2 = j - 1;
00284                     for (i = 1; i <= i__2; ++i) {
00285                         temp -= a[i + j * a_dim1] * x[i];
00286 /* L90: */
00287                     }
00288                     if (nounit) {
00289                         temp /= a[j + j * a_dim1];
00290                     }
00291                     x[j] = temp;
00292 /* L100: */
00293                 }
00294             } else {
00295                 jx = kx;
00296                 i__1 = *n;
00297                 for (j = 1; j <= i__1; ++j) {
00298                     temp = x[jx];
00299                     ix = kx;
00300                     i__2 = j - 1;
00301                     for (i = 1; i <= i__2; ++i) {
00302                         temp -= a[i + j * a_dim1] * x[ix];
00303                         ix += *incx;
00304 /* L110: */
00305                     }
00306                     if (nounit) {
00307                         temp /= a[j + j * a_dim1];
00308                     }
00309                     x[jx] = temp;
00310                     jx += *incx;
00311 /* L120: */
00312                 }
00313             }
00314         } else {
00315             if (*incx == 1) {
00316                 for (j = *n; j >= 1; --j) {
00317                     temp = x[j];
00318                     i__1 = j + 1;
00319                     for (i = *n; i >= i__1; --i) {
00320                         temp -= a[i + j * a_dim1] * x[i];
00321 /* L130: */
00322                     }
00323                     if (nounit) {
00324                         temp /= a[j + j * a_dim1];
00325                     }
00326                     x[j] = temp;
00327 /* L140: */
00328                 }
00329             } else {
00330                 kx += (*n - 1) * *incx;
00331                 jx = kx;
00332                 for (j = *n; j >= 1; --j) {
00333                     temp = x[jx];
00334                     ix = kx;
00335                     i__1 = j + 1;
00336                     for (i = *n; i >= i__1; --i) {
00337                         temp -= a[i + j * a_dim1] * x[ix];
00338                         ix -= *incx;
00339 /* L150: */
00340                     }
00341                     if (nounit) {
00342                         temp /= a[j + j * a_dim1];
00343                     }
00344                     x[jx] = temp;
00345                     jx -= *incx;
00346 /* L160: */
00347                 }
00348             }
00349         }
00350     }
00351 
00352     return 0;
00353 
00354 /*     End of DTRSV . */
00355 
00356 } /* dtrsv_ */
00357 

Generated on Wed Sep 5 12:54:20 2007 for DSACSS Operational Code by  doxygen 1.3.9.1