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