|
0
|
1 SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO ) |
|
|
2 * |
|
|
3 * -- LAPACK routine (version 3.0) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
|
6 * October 31, 1999 |
|
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 INTEGER INFO, LDB, N, NRHS |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * DGTSV solves the equation |
|
|
19 * |
|
|
20 * A*X = B, |
|
|
21 * |
|
|
22 * where A is an n by n tridiagonal matrix, by Gaussian elimination with |
|
|
23 * partial pivoting. |
|
|
24 * |
|
|
25 * Note that the equation A'*X = B may be solved by interchanging the |
|
|
26 * order of the arguments DU and DL. |
|
|
27 * |
|
|
28 * Arguments |
|
|
29 * ========= |
|
|
30 * |
|
|
31 * N (input) INTEGER |
|
|
32 * The order of the matrix A. N >= 0. |
|
|
33 * |
|
|
34 * NRHS (input) INTEGER |
|
|
35 * The number of right hand sides, i.e., the number of columns |
|
|
36 * of the matrix B. NRHS >= 0. |
|
|
37 * |
|
|
38 * DL (input/output) DOUBLE PRECISION array, dimension (N-1) |
|
|
39 * On entry, DL must contain the (n-1) sub-diagonal elements of |
|
|
40 * A. |
|
|
41 * |
|
|
42 * On exit, DL is overwritten by the (n-2) elements of the |
|
|
43 * second super-diagonal of the upper triangular matrix U from |
|
|
44 * the LU factorization of A, in DL(1), ..., DL(n-2). |
|
|
45 * |
|
|
46 * D (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
47 * On entry, D must contain the diagonal elements of A. |
|
|
48 * |
|
|
49 * On exit, D is overwritten by the n diagonal elements of U. |
|
|
50 * |
|
|
51 * DU (input/output) DOUBLE PRECISION array, dimension (N-1) |
|
|
52 * On entry, DU must contain the (n-1) super-diagonal elements |
|
|
53 * of A. |
|
|
54 * |
|
|
55 * On exit, DU is overwritten by the (n-1) elements of the first |
|
|
56 * super-diagonal of U. |
|
|
57 * |
|
|
58 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) |
|
|
59 * On entry, the N by NRHS matrix of right hand side matrix B. |
|
|
60 * On exit, if INFO = 0, the N by NRHS solution matrix X. |
|
|
61 * |
|
|
62 * LDB (input) INTEGER |
|
|
63 * The leading dimension of the array B. LDB >= max(1,N). |
|
|
64 * |
|
|
65 * INFO (output) INTEGER |
|
|
66 * = 0: successful exit |
|
|
67 * < 0: if INFO = -i, the i-th argument had an illegal value |
|
|
68 * > 0: if INFO = i, U(i,i) is exactly zero, and the solution |
|
|
69 * has not been computed. The factorization has not been |
|
|
70 * completed unless i = N. |
|
|
71 * |
|
|
72 * ===================================================================== |
|
|
73 * |
|
|
74 * .. Parameters .. |
|
|
75 DOUBLE PRECISION ZERO |
|
|
76 PARAMETER ( ZERO = 0.0D+0 ) |
|
|
77 * .. |
|
|
78 * .. Local Scalars .. |
|
|
79 INTEGER I, J |
|
|
80 DOUBLE PRECISION FACT, TEMP |
|
|
81 * .. |
|
|
82 * .. Intrinsic Functions .. |
|
|
83 INTRINSIC ABS, MAX |
|
|
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 ELSE IF( NRHS.LT.0 ) THEN |
|
|
94 INFO = -2 |
|
|
95 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN |
|
|
96 INFO = -7 |
|
|
97 END IF |
|
|
98 IF( INFO.NE.0 ) THEN |
|
|
99 CALL XERBLA( 'DGTSV ', -INFO ) |
|
|
100 RETURN |
|
|
101 END IF |
|
|
102 * |
|
|
103 IF( N.EQ.0 ) |
|
|
104 $ RETURN |
|
|
105 * |
|
|
106 IF( NRHS.EQ.1 ) THEN |
|
|
107 DO 10 I = 1, N - 2 |
|
|
108 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN |
|
|
109 * |
|
|
110 * No row interchange required |
|
|
111 * |
|
|
112 IF( D( I ).NE.ZERO ) THEN |
|
|
113 FACT = DL( I ) / D( I ) |
|
|
114 D( I+1 ) = D( I+1 ) - FACT*DU( I ) |
|
|
115 B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 ) |
|
|
116 ELSE |
|
|
117 INFO = I |
|
|
118 RETURN |
|
|
119 END IF |
|
|
120 DL( I ) = ZERO |
|
|
121 ELSE |
|
|
122 * |
|
|
123 * Interchange rows I and I+1 |
|
|
124 * |
|
|
125 FACT = D( I ) / DL( I ) |
|
|
126 D( I ) = DL( I ) |
|
|
127 TEMP = D( I+1 ) |
|
|
128 D( I+1 ) = DU( I ) - FACT*TEMP |
|
|
129 DL( I ) = DU( I+1 ) |
|
|
130 DU( I+1 ) = -FACT*DL( I ) |
|
|
131 DU( I ) = TEMP |
|
|
132 TEMP = B( I, 1 ) |
|
|
133 B( I, 1 ) = B( I+1, 1 ) |
|
|
134 B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 ) |
|
|
135 END IF |
|
|
136 10 CONTINUE |
|
|
137 IF( N.GT.1 ) THEN |
|
|
138 I = N - 1 |
|
|
139 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN |
|
|
140 IF( D( I ).NE.ZERO ) THEN |
|
|
141 FACT = DL( I ) / D( I ) |
|
|
142 D( I+1 ) = D( I+1 ) - FACT*DU( I ) |
|
|
143 B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 ) |
|
|
144 ELSE |
|
|
145 INFO = I |
|
|
146 RETURN |
|
|
147 END IF |
|
|
148 ELSE |
|
|
149 FACT = D( I ) / DL( I ) |
|
|
150 D( I ) = DL( I ) |
|
|
151 TEMP = D( I+1 ) |
|
|
152 D( I+1 ) = DU( I ) - FACT*TEMP |
|
|
153 DU( I ) = TEMP |
|
|
154 TEMP = B( I, 1 ) |
|
|
155 B( I, 1 ) = B( I+1, 1 ) |
|
|
156 B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 ) |
|
|
157 END IF |
|
|
158 END IF |
|
|
159 IF( D( N ).EQ.ZERO ) THEN |
|
|
160 INFO = N |
|
|
161 RETURN |
|
|
162 END IF |
|
|
163 ELSE |
|
|
164 DO 40 I = 1, N - 2 |
|
|
165 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN |
|
|
166 * |
|
|
167 * No row interchange required |
|
|
168 * |
|
|
169 IF( D( I ).NE.ZERO ) THEN |
|
|
170 FACT = DL( I ) / D( I ) |
|
|
171 D( I+1 ) = D( I+1 ) - FACT*DU( I ) |
|
|
172 DO 20 J = 1, NRHS |
|
|
173 B( I+1, J ) = B( I+1, J ) - FACT*B( I, J ) |
|
|
174 20 CONTINUE |
|
|
175 ELSE |
|
|
176 INFO = I |
|
|
177 RETURN |
|
|
178 END IF |
|
|
179 DL( I ) = ZERO |
|
|
180 ELSE |
|
|
181 * |
|
|
182 * Interchange rows I and I+1 |
|
|
183 * |
|
|
184 FACT = D( I ) / DL( I ) |
|
|
185 D( I ) = DL( I ) |
|
|
186 TEMP = D( I+1 ) |
|
|
187 D( I+1 ) = DU( I ) - FACT*TEMP |
|
|
188 DL( I ) = DU( I+1 ) |
|
|
189 DU( I+1 ) = -FACT*DL( I ) |
|
|
190 DU( I ) = TEMP |
|
|
191 DO 30 J = 1, NRHS |
|
|
192 TEMP = B( I, J ) |
|
|
193 B( I, J ) = B( I+1, J ) |
|
|
194 B( I+1, J ) = TEMP - FACT*B( I+1, J ) |
|
|
195 30 CONTINUE |
|
|
196 END IF |
|
|
197 40 CONTINUE |
|
|
198 IF( N.GT.1 ) THEN |
|
|
199 I = N - 1 |
|
|
200 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN |
|
|
201 IF( D( I ).NE.ZERO ) THEN |
|
|
202 FACT = DL( I ) / D( I ) |
|
|
203 D( I+1 ) = D( I+1 ) - FACT*DU( I ) |
|
|
204 DO 50 J = 1, NRHS |
|
|
205 B( I+1, J ) = B( I+1, J ) - FACT*B( I, J ) |
|
|
206 50 CONTINUE |
|
|
207 ELSE |
|
|
208 INFO = I |
|
|
209 RETURN |
|
|
210 END IF |
|
|
211 ELSE |
|
|
212 FACT = D( I ) / DL( I ) |
|
|
213 D( I ) = DL( I ) |
|
|
214 TEMP = D( I+1 ) |
|
|
215 D( I+1 ) = DU( I ) - FACT*TEMP |
|
|
216 DU( I ) = TEMP |
|
|
217 DO 60 J = 1, NRHS |
|
|
218 TEMP = B( I, J ) |
|
|
219 B( I, J ) = B( I+1, J ) |
|
|
220 B( I+1, J ) = TEMP - FACT*B( I+1, J ) |
|
|
221 60 CONTINUE |
|
|
222 END IF |
|
|
223 END IF |
|
|
224 IF( D( N ).EQ.ZERO ) THEN |
|
|
225 INFO = N |
|
|
226 RETURN |
|
|
227 END IF |
|
|
228 END IF |
|
|
229 * |
|
|
230 * Back solve with the matrix U from the factorization. |
|
|
231 * |
|
|
232 IF( NRHS.LE.2 ) THEN |
|
|
233 J = 1 |
|
|
234 70 CONTINUE |
|
|
235 B( N, J ) = B( N, J ) / D( N ) |
|
|
236 IF( N.GT.1 ) |
|
|
237 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 ) |
|
|
238 DO 80 I = N - 2, 1, -1 |
|
|
239 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )* |
|
|
240 $ B( I+2, J ) ) / D( I ) |
|
|
241 80 CONTINUE |
|
|
242 IF( J.LT.NRHS ) THEN |
|
|
243 J = J + 1 |
|
|
244 GO TO 70 |
|
|
245 END IF |
|
|
246 ELSE |
|
|
247 DO 100 J = 1, NRHS |
|
|
248 B( N, J ) = B( N, J ) / D( N ) |
|
|
249 IF( N.GT.1 ) |
|
|
250 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / |
|
|
251 $ D( N-1 ) |
|
|
252 DO 90 I = N - 2, 1, -1 |
|
|
253 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )* |
|
|
254 $ B( I+2, J ) ) / D( I ) |
|
|
255 90 CONTINUE |
|
|
256 100 CONTINUE |
|
|
257 END IF |
|
|
258 * |
|
|
259 RETURN |
|
|
260 * |
|
|
261 * End of DGTSV |
|
|
262 * |
|
|
263 END |
