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