1 SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 DOUBLE PRECISION D( * ), DL( * ), DU( * ), DU2( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DGTTRF computes an LU factorization of a real tridiagonal matrix A
20 * using elimination with partial pivoting and row interchanges.
21 *
22 * The factorization has the form
23 * A = L * U
24 * where L is a product of permutation and unit lower bidiagonal
25 * matrices and U is upper triangular with nonzeros in only the main
26 * diagonal and first two superdiagonals.
27 *
28 * Arguments
29 * =========
30 *
31 * N (input) INTEGER
32 * The order of the matrix A.
33 *
34 * DL (input/output) DOUBLE PRECISION array, dimension (N-1)
35 * On entry, DL must contain the (n-1) sub-diagonal elements of
36 * A.
37 *
38 * On exit, DL is overwritten by the (n-1) multipliers that
39 * define the matrix L from the LU factorization of A.
40 *
41 * D (input/output) DOUBLE PRECISION array, dimension (N)
42 * On entry, D must contain the diagonal elements of A.
43 *
44 * On exit, D is overwritten by the n diagonal elements of the
45 * upper triangular matrix U from the LU factorization of A.
46 *
47 * DU (input/output) DOUBLE PRECISION array, dimension (N-1)
48 * On entry, DU must contain the (n-1) super-diagonal elements
49 * of A.
50 *
51 * On exit, DU is overwritten by the (n-1) elements of the first
52 * super-diagonal of U.
53 *
54 * DU2 (output) DOUBLE PRECISION array, dimension (N-2)
55 * On exit, DU2 is overwritten by the (n-2) elements of the
56 * second super-diagonal of U.
57 *
58 * IPIV (output) INTEGER array, dimension (N)
59 * The pivot indices; for 1 <= i <= n, row i of the matrix was
60 * interchanged with row IPIV(i). IPIV(i) will always be either
61 * i or i+1; IPIV(i) = i indicates a row interchange was not
62 * required.
63 *
64 * INFO (output) INTEGER
65 * = 0: successful exit
66 * < 0: if INFO = -k, the k-th argument had an illegal value
67 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
68 * has been completed, but the factor U is exactly
69 * singular, and division by zero will occur if it is used
70 * to solve a system of equations.
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ZERO
76 PARAMETER ( ZERO = 0.0D+0 )
77 * ..
78 * .. Local Scalars ..
79 INTEGER I
80 DOUBLE PRECISION FACT, TEMP
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL XERBLA
87 * ..
88 * .. Executable Statements ..
89 *
90 INFO = 0
91 IF( N.LT.0 ) THEN
92 INFO = -1
93 CALL XERBLA( 'DGTTRF', -INFO )
94 RETURN
95 END IF
96 *
97 * Quick return if possible
98 *
99 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 DGTTRF
168 *
169 END
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 DOUBLE PRECISION D( * ), DL( * ), DU( * ), DU2( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DGTTRF computes an LU factorization of a real tridiagonal matrix A
20 * using elimination with partial pivoting and row interchanges.
21 *
22 * The factorization has the form
23 * A = L * U
24 * where L is a product of permutation and unit lower bidiagonal
25 * matrices and U is upper triangular with nonzeros in only the main
26 * diagonal and first two superdiagonals.
27 *
28 * Arguments
29 * =========
30 *
31 * N (input) INTEGER
32 * The order of the matrix A.
33 *
34 * DL (input/output) DOUBLE PRECISION array, dimension (N-1)
35 * On entry, DL must contain the (n-1) sub-diagonal elements of
36 * A.
37 *
38 * On exit, DL is overwritten by the (n-1) multipliers that
39 * define the matrix L from the LU factorization of A.
40 *
41 * D (input/output) DOUBLE PRECISION array, dimension (N)
42 * On entry, D must contain the diagonal elements of A.
43 *
44 * On exit, D is overwritten by the n diagonal elements of the
45 * upper triangular matrix U from the LU factorization of A.
46 *
47 * DU (input/output) DOUBLE PRECISION array, dimension (N-1)
48 * On entry, DU must contain the (n-1) super-diagonal elements
49 * of A.
50 *
51 * On exit, DU is overwritten by the (n-1) elements of the first
52 * super-diagonal of U.
53 *
54 * DU2 (output) DOUBLE PRECISION array, dimension (N-2)
55 * On exit, DU2 is overwritten by the (n-2) elements of the
56 * second super-diagonal of U.
57 *
58 * IPIV (output) INTEGER array, dimension (N)
59 * The pivot indices; for 1 <= i <= n, row i of the matrix was
60 * interchanged with row IPIV(i). IPIV(i) will always be either
61 * i or i+1; IPIV(i) = i indicates a row interchange was not
62 * required.
63 *
64 * INFO (output) INTEGER
65 * = 0: successful exit
66 * < 0: if INFO = -k, the k-th argument had an illegal value
67 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
68 * has been completed, but the factor U is exactly
69 * singular, and division by zero will occur if it is used
70 * to solve a system of equations.
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ZERO
76 PARAMETER ( ZERO = 0.0D+0 )
77 * ..
78 * .. Local Scalars ..
79 INTEGER I
80 DOUBLE PRECISION FACT, TEMP
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL XERBLA
87 * ..
88 * .. Executable Statements ..
89 *
90 INFO = 0
91 IF( N.LT.0 ) THEN
92 INFO = -1
93 CALL XERBLA( 'DGTTRF', -INFO )
94 RETURN
95 END IF
96 *
97 * Quick return if possible
98 *
99 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 DGTTRF
168 *
169 END