001:       SUBROUTINE SGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, N
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       REAL               D( * ), DL( * ), DU( * ), DU2( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  SGTTRF computes an LU factorization of a real tridiagonal matrix A
020: *  using elimination with partial pivoting and row interchanges.
021: *
022: *  The factorization has the form
023: *     A = L * U
024: *  where L is a product of permutation and unit lower bidiagonal
025: *  matrices and U is upper triangular with nonzeros in only the main
026: *  diagonal and first two superdiagonals.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  N       (input) INTEGER
032: *          The order of the matrix A.
033: *
034: *  DL      (input/output) REAL array, dimension (N-1)
035: *          On entry, DL must contain the (n-1) sub-diagonal elements of
036: *          A.
037: *
038: *          On exit, DL is overwritten by the (n-1) multipliers that
039: *          define the matrix L from the LU factorization of A.
040: *
041: *  D       (input/output) REAL array, dimension (N)
042: *          On entry, D must contain the diagonal elements of A.
043: *
044: *          On exit, D is overwritten by the n diagonal elements of the
045: *          upper triangular matrix U from the LU factorization of A.
046: *
047: *  DU      (input/output) REAL array, dimension (N-1)
048: *          On entry, DU must contain the (n-1) super-diagonal elements
049: *          of A.
050: *
051: *          On exit, DU is overwritten by the (n-1) elements of the first
052: *          super-diagonal of U.
053: *
054: *  DU2     (output) REAL array, dimension (N-2)
055: *          On exit, DU2 is overwritten by the (n-2) elements of the
056: *          second super-diagonal of U.
057: *
058: *  IPIV    (output) INTEGER array, dimension (N)
059: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
060: *          interchanged with row IPIV(i).  IPIV(i) will always be either
061: *          i or i+1; IPIV(i) = i indicates a row interchange was not
062: *          required.
063: *
064: *  INFO    (output) INTEGER
065: *          = 0:  successful exit
066: *          < 0:  if INFO = -k, the k-th argument had an illegal value
067: *          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
068: *                has been completed, but the factor U is exactly
069: *                singular, and division by zero will occur if it is used
070: *                to solve a system of equations.
071: *
072: *  =====================================================================
073: *
074: *     .. Parameters ..
075:       REAL               ZERO
076:       PARAMETER          ( ZERO = 0.0E+0 )
077: *     ..
078: *     .. Local Scalars ..
079:       INTEGER            I
080:       REAL               FACT, TEMP
081: *     ..
082: *     .. Intrinsic Functions ..
083:       INTRINSIC          ABS
084: *     ..
085: *     .. External Subroutines ..
086:       EXTERNAL           XERBLA
087: *     ..
088: *     .. Executable Statements ..
089: *
090:       INFO = 0
091:       IF( N.LT.0 ) THEN
092:          INFO = -1
093:          CALL XERBLA( 'SGTTRF', -INFO )
094:          RETURN
095:       END IF
096: *
097: *     Quick return if possible
098: *
099:       IF( N.EQ.0 )
100:      $   RETURN
101: *
102: *     Initialize IPIV(i) = i and DU2(I) = 0
103: *
104:       DO 10 I = 1, N
105:          IPIV( I ) = I
106:    10 CONTINUE
107:       DO 20 I = 1, N - 2
108:          DU2( I ) = ZERO
109:    20 CONTINUE
110: *
111:       DO 30 I = 1, N - 2
112:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
113: *
114: *           No row interchange required, eliminate DL(I)
115: *
116:             IF( D( I ).NE.ZERO ) THEN
117:                FACT = DL( I ) / D( I )
118:                DL( I ) = FACT
119:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
120:             END IF
121:          ELSE
122: *
123: *           Interchange rows I and I+1, eliminate DL(I)
124: *
125:             FACT = D( I ) / DL( I )
126:             D( I ) = DL( I )
127:             DL( I ) = FACT
128:             TEMP = DU( I )
129:             DU( I ) = D( I+1 )
130:             D( I+1 ) = TEMP - FACT*D( I+1 )
131:             DU2( I ) = DU( I+1 )
132:             DU( I+1 ) = -FACT*DU( I+1 )
133:             IPIV( I ) = I + 1
134:          END IF
135:    30 CONTINUE
136:       IF( N.GT.1 ) THEN
137:          I = N - 1
138:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
139:             IF( D( I ).NE.ZERO ) THEN
140:                FACT = DL( I ) / D( I )
141:                DL( I ) = FACT
142:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
143:             END IF
144:          ELSE
145:             FACT = D( I ) / DL( I )
146:             D( I ) = DL( I )
147:             DL( I ) = FACT
148:             TEMP = DU( I )
149:             DU( I ) = D( I+1 )
150:             D( I+1 ) = TEMP - FACT*D( I+1 )
151:             IPIV( I ) = I + 1
152:          END IF
153:       END IF
154: *
155: *     Check for a zero on the diagonal of U.
156: *
157:       DO 40 I = 1, N
158:          IF( D( I ).EQ.ZERO ) THEN
159:             INFO = I
160:             GO TO 50
161:          END IF
162:    40 CONTINUE
163:    50 CONTINUE
164: *
165:       RETURN
166: *
167: *     End of SGTTRF
168: *
169:       END
170: