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