Mercurial > octave-nkf

comparison libcruft/lapack/zgtts2.f @ 7053:570a382ce556

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-23 23:17:36 by jwe]
author jwe
date Tue, 23 Oct 2007 23:17:36 +0000
parents
children
comparison
equal deleted inserted replaced
7052:ee70ac66041f 7053:570a382ce556
1 SUBROUTINE ZGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER ITRANS, LDB, N, NRHS
9 * ..
10 * .. Array Arguments ..
11 INTEGER IPIV( * )
12 COMPLEX*16 B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * ZGTTS2 solves one of the systems of equations
19 * A * X = B, A**T * X = B, or A**H * X = B,
20 * with a tridiagonal matrix A using the LU factorization computed
21 * by ZGTTRF.
22 *
23 * Arguments
24 * =========
25 *
26 * ITRANS (input) INTEGER
27 * Specifies the form of the system of equations.
28 * = 0: A * X = B (No transpose)
29 * = 1: A**T * X = B (Transpose)
30 * = 2: A**H * X = B (Conjugate transpose)
31 *
32 * N (input) INTEGER
33 * The order of the matrix A.
34 *
35 * NRHS (input) INTEGER
36 * The number of right hand sides, i.e., the number of columns
37 * of the matrix B. NRHS >= 0.
38 *
39 * DL (input) COMPLEX*16 array, dimension (N-1)
40 * The (n-1) multipliers that define the matrix L from the
41 * LU factorization of A.
42 *
43 * D (input) COMPLEX*16 array, dimension (N)
44 * The n diagonal elements of the upper triangular matrix U from
45 * the LU factorization of A.
46 *
47 * DU (input) COMPLEX*16 array, dimension (N-1)
48 * The (n-1) elements of the first super-diagonal of U.
49 *
50 * DU2 (input) COMPLEX*16 array, dimension (N-2)
51 * The (n-2) elements of the second super-diagonal of U.
52 *
53 * IPIV (input) INTEGER array, dimension (N)
54 * The pivot indices; for 1 <= i <= n, row i of the matrix was
55 * interchanged with row IPIV(i). IPIV(i) will always be either
56 * i or i+1; IPIV(i) = i indicates a row interchange was not
57 * required.
58 *
59 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
60 * On entry, the matrix of right hand side vectors B.
61 * On exit, B is overwritten by the solution vectors X.
62 *
63 * LDB (input) INTEGER
64 * The leading dimension of the array B. LDB >= max(1,N).
65 *
66 * =====================================================================
67 *
68 * .. Local Scalars ..
69 INTEGER I, J
70 COMPLEX*16 TEMP
71 * ..
72 * .. Intrinsic Functions ..
73 INTRINSIC DCONJG
74 * ..
75 * .. Executable Statements ..
76 *
77 * Quick return if possible
78 *
79 IF( N.EQ.0 .OR. NRHS.EQ.0 )
80 $ RETURN
81 *
82 IF( ITRANS.EQ.0 ) THEN
83 *
84 * Solve A*X = B using the LU factorization of A,
85 * overwriting each right hand side vector with its solution.
86 *
87 IF( NRHS.LE.1 ) THEN
88 J = 1
89 10 CONTINUE
90 *
91 * Solve L*x = b.
92 *
93 DO 20 I = 1, N - 1
94 IF( IPIV( I ).EQ.I ) THEN
95 B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
96 ELSE
97 TEMP = B( I, J )
98 B( I, J ) = B( I+1, J )
99 B( I+1, J ) = TEMP - DL( I )*B( I, J )
100 END IF
101 20 CONTINUE
102 *
103 * Solve U*x = b.
104 *
105 B( N, J ) = B( N, J ) / D( N )
106 IF( N.GT.1 )
107 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
108 $ D( N-1 )
109 DO 30 I = N - 2, 1, -1
110 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
111 $ B( I+2, J ) ) / D( I )
112 30 CONTINUE
113 IF( J.LT.NRHS ) THEN
114 J = J + 1
115 GO TO 10
116 END IF
117 ELSE
118 DO 60 J = 1, NRHS
119 *
120 * Solve L*x = b.
121 *
122 DO 40 I = 1, N - 1
123 IF( IPIV( I ).EQ.I ) THEN
124 B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
125 ELSE
126 TEMP = B( I, J )
127 B( I, J ) = B( I+1, J )
128 B( I+1, J ) = TEMP - DL( I )*B( I, J )
129 END IF
130 40 CONTINUE
131 *
132 * Solve U*x = b.
133 *
134 B( N, J ) = B( N, J ) / D( N )
135 IF( N.GT.1 )
136 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
137 $ D( N-1 )
138 DO 50 I = N - 2, 1, -1
139 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
140 $ B( I+2, J ) ) / D( I )
141 50 CONTINUE
142 60 CONTINUE
143 END IF
144 ELSE IF( ITRANS.EQ.1 ) THEN
145 *
146 * Solve A**T * X = B.
147 *
148 IF( NRHS.LE.1 ) THEN
149 J = 1
150 70 CONTINUE
151 *
152 * Solve U**T * x = b.
153 *
154 B( 1, J ) = B( 1, J ) / D( 1 )
155 IF( N.GT.1 )
156 $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
157 DO 80 I = 3, N
158 B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
159 $ B( I-2, J ) ) / D( I )
160 80 CONTINUE
161 *
162 * Solve L**T * x = b.
163 *
164 DO 90 I = N - 1, 1, -1
165 IF( IPIV( I ).EQ.I ) THEN
166 B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
167 ELSE
168 TEMP = B( I+1, J )
169 B( I+1, J ) = B( I, J ) - DL( I )*TEMP
170 B( I, J ) = TEMP
171 END IF
172 90 CONTINUE
173 IF( J.LT.NRHS ) THEN
174 J = J + 1
175 GO TO 70
176 END IF
177 ELSE
178 DO 120 J = 1, NRHS
179 *
180 * Solve U**T * x = b.
181 *
182 B( 1, J ) = B( 1, J ) / D( 1 )
183 IF( N.GT.1 )
184 $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
185 DO 100 I = 3, N
186 B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
187 $ DU2( I-2 )*B( I-2, J ) ) / D( I )
188 100 CONTINUE
189 *
190 * Solve L**T * x = b.
191 *
192 DO 110 I = N - 1, 1, -1
193 IF( IPIV( I ).EQ.I ) THEN
194 B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
195 ELSE
196 TEMP = B( I+1, J )
197 B( I+1, J ) = B( I, J ) - DL( I )*TEMP
198 B( I, J ) = TEMP
199 END IF
200 110 CONTINUE
201 120 CONTINUE
202 END IF
203 ELSE
204 *
205 * Solve A**H * X = B.
206 *
207 IF( NRHS.LE.1 ) THEN
208 J = 1
209 130 CONTINUE
210 *
211 * Solve U**H * x = b.
212 *
213 B( 1, J ) = B( 1, J ) / DCONJG( D( 1 ) )
214 IF( N.GT.1 )
215 $ B( 2, J ) = ( B( 2, J )-DCONJG( DU( 1 ) )*B( 1, J ) ) /
216 $ DCONJG( D( 2 ) )
217 DO 140 I = 3, N
218 B( I, J ) = ( B( I, J )-DCONJG( DU( I-1 ) )*B( I-1, J )-
219 $ DCONJG( DU2( I-2 ) )*B( I-2, J ) ) /
220 $ DCONJG( D( I ) )
221 140 CONTINUE
222 *
223 * Solve L**H * x = b.
224 *
225 DO 150 I = N - 1, 1, -1
226 IF( IPIV( I ).EQ.I ) THEN
227 B( I, J ) = B( I, J ) - DCONJG( DL( I ) )*B( I+1, J )
228 ELSE
229 TEMP = B( I+1, J )
230 B( I+1, J ) = B( I, J ) - DCONJG( DL( I ) )*TEMP
231 B( I, J ) = TEMP
232 END IF
233 150 CONTINUE
234 IF( J.LT.NRHS ) THEN
235 J = J + 1
236 GO TO 130
237 END IF
238 ELSE
239 DO 180 J = 1, NRHS
240 *
241 * Solve U**H * x = b.
242 *
243 B( 1, J ) = B( 1, J ) / DCONJG( D( 1 ) )
244 IF( N.GT.1 )
245 $ B( 2, J ) = ( B( 2, J )-DCONJG( DU( 1 ) )*B( 1, J ) )
246 $ / DCONJG( D( 2 ) )
247 DO 160 I = 3, N
248 B( I, J ) = ( B( I, J )-DCONJG( DU( I-1 ) )*
249 $ B( I-1, J )-DCONJG( DU2( I-2 ) )*
250 $ B( I-2, J ) ) / DCONJG( D( I ) )
251 160 CONTINUE
252 *
253 * Solve L**H * x = b.
254 *
255 DO 170 I = N - 1, 1, -1
256 IF( IPIV( I ).EQ.I ) THEN
257 B( I, J ) = B( I, J ) - DCONJG( DL( I ) )*
258 $ B( I+1, J )
259 ELSE
260 TEMP = B( I+1, J )
261 B( I+1, J ) = B( I, J ) - DCONJG( DL( I ) )*TEMP
262 B( I, J ) = TEMP
263 END IF
264 170 CONTINUE
265 180 CONTINUE
266 END IF
267 END IF
268 *
269 * End of ZGTTS2
270 *
271 END