sgtts2.c
Go to the documentation of this file.
00001 /* sgtts2.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 sgtts2_(integer *itrans, integer *n, integer *nrhs, real 
00017         *dl, real *d__, real *du, real *du2, integer *ipiv, real *b, integer *
00018         ldb)
00019 {
00020     /* System generated locals */
00021     integer b_dim1, b_offset, i__1, i__2;
00022 
00023     /* Local variables */
00024     integer i__, j, ip;
00025     real temp;
00026 
00027 
00028 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00029 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00030 /*     November 2006 */
00031 
00032 /*     .. Scalar Arguments .. */
00033 /*     .. */
00034 /*     .. Array Arguments .. */
00035 /*     .. */
00036 
00037 /*  Purpose */
00038 /*  ======= */
00039 
00040 /*  SGTTS2 solves one of the systems of equations */
00041 /*     A*X = B  or  A'*X = B, */
00042 /*  with a tridiagonal matrix A using the LU factorization computed */
00043 /*  by SGTTRF. */
00044 
00045 /*  Arguments */
00046 /*  ========= */
00047 
00048 /*  ITRANS  (input) INTEGER */
00049 /*          Specifies the form of the system of equations. */
00050 /*          = 0:  A * X = B  (No transpose) */
00051 /*          = 1:  A'* X = B  (Transpose) */
00052 /*          = 2:  A'* X = B  (Conjugate transpose = Transpose) */
00053 
00054 /*  N       (input) INTEGER */
00055 /*          The order of the matrix A. */
00056 
00057 /*  NRHS    (input) INTEGER */
00058 /*          The number of right hand sides, i.e., the number of columns */
00059 /*          of the matrix B.  NRHS >= 0. */
00060 
00061 /*  DL      (input) REAL array, dimension (N-1) */
00062 /*          The (n-1) multipliers that define the matrix L from the */
00063 /*          LU factorization of A. */
00064 
00065 /*  D       (input) REAL array, dimension (N) */
00066 /*          The n diagonal elements of the upper triangular matrix U from */
00067 /*          the LU factorization of A. */
00068 
00069 /*  DU      (input) REAL array, dimension (N-1) */
00070 /*          The (n-1) elements of the first super-diagonal of U. */
00071 
00072 /*  DU2     (input) REAL array, dimension (N-2) */
00073 /*          The (n-2) elements of the second super-diagonal of U. */
00074 
00075 /*  IPIV    (input) INTEGER array, dimension (N) */
00076 /*          The pivot indices; for 1 <= i <= n, row i of the matrix was */
00077 /*          interchanged with row IPIV(i).  IPIV(i) will always be either */
00078 /*          i or i+1; IPIV(i) = i indicates a row interchange was not */
00079 /*          required. */
00080 
00081 /*  B       (input/output) REAL array, dimension (LDB,NRHS) */
00082 /*          On entry, the matrix of right hand side vectors B. */
00083 /*          On exit, B is overwritten by the solution vectors X. */
00084 
00085 /*  LDB     (input) INTEGER */
00086 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00087 
00088 /*  ===================================================================== */
00089 
00090 /*     .. Local Scalars .. */
00091 /*     .. */
00092 /*     .. Executable Statements .. */
00093 
00094 /*     Quick return if possible */
00095 
00096     /* Parameter adjustments */
00097     --dl;
00098     --d__;
00099     --du;
00100     --du2;
00101     --ipiv;
00102     b_dim1 = *ldb;
00103     b_offset = 1 + b_dim1;
00104     b -= b_offset;
00105 
00106     /* Function Body */
00107     if (*n == 0 || *nrhs == 0) {
00108         return 0;
00109     }
00110 
00111     if (*itrans == 0) {
00112 
00113 /*        Solve A*X = B using the LU factorization of A, */
00114 /*        overwriting each right hand side vector with its solution. */
00115 
00116         if (*nrhs <= 1) {
00117             j = 1;
00118 L10:
00119 
00120 /*           Solve L*x = b. */
00121 
00122             i__1 = *n - 1;
00123             for (i__ = 1; i__ <= i__1; ++i__) {
00124                 ip = ipiv[i__];
00125                 temp = b[i__ + 1 - ip + i__ + j * b_dim1] - dl[i__] * b[ip + 
00126                         j * b_dim1];
00127                 b[i__ + j * b_dim1] = b[ip + j * b_dim1];
00128                 b[i__ + 1 + j * b_dim1] = temp;
00129 /* L20: */
00130             }
00131 
00132 /*           Solve U*x = b. */
00133 
00134             b[*n + j * b_dim1] /= d__[*n];
00135             if (*n > 1) {
00136                 b[*n - 1 + j * b_dim1] = (b[*n - 1 + j * b_dim1] - du[*n - 1] 
00137                         * b[*n + j * b_dim1]) / d__[*n - 1];
00138             }
00139             for (i__ = *n - 2; i__ >= 1; --i__) {
00140                 b[i__ + j * b_dim1] = (b[i__ + j * b_dim1] - du[i__] * b[i__ 
00141                         + 1 + j * b_dim1] - du2[i__] * b[i__ + 2 + j * b_dim1]
00142                         ) / d__[i__];
00143 /* L30: */
00144             }
00145             if (j < *nrhs) {
00146                 ++j;
00147                 goto L10;
00148             }
00149         } else {
00150             i__1 = *nrhs;
00151             for (j = 1; j <= i__1; ++j) {
00152 
00153 /*              Solve L*x = b. */
00154 
00155                 i__2 = *n - 1;
00156                 for (i__ = 1; i__ <= i__2; ++i__) {
00157                     if (ipiv[i__] == i__) {
00158                         b[i__ + 1 + j * b_dim1] -= dl[i__] * b[i__ + j * 
00159                                 b_dim1];
00160                     } else {
00161                         temp = b[i__ + j * b_dim1];
00162                         b[i__ + j * b_dim1] = b[i__ + 1 + j * b_dim1];
00163                         b[i__ + 1 + j * b_dim1] = temp - dl[i__] * b[i__ + j *
00164                                  b_dim1];
00165                     }
00166 /* L40: */
00167                 }
00168 
00169 /*              Solve U*x = b. */
00170 
00171                 b[*n + j * b_dim1] /= d__[*n];
00172                 if (*n > 1) {
00173                     b[*n - 1 + j * b_dim1] = (b[*n - 1 + j * b_dim1] - du[*n 
00174                             - 1] * b[*n + j * b_dim1]) / d__[*n - 1];
00175                 }
00176                 for (i__ = *n - 2; i__ >= 1; --i__) {
00177                     b[i__ + j * b_dim1] = (b[i__ + j * b_dim1] - du[i__] * b[
00178                             i__ + 1 + j * b_dim1] - du2[i__] * b[i__ + 2 + j *
00179                              b_dim1]) / d__[i__];
00180 /* L50: */
00181                 }
00182 /* L60: */
00183             }
00184         }
00185     } else {
00186 
00187 /*        Solve A' * X = B. */
00188 
00189         if (*nrhs <= 1) {
00190 
00191 /*           Solve U'*x = b. */
00192 
00193             j = 1;
00194 L70:
00195             b[j * b_dim1 + 1] /= d__[1];
00196             if (*n > 1) {
00197                 b[j * b_dim1 + 2] = (b[j * b_dim1 + 2] - du[1] * b[j * b_dim1 
00198                         + 1]) / d__[2];
00199             }
00200             i__1 = *n;
00201             for (i__ = 3; i__ <= i__1; ++i__) {
00202                 b[i__ + j * b_dim1] = (b[i__ + j * b_dim1] - du[i__ - 1] * b[
00203                         i__ - 1 + j * b_dim1] - du2[i__ - 2] * b[i__ - 2 + j *
00204                          b_dim1]) / d__[i__];
00205 /* L80: */
00206             }
00207 
00208 /*           Solve L'*x = b. */
00209 
00210             for (i__ = *n - 1; i__ >= 1; --i__) {
00211                 ip = ipiv[i__];
00212                 temp = b[i__ + j * b_dim1] - dl[i__] * b[i__ + 1 + j * b_dim1]
00213                         ;
00214                 b[i__ + j * b_dim1] = b[ip + j * b_dim1];
00215                 b[ip + j * b_dim1] = temp;
00216 /* L90: */
00217             }
00218             if (j < *nrhs) {
00219                 ++j;
00220                 goto L70;
00221             }
00222 
00223         } else {
00224             i__1 = *nrhs;
00225             for (j = 1; j <= i__1; ++j) {
00226 
00227 /*              Solve U'*x = b. */
00228 
00229                 b[j * b_dim1 + 1] /= d__[1];
00230                 if (*n > 1) {
00231                     b[j * b_dim1 + 2] = (b[j * b_dim1 + 2] - du[1] * b[j * 
00232                             b_dim1 + 1]) / d__[2];
00233                 }
00234                 i__2 = *n;
00235                 for (i__ = 3; i__ <= i__2; ++i__) {
00236                     b[i__ + j * b_dim1] = (b[i__ + j * b_dim1] - du[i__ - 1] *
00237                              b[i__ - 1 + j * b_dim1] - du2[i__ - 2] * b[i__ - 
00238                             2 + j * b_dim1]) / d__[i__];
00239 /* L100: */
00240                 }
00241                 for (i__ = *n - 1; i__ >= 1; --i__) {
00242                     if (ipiv[i__] == i__) {
00243                         b[i__ + j * b_dim1] -= dl[i__] * b[i__ + 1 + j * 
00244                                 b_dim1];
00245                     } else {
00246                         temp = b[i__ + 1 + j * b_dim1];
00247                         b[i__ + 1 + j * b_dim1] = b[i__ + j * b_dim1] - dl[
00248                                 i__] * temp;
00249                         b[i__ + j * b_dim1] = temp;
00250                     }
00251 /* L110: */
00252                 }
00253 /* L120: */
00254             }
00255         }
00256     }
00257 
00258 /*     End of SGTTS2 */
00259 
00260     return 0;
00261 } /* sgtts2_ */


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