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