dgttrf.c
Go to the documentation of this file.
00001 /* dgttrf.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 dgttrf_(integer *n, doublereal *dl, doublereal *d__, 
00017         doublereal *du, doublereal *du2, integer *ipiv, integer *info)
00018 {
00019     /* System generated locals */
00020     integer i__1;
00021     doublereal d__1, d__2;
00022 
00023     /* Local variables */
00024     integer i__;
00025     doublereal fact, temp;
00026     extern /* Subroutine */ int xerbla_(char *, integer *);
00027 
00028 
00029 /*  -- LAPACK routine (version 3.2) -- */
00030 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00031 /*     November 2006 */
00032 
00033 /*     .. Scalar Arguments .. */
00034 /*     .. */
00035 /*     .. Array Arguments .. */
00036 /*     .. */
00037 
00038 /*  Purpose */
00039 /*  ======= */
00040 
00041 /*  DGTTRF computes an LU factorization of a real tridiagonal matrix A */
00042 /*  using elimination with partial pivoting and row interchanges. */
00043 
00044 /*  The factorization has the form */
00045 /*     A = L * U */
00046 /*  where L is a product of permutation and unit lower bidiagonal */
00047 /*  matrices and U is upper triangular with nonzeros in only the main */
00048 /*  diagonal and first two superdiagonals. */
00049 
00050 /*  Arguments */
00051 /*  ========= */
00052 
00053 /*  N       (input) INTEGER */
00054 /*          The order of the matrix A. */
00055 
00056 /*  DL      (input/output) DOUBLE PRECISION array, dimension (N-1) */
00057 /*          On entry, DL must contain the (n-1) sub-diagonal elements of */
00058 /*          A. */
00059 
00060 /*          On exit, DL is overwritten by the (n-1) multipliers that */
00061 /*          define the matrix L from the LU factorization of A. */
00062 
00063 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
00064 /*          On entry, D must contain the diagonal elements of A. */
00065 
00066 /*          On exit, D is overwritten by the n diagonal elements of the */
00067 /*          upper triangular matrix U from the LU factorization of A. */
00068 
00069 /*  DU      (input/output) DOUBLE PRECISION array, dimension (N-1) */
00070 /*          On entry, DU must contain the (n-1) super-diagonal elements */
00071 /*          of A. */
00072 
00073 /*          On exit, DU is overwritten by the (n-1) elements of the first */
00074 /*          super-diagonal of U. */
00075 
00076 /*  DU2     (output) DOUBLE PRECISION array, dimension (N-2) */
00077 /*          On exit, DU2 is overwritten by the (n-2) elements of the */
00078 /*          second super-diagonal of U. */
00079 
00080 /*  IPIV    (output) INTEGER array, dimension (N) */
00081 /*          The pivot indices; for 1 <= i <= n, row i of the matrix was */
00082 /*          interchanged with row IPIV(i).  IPIV(i) will always be either */
00083 /*          i or i+1; IPIV(i) = i indicates a row interchange was not */
00084 /*          required. */
00085 
00086 /*  INFO    (output) INTEGER */
00087 /*          = 0:  successful exit */
00088 /*          < 0:  if INFO = -k, the k-th argument had an illegal value */
00089 /*          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization */
00090 /*                has been completed, but the factor U is exactly */
00091 /*                singular, and division by zero will occur if it is used */
00092 /*                to solve a system of equations. */
00093 
00094 /*  ===================================================================== */
00095 
00096 /*     .. Parameters .. */
00097 /*     .. */
00098 /*     .. Local Scalars .. */
00099 /*     .. */
00100 /*     .. Intrinsic Functions .. */
00101 /*     .. */
00102 /*     .. External Subroutines .. */
00103 /*     .. */
00104 /*     .. Executable Statements .. */
00105 
00106     /* Parameter adjustments */
00107     --ipiv;
00108     --du2;
00109     --du;
00110     --d__;
00111     --dl;
00112 
00113     /* Function Body */
00114     *info = 0;
00115     if (*n < 0) {
00116         *info = -1;
00117         i__1 = -(*info);
00118         xerbla_("DGTTRF", &i__1);
00119         return 0;
00120     }
00121 
00122 /*     Quick return if possible */
00123 
00124     if (*n == 0) {
00125         return 0;
00126     }
00127 
00128 /*     Initialize IPIV(i) = i and DU2(I) = 0 */
00129 
00130     i__1 = *n;
00131     for (i__ = 1; i__ <= i__1; ++i__) {
00132         ipiv[i__] = i__;
00133 /* L10: */
00134     }
00135     i__1 = *n - 2;
00136     for (i__ = 1; i__ <= i__1; ++i__) {
00137         du2[i__] = 0.;
00138 /* L20: */
00139     }
00140 
00141     i__1 = *n - 2;
00142     for (i__ = 1; i__ <= i__1; ++i__) {
00143         if ((d__1 = d__[i__], abs(d__1)) >= (d__2 = dl[i__], abs(d__2))) {
00144 
00145 /*           No row interchange required, eliminate DL(I) */
00146 
00147             if (d__[i__] != 0.) {
00148                 fact = dl[i__] / d__[i__];
00149                 dl[i__] = fact;
00150                 d__[i__ + 1] -= fact * du[i__];
00151             }
00152         } else {
00153 
00154 /*           Interchange rows I and I+1, eliminate DL(I) */
00155 
00156             fact = d__[i__] / dl[i__];
00157             d__[i__] = dl[i__];
00158             dl[i__] = fact;
00159             temp = du[i__];
00160             du[i__] = d__[i__ + 1];
00161             d__[i__ + 1] = temp - fact * d__[i__ + 1];
00162             du2[i__] = du[i__ + 1];
00163             du[i__ + 1] = -fact * du[i__ + 1];
00164             ipiv[i__] = i__ + 1;
00165         }
00166 /* L30: */
00167     }
00168     if (*n > 1) {
00169         i__ = *n - 1;
00170         if ((d__1 = d__[i__], abs(d__1)) >= (d__2 = dl[i__], abs(d__2))) {
00171             if (d__[i__] != 0.) {
00172                 fact = dl[i__] / d__[i__];
00173                 dl[i__] = fact;
00174                 d__[i__ + 1] -= fact * du[i__];
00175             }
00176         } else {
00177             fact = d__[i__] / dl[i__];
00178             d__[i__] = dl[i__];
00179             dl[i__] = fact;
00180             temp = du[i__];
00181             du[i__] = d__[i__ + 1];
00182             d__[i__ + 1] = temp - fact * d__[i__ + 1];
00183             ipiv[i__] = i__ + 1;
00184         }
00185     }
00186 
00187 /*     Check for a zero on the diagonal of U. */
00188 
00189     i__1 = *n;
00190     for (i__ = 1; i__ <= i__1; ++i__) {
00191         if (d__[i__] == 0.) {
00192             *info = i__;
00193             goto L50;
00194         }
00195 /* L40: */
00196     }
00197 L50:
00198 
00199     return 0;
00200 
00201 /*     End of DGTTRF */
00202 
00203 } /* dgttrf_ */


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