1 SUBROUTINE ZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF,
2 $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
3 $ WORK, RWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER FACT, TRANS
12 INTEGER INFO, LDB, LDX, N, NRHS
13 DOUBLE PRECISION RCOND
14 * ..
15 * .. Array Arguments ..
16 INTEGER IPIV( * )
17 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
18 COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ),
19 $ DLF( * ), DU( * ), DU2( * ), DUF( * ),
20 $ WORK( * ), X( LDX, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZGTSVX uses the LU factorization to compute the solution to a complex
27 * system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
28 * where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
29 * matrices.
30 *
31 * Error bounds on the solution and a condition estimate are also
32 * provided.
33 *
34 * Description
35 * ===========
36 *
37 * The following steps are performed:
38 *
39 * 1. If FACT = 'N', the LU decomposition is used to factor the matrix A
40 * as A = L * U, where L is a product of permutation and unit lower
41 * bidiagonal matrices and U is upper triangular with nonzeros in
42 * only the main diagonal and first two superdiagonals.
43 *
44 * 2. If some U(i,i)=0, so that U is exactly singular, then the routine
45 * returns with INFO = i. Otherwise, the factored form of A is used
46 * to estimate the condition number of the matrix A. If the
47 * reciprocal of the condition number is less than machine precision,
48 * INFO = N+1 is returned as a warning, but the routine still goes on
49 * to solve for X and compute error bounds as described below.
50 *
51 * 3. The system of equations is solved for X using the factored form
52 * of A.
53 *
54 * 4. Iterative refinement is applied to improve the computed solution
55 * matrix and calculate error bounds and backward error estimates
56 * for it.
57 *
58 * Arguments
59 * =========
60 *
61 * FACT (input) CHARACTER*1
62 * Specifies whether or not the factored form of A has been
63 * supplied on entry.
64 * = 'F': DLF, DF, DUF, DU2, and IPIV contain the factored form
65 * of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not
66 * be modified.
67 * = 'N': The matrix will be copied to DLF, DF, and DUF
68 * and factored.
69 *
70 * TRANS (input) CHARACTER*1
71 * Specifies the form of the system of equations:
72 * = 'N': A * X = B (No transpose)
73 * = 'T': A**T * X = B (Transpose)
74 * = 'C': A**H * X = B (Conjugate transpose)
75 *
76 * N (input) INTEGER
77 * The order of the matrix A. N >= 0.
78 *
79 * NRHS (input) INTEGER
80 * The number of right hand sides, i.e., the number of columns
81 * of the matrix B. NRHS >= 0.
82 *
83 * DL (input) COMPLEX*16 array, dimension (N-1)
84 * The (n-1) subdiagonal elements of A.
85 *
86 * D (input) COMPLEX*16 array, dimension (N)
87 * The n diagonal elements of A.
88 *
89 * DU (input) COMPLEX*16 array, dimension (N-1)
90 * The (n-1) superdiagonal elements of A.
91 *
92 * DLF (input or output) COMPLEX*16 array, dimension (N-1)
93 * If FACT = 'F', then DLF is an input argument and on entry
94 * contains the (n-1) multipliers that define the matrix L from
95 * the LU factorization of A as computed by ZGTTRF.
96 *
97 * If FACT = 'N', then DLF is an output argument and on exit
98 * contains the (n-1) multipliers that define the matrix L from
99 * the LU factorization of A.
100 *
101 * DF (input or output) COMPLEX*16 array, dimension (N)
102 * If FACT = 'F', then DF is an input argument and on entry
103 * contains the n diagonal elements of the upper triangular
104 * matrix U from the LU factorization of A.
105 *
106 * If FACT = 'N', then DF is an output argument and on exit
107 * contains the n diagonal elements of the upper triangular
108 * matrix U from the LU factorization of A.
109 *
110 * DUF (input or output) COMPLEX*16 array, dimension (N-1)
111 * If FACT = 'F', then DUF is an input argument and on entry
112 * contains the (n-1) elements of the first superdiagonal of U.
113 *
114 * If FACT = 'N', then DUF is an output argument and on exit
115 * contains the (n-1) elements of the first superdiagonal of U.
116 *
117 * DU2 (input or output) COMPLEX*16 array, dimension (N-2)
118 * If FACT = 'F', then DU2 is an input argument and on entry
119 * contains the (n-2) elements of the second superdiagonal of
120 * U.
121 *
122 * If FACT = 'N', then DU2 is an output argument and on exit
123 * contains the (n-2) elements of the second superdiagonal of
124 * U.
125 *
126 * IPIV (input or output) INTEGER array, dimension (N)
127 * If FACT = 'F', then IPIV is an input argument and on entry
128 * contains the pivot indices from the LU factorization of A as
129 * computed by ZGTTRF.
130 *
131 * If FACT = 'N', then IPIV is an output argument and on exit
132 * contains the pivot indices from the LU factorization of A;
133 * row i of the matrix was interchanged with row IPIV(i).
134 * IPIV(i) will always be either i or i+1; IPIV(i) = i indicates
135 * a row interchange was not required.
136 *
137 * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
138 * The N-by-NRHS right hand side matrix B.
139 *
140 * LDB (input) INTEGER
141 * The leading dimension of the array B. LDB >= max(1,N).
142 *
143 * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
144 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
145 *
146 * LDX (input) INTEGER
147 * The leading dimension of the array X. LDX >= max(1,N).
148 *
149 * RCOND (output) DOUBLE PRECISION
150 * The estimate of the reciprocal condition number of the matrix
151 * A. If RCOND is less than the machine precision (in
152 * particular, if RCOND = 0), the matrix is singular to working
153 * precision. This condition is indicated by a return code of
154 * INFO > 0.
155 *
156 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
157 * The estimated forward error bound for each solution vector
158 * X(j) (the j-th column of the solution matrix X).
159 * If XTRUE is the true solution corresponding to X(j), FERR(j)
160 * is an estimated upper bound for the magnitude of the largest
161 * element in (X(j) - XTRUE) divided by the magnitude of the
162 * largest element in X(j). The estimate is as reliable as
163 * the estimate for RCOND, and is almost always a slight
164 * overestimate of the true error.
165 *
166 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
167 * The componentwise relative backward error of each solution
168 * vector X(j) (i.e., the smallest relative change in
169 * any element of A or B that makes X(j) an exact solution).
170 *
171 * WORK (workspace) COMPLEX*16 array, dimension (2*N)
172 *
173 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
174 *
175 * INFO (output) INTEGER
176 * = 0: successful exit
177 * < 0: if INFO = -i, the i-th argument had an illegal value
178 * > 0: if INFO = i, and i is
179 * <= N: U(i,i) is exactly zero. The factorization
180 * has not been completed unless i = N, but the
181 * factor U is exactly singular, so the solution
182 * and error bounds could not be computed.
183 * RCOND = 0 is returned.
184 * = N+1: U is nonsingular, but RCOND is less than machine
185 * precision, meaning that the matrix is singular
186 * to working precision. Nevertheless, the
187 * solution and error bounds are computed because
188 * there are a number of situations where the
189 * computed solution can be more accurate than the
190 * value of RCOND would suggest.
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195 DOUBLE PRECISION ZERO
196 PARAMETER ( ZERO = 0.0D+0 )
197 * ..
198 * .. Local Scalars ..
199 LOGICAL NOFACT, NOTRAN
200 CHARACTER NORM
201 DOUBLE PRECISION ANORM
202 * ..
203 * .. External Functions ..
204 LOGICAL LSAME
205 DOUBLE PRECISION DLAMCH, ZLANGT
206 EXTERNAL LSAME, DLAMCH, ZLANGT
207 * ..
208 * .. External Subroutines ..
209 EXTERNAL XERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS,
210 $ ZLACPY
211 * ..
212 * .. Intrinsic Functions ..
213 INTRINSIC MAX
214 * ..
215 * .. Executable Statements ..
216 *
217 INFO = 0
218 NOFACT = LSAME( FACT, 'N' )
219 NOTRAN = LSAME( TRANS, 'N' )
220 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
221 INFO = -1
222 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
223 $ LSAME( TRANS, 'C' ) ) THEN
224 INFO = -2
225 ELSE IF( N.LT.0 ) THEN
226 INFO = -3
227 ELSE IF( NRHS.LT.0 ) THEN
228 INFO = -4
229 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
230 INFO = -14
231 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
232 INFO = -16
233 END IF
234 IF( INFO.NE.0 ) THEN
235 CALL XERBLA( 'ZGTSVX', -INFO )
236 RETURN
237 END IF
238 *
239 IF( NOFACT ) THEN
240 *
241 * Compute the LU factorization of A.
242 *
243 CALL ZCOPY( N, D, 1, DF, 1 )
244 IF( N.GT.1 ) THEN
245 CALL ZCOPY( N-1, DL, 1, DLF, 1 )
246 CALL ZCOPY( N-1, DU, 1, DUF, 1 )
247 END IF
248 CALL ZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
249 *
250 * Return if INFO is non-zero.
251 *
252 IF( INFO.GT.0 )THEN
253 RCOND = ZERO
254 RETURN
255 END IF
256 END IF
257 *
258 * Compute the norm of the matrix A.
259 *
260 IF( NOTRAN ) THEN
261 NORM = '1'
262 ELSE
263 NORM = 'I'
264 END IF
265 ANORM = ZLANGT( NORM, N, DL, D, DU )
266 *
267 * Compute the reciprocal of the condition number of A.
268 *
269 CALL ZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
270 $ INFO )
271 *
272 * Compute the solution vectors X.
273 *
274 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
275 CALL ZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
276 $ INFO )
277 *
278 * Use iterative refinement to improve the computed solutions and
279 * compute error bounds and backward error estimates for them.
280 *
281 CALL ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
282 $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
283 *
284 * Set INFO = N+1 if the matrix is singular to working precision.
285 *
286 IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
287 $ INFO = N + 1
288 *
289 RETURN
290 *
291 * End of ZGTSVX
292 *
293 END
2 $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR,
3 $ WORK, RWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER FACT, TRANS
12 INTEGER INFO, LDB, LDX, N, NRHS
13 DOUBLE PRECISION RCOND
14 * ..
15 * .. Array Arguments ..
16 INTEGER IPIV( * )
17 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
18 COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ),
19 $ DLF( * ), DU( * ), DU2( * ), DUF( * ),
20 $ WORK( * ), X( LDX, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZGTSVX uses the LU factorization to compute the solution to a complex
27 * system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
28 * where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
29 * matrices.
30 *
31 * Error bounds on the solution and a condition estimate are also
32 * provided.
33 *
34 * Description
35 * ===========
36 *
37 * The following steps are performed:
38 *
39 * 1. If FACT = 'N', the LU decomposition is used to factor the matrix A
40 * as A = L * U, where L is a product of permutation and unit lower
41 * bidiagonal matrices and U is upper triangular with nonzeros in
42 * only the main diagonal and first two superdiagonals.
43 *
44 * 2. If some U(i,i)=0, so that U is exactly singular, then the routine
45 * returns with INFO = i. Otherwise, the factored form of A is used
46 * to estimate the condition number of the matrix A. If the
47 * reciprocal of the condition number is less than machine precision,
48 * INFO = N+1 is returned as a warning, but the routine still goes on
49 * to solve for X and compute error bounds as described below.
50 *
51 * 3. The system of equations is solved for X using the factored form
52 * of A.
53 *
54 * 4. Iterative refinement is applied to improve the computed solution
55 * matrix and calculate error bounds and backward error estimates
56 * for it.
57 *
58 * Arguments
59 * =========
60 *
61 * FACT (input) CHARACTER*1
62 * Specifies whether or not the factored form of A has been
63 * supplied on entry.
64 * = 'F': DLF, DF, DUF, DU2, and IPIV contain the factored form
65 * of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not
66 * be modified.
67 * = 'N': The matrix will be copied to DLF, DF, and DUF
68 * and factored.
69 *
70 * TRANS (input) CHARACTER*1
71 * Specifies the form of the system of equations:
72 * = 'N': A * X = B (No transpose)
73 * = 'T': A**T * X = B (Transpose)
74 * = 'C': A**H * X = B (Conjugate transpose)
75 *
76 * N (input) INTEGER
77 * The order of the matrix A. N >= 0.
78 *
79 * NRHS (input) INTEGER
80 * The number of right hand sides, i.e., the number of columns
81 * of the matrix B. NRHS >= 0.
82 *
83 * DL (input) COMPLEX*16 array, dimension (N-1)
84 * The (n-1) subdiagonal elements of A.
85 *
86 * D (input) COMPLEX*16 array, dimension (N)
87 * The n diagonal elements of A.
88 *
89 * DU (input) COMPLEX*16 array, dimension (N-1)
90 * The (n-1) superdiagonal elements of A.
91 *
92 * DLF (input or output) COMPLEX*16 array, dimension (N-1)
93 * If FACT = 'F', then DLF is an input argument and on entry
94 * contains the (n-1) multipliers that define the matrix L from
95 * the LU factorization of A as computed by ZGTTRF.
96 *
97 * If FACT = 'N', then DLF is an output argument and on exit
98 * contains the (n-1) multipliers that define the matrix L from
99 * the LU factorization of A.
100 *
101 * DF (input or output) COMPLEX*16 array, dimension (N)
102 * If FACT = 'F', then DF is an input argument and on entry
103 * contains the n diagonal elements of the upper triangular
104 * matrix U from the LU factorization of A.
105 *
106 * If FACT = 'N', then DF is an output argument and on exit
107 * contains the n diagonal elements of the upper triangular
108 * matrix U from the LU factorization of A.
109 *
110 * DUF (input or output) COMPLEX*16 array, dimension (N-1)
111 * If FACT = 'F', then DUF is an input argument and on entry
112 * contains the (n-1) elements of the first superdiagonal of U.
113 *
114 * If FACT = 'N', then DUF is an output argument and on exit
115 * contains the (n-1) elements of the first superdiagonal of U.
116 *
117 * DU2 (input or output) COMPLEX*16 array, dimension (N-2)
118 * If FACT = 'F', then DU2 is an input argument and on entry
119 * contains the (n-2) elements of the second superdiagonal of
120 * U.
121 *
122 * If FACT = 'N', then DU2 is an output argument and on exit
123 * contains the (n-2) elements of the second superdiagonal of
124 * U.
125 *
126 * IPIV (input or output) INTEGER array, dimension (N)
127 * If FACT = 'F', then IPIV is an input argument and on entry
128 * contains the pivot indices from the LU factorization of A as
129 * computed by ZGTTRF.
130 *
131 * If FACT = 'N', then IPIV is an output argument and on exit
132 * contains the pivot indices from the LU factorization of A;
133 * row i of the matrix was interchanged with row IPIV(i).
134 * IPIV(i) will always be either i or i+1; IPIV(i) = i indicates
135 * a row interchange was not required.
136 *
137 * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
138 * The N-by-NRHS right hand side matrix B.
139 *
140 * LDB (input) INTEGER
141 * The leading dimension of the array B. LDB >= max(1,N).
142 *
143 * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
144 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
145 *
146 * LDX (input) INTEGER
147 * The leading dimension of the array X. LDX >= max(1,N).
148 *
149 * RCOND (output) DOUBLE PRECISION
150 * The estimate of the reciprocal condition number of the matrix
151 * A. If RCOND is less than the machine precision (in
152 * particular, if RCOND = 0), the matrix is singular to working
153 * precision. This condition is indicated by a return code of
154 * INFO > 0.
155 *
156 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
157 * The estimated forward error bound for each solution vector
158 * X(j) (the j-th column of the solution matrix X).
159 * If XTRUE is the true solution corresponding to X(j), FERR(j)
160 * is an estimated upper bound for the magnitude of the largest
161 * element in (X(j) - XTRUE) divided by the magnitude of the
162 * largest element in X(j). The estimate is as reliable as
163 * the estimate for RCOND, and is almost always a slight
164 * overestimate of the true error.
165 *
166 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
167 * The componentwise relative backward error of each solution
168 * vector X(j) (i.e., the smallest relative change in
169 * any element of A or B that makes X(j) an exact solution).
170 *
171 * WORK (workspace) COMPLEX*16 array, dimension (2*N)
172 *
173 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
174 *
175 * INFO (output) INTEGER
176 * = 0: successful exit
177 * < 0: if INFO = -i, the i-th argument had an illegal value
178 * > 0: if INFO = i, and i is
179 * <= N: U(i,i) is exactly zero. The factorization
180 * has not been completed unless i = N, but the
181 * factor U is exactly singular, so the solution
182 * and error bounds could not be computed.
183 * RCOND = 0 is returned.
184 * = N+1: U is nonsingular, but RCOND is less than machine
185 * precision, meaning that the matrix is singular
186 * to working precision. Nevertheless, the
187 * solution and error bounds are computed because
188 * there are a number of situations where the
189 * computed solution can be more accurate than the
190 * value of RCOND would suggest.
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195 DOUBLE PRECISION ZERO
196 PARAMETER ( ZERO = 0.0D+0 )
197 * ..
198 * .. Local Scalars ..
199 LOGICAL NOFACT, NOTRAN
200 CHARACTER NORM
201 DOUBLE PRECISION ANORM
202 * ..
203 * .. External Functions ..
204 LOGICAL LSAME
205 DOUBLE PRECISION DLAMCH, ZLANGT
206 EXTERNAL LSAME, DLAMCH, ZLANGT
207 * ..
208 * .. External Subroutines ..
209 EXTERNAL XERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS,
210 $ ZLACPY
211 * ..
212 * .. Intrinsic Functions ..
213 INTRINSIC MAX
214 * ..
215 * .. Executable Statements ..
216 *
217 INFO = 0
218 NOFACT = LSAME( FACT, 'N' )
219 NOTRAN = LSAME( TRANS, 'N' )
220 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
221 INFO = -1
222 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
223 $ LSAME( TRANS, 'C' ) ) THEN
224 INFO = -2
225 ELSE IF( N.LT.0 ) THEN
226 INFO = -3
227 ELSE IF( NRHS.LT.0 ) THEN
228 INFO = -4
229 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
230 INFO = -14
231 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
232 INFO = -16
233 END IF
234 IF( INFO.NE.0 ) THEN
235 CALL XERBLA( 'ZGTSVX', -INFO )
236 RETURN
237 END IF
238 *
239 IF( NOFACT ) THEN
240 *
241 * Compute the LU factorization of A.
242 *
243 CALL ZCOPY( N, D, 1, DF, 1 )
244 IF( N.GT.1 ) THEN
245 CALL ZCOPY( N-1, DL, 1, DLF, 1 )
246 CALL ZCOPY( N-1, DU, 1, DUF, 1 )
247 END IF
248 CALL ZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO )
249 *
250 * Return if INFO is non-zero.
251 *
252 IF( INFO.GT.0 )THEN
253 RCOND = ZERO
254 RETURN
255 END IF
256 END IF
257 *
258 * Compute the norm of the matrix A.
259 *
260 IF( NOTRAN ) THEN
261 NORM = '1'
262 ELSE
263 NORM = 'I'
264 END IF
265 ANORM = ZLANGT( NORM, N, DL, D, DU )
266 *
267 * Compute the reciprocal of the condition number of A.
268 *
269 CALL ZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK,
270 $ INFO )
271 *
272 * Compute the solution vectors X.
273 *
274 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
275 CALL ZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX,
276 $ INFO )
277 *
278 * Use iterative refinement to improve the computed solutions and
279 * compute error bounds and backward error estimates for them.
280 *
281 CALL ZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV,
282 $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
283 *
284 * Set INFO = N+1 if the matrix is singular to working precision.
285 *
286 IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
287 $ INFO = N + 1
288 *
289 RETURN
290 *
291 * End of ZGTSVX
292 *
293 END