1ZGESVX(1) LAPACK driver routine (version 3.1) ZGESVX(1)
2
6 ZGESVX - the LU factorization to compute the solution to a complex sys‐
7 tem of linear equations A * X = B,
8
10 SUBROUTINE ZGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED,
11 R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
12 RWORK, INFO )
13 CHARACTER EQUED, FACT, TRANS
15 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
17 DOUBLE PRECISION RCOND
19 INTEGER IPIV( * )
21 DOUBLE PRECISION BERR( * ), C( * ), FERR( * ), R( * ),
23 RWORK( * )
24 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), WORK( * ),
26 X( LDX, * )
27
29 ZGESVX uses the LU factorization to compute the solution to a complex
30 system of linear equations
31 A * X = B, where A is an N-by-N matrix and X and B are N-by-NRHS
32 matrices.
33 Error bounds on the solution and a condition estimate are also pro‐
35 vided.
36
39 The following steps are performed:
40 1. If FACT = 'E', real scaling factors are computed to equilibrate
42 the system:
43 TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
44 TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
45 TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
46 Whether or not the system will be equilibrated depends on the
47 scaling of the matrix A, but if equilibration is used, A is
48 overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
49 or diag(C)*B (if TRANS = 'T' or 'C').
50 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
52 matrix A (after equilibration if FACT = 'E') as
53 A = P * L * U,
54 where P is a permutation matrix, L is a unit lower triangular
55 matrix, and U is upper triangular.
56 3. If some U(i,i)=0, so that U is exactly singular, then the routine
58 returns with INFO = i. Otherwise, the factored form of A is used
59 to estimate the condition number of the matrix A. If the
60 reciprocal of the condition number is less than machine precision,
61 INFO = N+1 is returned as a warning, but the routine still goes on
62 to solve for X and compute error bounds as described below.
63 4. The system of equations is solved for X using the factored form
65 of A.
66 5. Iterative refinement is applied to improve the computed solution
68 matrix and calculate error bounds and backward error estimates
69 for it.
70 6. If equilibration was used, the matrix X is premultiplied by
72 diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
73 that it solves the original system before equilibration.
74
77 FACT (input) CHARACTER*1
78 Specifies whether or not the factored form of the matrix A is
79 supplied on entry, and if not, whether the matrix A should be
80 equilibrated before it is factored. = 'F': On entry, AF and
81 IPIV contain the factored form of A. If EQUED is not 'N', the
82 matrix A has been equilibrated with scaling factors given by R
83 and C. A, AF, and IPIV are not modified. = 'N': The matrix A
84 will be copied to AF and factored.
85 = 'E': The matrix A will be equilibrated if necessary, then
86 copied to AF and factored.
87 TRANS (input) CHARACTER*1
89 Specifies the form of the system of equations:
90 = 'N': A * X = B (No transpose)
91 = 'T': A**T * X = B (Transpose)
92 = 'C': A**H * X = B (Conjugate transpose)
93 N (input) INTEGER
95 The number of linear equations, i.e., the order of the matrix
96 A. N >= 0.
97 NRHS (input) INTEGER
99 The number of right hand sides, i.e., the number of columns of
100 the matrices B and X. NRHS >= 0.
101 A (input/output) COMPLEX*16 array, dimension (LDA,N)
103 On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is not
104 'N', then A must have been equilibrated by the scaling factors
105 in R and/or C. A is not modified if FACT = 'F' or 'N', or if
106 FACT = 'E' and EQUED = 'N' on exit.
107 On exit, if EQUED .ne. 'N', A is scaled as follows: EQUED =
109 'R': A := diag(R) * A
110 EQUED = 'C': A := A * diag(C)
111 EQUED = 'B': A := diag(R) * A * diag(C).
112 LDA (input) INTEGER
114 The leading dimension of the array A. LDA >= max(1,N).
115 AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
117 If FACT = 'F', then AF is an input argument and on entry con‐
118 tains the factors L and U from the factorization A = P*L*U as
119 computed by ZGETRF. If EQUED .ne. 'N', then AF is the factored
120 form of the equilibrated matrix A.
121 If FACT = 'N', then AF is an output argument and on exit
123 returns the factors L and U from the factorization A = P*L*U of
124 the original matrix A.
125 If FACT = 'E', then AF is an output argument and on exit
127 returns the factors L and U from the factorization A = P*L*U of
128 the equilibrated matrix A (see the description of A for the
129 form of the equilibrated matrix).
130 LDAF (input) INTEGER
132 The leading dimension of the array AF. LDAF >= max(1,N).
133 IPIV (input or output) INTEGER array, dimension (N)
135 If FACT = 'F', then IPIV is an input argument and on entry con‐
136 tains the pivot indices from the factorization A = P*L*U as
137 computed by ZGETRF; row i of the matrix was interchanged with
138 row IPIV(i).
139 If FACT = 'N', then IPIV is an output argument and on exit con‐
141 tains the pivot indices from the factorization A = P*L*U of the
142 original matrix A.
143 If FACT = 'E', then IPIV is an output argument and on exit con‐
145 tains the pivot indices from the factorization A = P*L*U of the
146 equilibrated matrix A.
147 EQUED (input or output) CHARACTER*1
149 Specifies the form of equilibration that was done. = 'N': No
150 equilibration (always true if FACT = 'N').
151 = 'R': Row equilibration, i.e., A has been premultiplied by
152 diag(R). = 'C': Column equilibration, i.e., A has been post‐
153 multiplied by diag(C). = 'B': Both row and column equilibra‐
154 tion, i.e., A has been replaced by diag(R) * A * diag(C).
155 EQUED is an input argument if FACT = 'F'; otherwise, it is an
156 output argument.
157 R (input or output) DOUBLE PRECISION array, dimension (N)
159 The row scale factors for A. If EQUED = 'R' or 'B', A is mul‐
160 tiplied on the left by diag(R); if EQUED = 'N' or 'C', R is not
161 accessed. R is an input argument if FACT = 'F'; otherwise, R
162 is an output argument. If FACT = 'F' and EQUED = 'R' or 'B',
163 each element of R must be positive.
164 C (input or output) DOUBLE PRECISION array, dimension (N)
166 The column scale factors for A. If EQUED = 'C' or 'B', A is
167 multiplied on the right by diag(C); if EQUED = 'N' or 'R', C is
168 not accessed. C is an input argument if FACT = 'F'; otherwise,
169 C is an output argument. If FACT = 'F' and EQUED = 'C' or 'B',
170 each element of C must be positive.
171 B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
173 On entry, the N-by-NRHS right hand side matrix B. On exit, if
174 EQUED = 'N', B is not modified; if TRANS = 'N' and EQUED = 'R'
175 or 'B', B is overwritten by diag(R)*B; if TRANS = 'T' or 'C'
176 and EQUED = 'C' or 'B', B is overwritten by diag(C)*B.
177 LDB (input) INTEGER
179 The leading dimension of the array B. LDB >= max(1,N).
180 X (output) COMPLEX*16 array, dimension (LDX,NRHS)
182 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
183 the original system of equations. Note that A and B are modi‐
184 fied on exit if EQUED .ne. 'N', and the solution to the equili‐
185 brated system is inv(diag(C))*X if TRANS = 'N' and EQUED = 'C'
186 or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R'
187 or 'B'.
188 LDX (input) INTEGER
190 The leading dimension of the array X. LDX >= max(1,N).
191 RCOND (output) DOUBLE PRECISION
193 The estimate of the reciprocal condition number of the matrix A
194 after equilibration (if done). If RCOND is less than the
195 machine precision (in particular, if RCOND = 0), the matrix is
196 singular to working precision. This condition is indicated by
197 a return code of INFO > 0.
198 FERR (output) DOUBLE PRECISION array, dimension (NRHS)
200 The estimated forward error bound for each solution vector X(j)
201 (the j-th column of the solution matrix X). If XTRUE is the
202 true solution corresponding to X(j), FERR(j) is an estimated
203 upper bound for the magnitude of the largest element in (X(j) -
204 XTRUE) divided by the magnitude of the largest element in X(j).
205 The estimate is as reliable as the estimate for RCOND, and is
206 almost always a slight overestimate of the true error.
207 BERR (output) DOUBLE PRECISION array, dimension (NRHS)
209 The componentwise relative backward error of each solution vec‐
210 tor X(j) (i.e., the smallest relative change in any element of
211 A or B that makes X(j) an exact solution).
212 WORK (workspace) COMPLEX*16 array, dimension (2*N)
214 RWORK (workspace/output) DOUBLE PRECISION array, dimension (2*N)
216 On exit, RWORK(1) contains the reciprocal pivot growth factor
217 norm(A)/norm(U). The "max absolute element" norm is used. If
218 RWORK(1) is much less than 1, then the stability of the LU fac‐
219 torization of the (equilibrated) matrix A could be poor. This
220 also means that the solution X, condition estimator RCOND, and
221 forward error bound FERR could be unreliable. If factorization
222 fails with 0<INFO<=N, then RWORK(1) contains the reciprocal
223 pivot growth factor for the leading INFO columns of A.
224 INFO (output) INTEGER
226 = 0: successful exit
227 < 0: if INFO = -i, the i-th argument had an illegal value
228 > 0: if INFO = i, and i is
229 <= N: U(i,i) is exactly zero. The factorization has been com‐
230 pleted, but the factor U is exactly singular, so the solution
231 and error bounds could not be computed. RCOND = 0 is returned.
232 = N+1: U is nonsingular, but RCOND is less than machine preci‐
233 sion, meaning that the matrix is singular to working precision.
234 Nevertheless, the solution and error bounds are computed
235 because there are a number of situations where the computed
236 solution can be more accurate than the value of RCOND would
237 suggest.
238 LAPACK driver routine (version 3.N1o)vember 2006 ZGESVX(1)