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