Mercurial > forge

comparison main/splines/dgtsv.f @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
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