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