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