1SGESVXX(1) LAPACK driver routine (version 3.2) SGESVXX(1)
2
6 SGESVXX - SGESVXX use the LU factorization to compute the solution to a
7 real system of linear equations A * X = B, where A is an N-by-N matrix
8 and X and B are N-by-NRHS matrices
9
11 SUBROUTINE SGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
12 EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR,
13 N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS,
14 PARAMS, WORK, IWORK, INFO )
15 IMPLICIT NONE
17 CHARACTER EQUED, FACT, TRANS
19 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
21 N_ERR_BNDS
22 REAL RCOND, RPVGRW
24 INTEGER IPIV( * ), IWORK( * )
26 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ), X( LDX , *
28 ),WORK( * )
29 REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
31 ERR_BNDS_NORM( NRHS, * ), ERR_BNDS_COMP( NRHS, * )
32
34 SGESVXX uses the LU factorization to compute the solution to a
35 real system of linear equations A * X = B, where A is an
36 N-by-N matrix and X and B are N-by-NRHS matrices.
37 If requested, both normwise and maximum componentwise error bounds
38 are returned. SGESVXX will return a solution with a tiny
39 guaranteed error (O(eps) where eps is the working machine
40 precision) unless the matrix is very ill-conditioned, in which
41 case a warning is returned. Relevant condition numbers also are
42 calculated and returned.
43 SGESVXX accepts user-provided factorizations and equilibration
44 factors; see the definitions of the FACT and EQUED options.
45 Solving with refinement and using a factorization from a previous
46 SGESVXX call will also produce a solution with either O(eps)
47 errors or warnings, but we cannot make that claim for general
48 user-provided factorizations and equilibration factors if they
49 differ from what SGESVXX would itself produce.
50
52 The following steps are performed:
53 1. If FACT = 'E', real scaling factors are computed to equilibrate
54 the system:
55 TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
56 TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
57 TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
58 Whether or not the system will be equilibrated depends on the
59 scaling of the matrix A, but if equilibration is used, A is
60 overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
61 or diag(C)*B (if TRANS = 'T' or 'C').
62 2. If FACT = 'N' or 'E', the LU decomposition is used to factor
63 the matrix A (after equilibration if FACT = 'E') as
64 A = P * L * U,
65 where P is a permutation matrix, L is a unit lower triangular
66 matrix, and U is upper triangular.
67 3. If some U(i,i)=0, so that U is exactly singular, then the
68 routine returns with INFO = i. Otherwise, the factored form of A
69 is used to estimate the condition number of the matrix A (see
70 argument RCOND). If the reciprocal of the condition number is less
71 than machine precision, the routine still goes on to solve for X
72 and compute error bounds as described below.
73 4. The system of equations is solved for X using the factored form
74 of A.
75 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
76 the routine will use iterative refinement to try to get a small
77 error and error bounds. Refinement calculates the residual to at
78 least twice the working precision.
79 6. If equilibration was used, the matrix X is premultiplied by
80 diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
81 that it solves the original system before equilibration.
82
84 Some optional parameters are bundled in the PARAMS array. These set‐
85 tings determine how refinement is performed, but often the defaults are
86 acceptable. If the defaults are acceptable, users can pass NPARAMS = 0
87 which prevents the source code from accessing the PARAMS argument.
88 FACT (input) CHARACTER*1
90 Specifies whether or not the factored form of the matrix A is
91 supplied on entry, and if not, whether the matrix A should be
92 equilibrated before it is factored. = 'F': On entry, AF and
93 IPIV contain the factored form of A. If EQUED is not 'N', the
94 matrix A has been equilibrated with scaling factors given by R
95 and C. A, AF, and IPIV are not modified. = 'N': The matrix A
96 will be copied to AF and factored.
97 = 'E': The matrix A will be equilibrated if necessary, then
98 copied to AF and factored.
99 TRANS (input) CHARACTER*1
101 Specifies the form of the system of equations:
102 = 'N': A * X = B (No transpose)
103 = 'T': A**T * X = B (Transpose)
104 = 'C': A**H * X = B (Conjugate Transpose = Transpose)
105 N (input) INTEGER
107 The number of linear equations, i.e., the order of the matrix
108 A. N >= 0.
109 NRHS (input) INTEGER
111 The number of right hand sides, i.e., the number of columns of
112 the matrices B and X. NRHS >= 0.
113 A (input/output) REAL array, dimension (LDA,N)
115 On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is not
116 'N', then A must have been equilibrated by the scaling factors
117 in R and/or C. A is not modified if FACT = 'F' or 'N', or if
118 FACT = 'E' and EQUED = 'N' on exit. On exit, if EQUED .ne.
119 'N', A is scaled as follows: EQUED = 'R': A := diag(R) * A
120 EQUED = 'C': A := A * diag(C)
121 EQUED = 'B': A := diag(R) * A * diag(C).
122 LDA (input) INTEGER
124 The leading dimension of the array A. LDA >= max(1,N).
125 AF (input or output) REAL array, dimension (LDAF,N)
127 If FACT = 'F', then AF is an input argument and on entry con‐
128 tains the factors L and U from the factorization A = P*L*U as
129 computed by SGETRF. If EQUED .ne. 'N', then AF is the factored
130 form of the equilibrated matrix A. If FACT = 'N', then AF is
131 an output argument and on exit returns the factors L and U from
132 the factorization A = P*L*U of the original matrix A. If FACT
133 = 'E', then AF is an output argument and on exit returns the
134 factors L and U from the factorization A = P*L*U of the equili‐
135 brated matrix A (see the description of A for the form of the
136 equilibrated matrix).
137 LDAF (input) INTEGER
139 The leading dimension of the array AF. LDAF >= max(1,N).
140 IPIV (input or output) INTEGER array, dimension (N)
142 If FACT = 'F', then IPIV is an input argument and on entry con‐
143 tains the pivot indices from the factorization A = P*L*U as
144 computed by SGETRF; row i of the matrix was interchanged with
145 row IPIV(i). If FACT = 'N', then IPIV is an output argument
146 and on exit contains the pivot indices from the factorization A
147 = P*L*U of the original matrix A. If FACT = 'E', then IPIV is
148 an output argument and on exit contains the pivot indices from
149 the factorization A = P*L*U of the equilibrated matrix A.
150 EQUED (input or output) CHARACTER*1
152 Specifies the form of equilibration that was done. = 'N': No
153 equilibration (always true if FACT = 'N').
154 = 'R': Row equilibration, i.e., A has been premultiplied by
155 diag(R). = 'C': Column equilibration, i.e., A has been post‐
156 multiplied by diag(C). = 'B': Both row and column equilibra‐
157 tion, i.e., A has been replaced by diag(R) * A * diag(C).
158 EQUED is an input argument if FACT = 'F'; otherwise, it is an
159 output argument.
160 R (input or output) REAL array, dimension (N)
162 The row scale factors for A. If EQUED = 'R' or 'B', A is mul‐
163 tiplied on the left by diag(R); if EQUED = 'N' or 'C', R is not
164 accessed. R is an input argument if FACT = 'F'; otherwise, R
165 is an output argument. If FACT = 'F' and EQUED = 'R' or 'B',
166 each element of R must be positive. If R is output, each ele‐
167 ment of R is a power of the radix. If R is input, each element
168 of R should be a power of the radix to ensure a reliable solu‐
169 tion and error estimates. Scaling by powers of the radix does
170 not cause rounding errors unless the result underflows or over‐
171 flows. Rounding errors during scaling lead to refining with a
172 matrix that is not equivalent to the input matrix, producing
173 error estimates that may not be reliable.
174 C (input or output) REAL array, dimension (N)
176 The column scale factors for A. If EQUED = 'C' or 'B', A is
177 multiplied on the right by diag(C); if EQUED = 'N' or 'R', C is
178 not accessed. C is an input argument if FACT = 'F'; otherwise,
179 C is an output argument. If FACT = 'F' and EQUED = 'C' or 'B',
180 each element of C must be positive. If C is output, each ele‐
181 ment of C is a power of the radix. If C is input, each element
182 of C should be a power of the radix to ensure a reliable solu‐
183 tion and error estimates. Scaling by powers of the radix does
184 not cause rounding errors unless the result underflows or over‐
185 flows. Rounding errors during scaling lead to refining with a
186 matrix that is not equivalent to the input matrix, producing
187 error estimates that may not be reliable.
188 B (input/output) REAL array, dimension (LDB,NRHS)
190 On entry, the N-by-NRHS right hand side matrix B. On exit, if
191 EQUED = 'N', B is not modified; if TRANS = 'N' and EQUED = 'R'
192 or 'B', B is overwritten by diag(R)*B; if TRANS = 'T' or 'C'
193 and EQUED = 'C' or 'B', B is overwritten by diag(C)*B.
194 LDB (input) INTEGER
196 The leading dimension of the array B. LDB >= max(1,N).
197 X (output) REAL array, dimension (LDX,NRHS)
199 If INFO = 0, the N-by-NRHS solution matrix X to the original
200 system of equations. Note that A and B are modified on exit if
201 EQUED .ne. 'N', and the solution to the equilibrated system is
202 inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
203 inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
204 LDX (input) INTEGER
206 The leading dimension of the array X. LDX >= max(1,N).
207 RCOND (output) REAL
209 Reciprocal scaled condition number. This is an estimate of the
210 reciprocal Skeel condition number of the matrix A after equili‐
211 bration (if done). If this is less than the machine precision
212 (in particular, if it is zero), the matrix is singular to work‐
213 ing precision. Note that the error may still be small even if
214 this number is very small and the matrix appears ill- condi‐
215 tioned.
216 RPVGRW (output) REAL
218 Reciprocal pivot growth. On exit, this contains the reciprocal
219 pivot growth factor norm(A)/norm(U). The "max absolute element"
220 norm is used. If this is much less than 1, then the stability
221 of the LU factorization of the (equilibrated) matrix A could be
222 poor. This also means that the solution X, estimated condition
223 numbers, and error bounds could be unreliable. If factorization
224 fails with 0<INFO<=N, then this contains the reciprocal pivot
225 growth factor for the leading INFO columns of A. In SGESVX,
226 this quantity is returned in WORK(1).
227 BERR (output) REAL array, dimension (NRHS)
229 Componentwise relative backward error. This is the component‐
230 wise relative backward error of each solution vector X(j)
231 (i.e., the smallest relative change in any element of A or B
232 that makes X(j) an exact solution). N_ERR_BNDS (input) INTEGER
233 Number of error bounds to return for each right hand side and
234 each type (normwise or componentwise). See ERR_BNDS_NORM and
235 ERR_BNDS_COMP below.
236 ERR_BNDS_NORM (output) REAL array, dimension (NRHS, N_ERR_BNDS)
238 For each right-hand side, this array contains informa‐
239 tion about various error bounds and condition numbers
240 corresponding to the normwise relative error, which is
241 defined as follows: Normwise relative error in the ith
242 solution vector: max_j (abs(XTRUE(j,i) - X(j,i)))
243 ------------------------------ max_j abs(X(j,i)) The
244 array is indexed by the type of error information as
245 described below. There currently are up to three pieces
246 of information returned. The first index in
247 ERR_BNDS_NORM(i,:) corresponds to the ith right-hand
248 side. The second index in ERR_BNDS_NORM(:,err) contains
249 the following three fields: err = 1 "Trust/don't trust"
250 boolean. Trust the answer if the reciprocal condition
251 number is less than the threshold sqrt(n) *
252 slamch('Epsilon'). err = 2 "Guaranteed" error bound:
253 The estimated forward error, almost certainly within a
254 factor of 10 of the true error so long as the next entry
255 is greater than the threshold sqrt(n) *
256 slamch('Epsilon'). This error bound should only be
257 trusted if the previous boolean is true. err = 3
258 Reciprocal condition number: Estimated normwise recipro‐
259 cal condition number. Compared with the threshold
260 sqrt(n) * slamch('Epsilon') to determine if the error
261 estimate is "guaranteed". These reciprocal condition
262 numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for
263 some appropriately scaled matrix Z. Let Z = S*A, where
264 S scales each row by a power of the radix so all abso‐
265 lute row sums of Z are approximately 1. See Lapack
266 Working Note 165 for further details and extra cautions.
267 ERR_BNDS_COMP (output) REAL array, dimension (NRHS, N_ERR_BNDS)
269 For each right-hand side, this array contains informa‐
270 tion about various error bounds and condition numbers
271 corresponding to the componentwise relative error, which
272 is defined as follows: Componentwise relative error in
273 the ith solution vector: abs(XTRUE(j,i) - X(j,i)) max_j
274 ---------------------- abs(X(j,i)) The array is indexed
275 by the right-hand side i (on which the componentwise
276 relative error depends), and the type of error informa‐
277 tion as described below. There currently are up to three
278 pieces of information returned for each right-hand side.
279 If componentwise accuracy is not requested (PARAMS(3) =
280 0.0), then ERR_BNDS_COMP is not accessed. If N_ERR_BNDS
281 .LT. 3, then at most the first (:,N_ERR_BNDS) entries
282 are returned. The first index in ERR_BNDS_COMP(i,:)
283 corresponds to the ith right-hand side. The second
284 index in ERR_BNDS_COMP(:,err) contains the following
285 three fields: err = 1 "Trust/don't trust" boolean. Trust
286 the answer if the reciprocal condition number is less
287 than the threshold sqrt(n) * slamch('Epsilon'). err = 2
288 "Guaranteed" error bound: The estimated forward error,
289 almost certainly within a factor of 10 of the true error
290 so long as the next entry is greater than the threshold
291 sqrt(n) * slamch('Epsilon'). This error bound should
292 only be trusted if the previous boolean is true. err =
293 3 Reciprocal condition number: Estimated componentwise
294 reciprocal condition number. Compared with the thresh‐
295 old sqrt(n) * slamch('Epsilon') to determine if the
296 error estimate is "guaranteed". These reciprocal condi‐
297 tion numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf))
298 for some appropriately scaled matrix Z. Let Z =
299 S*(A*diag(x)), where x is the solution for the current
300 right-hand side and S scales each row of A*diag(x) by a
301 power of the radix so all absolute row sums of Z are
302 approximately 1. See Lapack Working Note 165 for fur‐
303 ther details and extra cautions. NPARAMS (input) INTE‐
304 GER Specifies the number of parameters set in PARAMS.
305 If .LE. 0, the PARAMS array is never referenced and
306 default values are used.
307 PARAMS (input / output) REAL array, dimension NPARAMS
309 Specifies algorithm parameters. If an entry is .LT. 0.0, then
310 that entry will be filled with default value used for that
311 parameter. Only positions up to NPARAMS are accessed; defaults
312 are used for higher-numbered parameters.
313 PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
314 refinement or not. Default: 1.0
315 = 0.0 : No refinement is performed, and no error bounds are
316 computed. = 1.0 : Use the double-precision refinement algo‐
317 rithm, possibly with doubled-single computations if the compi‐
318 lation environment does not support DOUBLE PRECISION. (other
319 values are reserved for future use) PARAMS(LA_LINRX_ITHRESH_I =
320 2) : Maximum number of residual computations allowed for
321 refinement. Default: 10
322 Aggressive: Set to 100 to permit convergence using approximate
323 factorizations or factorizations other than LU. If the factor‐
324 ization uses a technique other than Gaussian elimination, the
325 guarantees in err_bnds_norm and err_bnds_comp may no longer be
326 trustworthy. PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining
327 if the code will attempt to find a solution with small compo‐
328 nentwise relative error in the double-precision algorithm.
329 Positive is true, 0.0 is false. Default: 1.0 (attempt compo‐
330 nentwise convergence)
331 WORK (workspace) REAL array, dimension (4*N)
333 IWORK (workspace) INTEGER array, dimension (N)
335 INFO (output) INTEGER
337 = 0: Successful exit. The solution to every right-hand side is
338 guaranteed. < 0: If INFO = -i, the i-th argument had an ille‐
339 gal value
340 > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
341 has been completed, but the factor U is exactly singular, so
342 the solution and error bounds could not be computed. RCOND = 0
343 is returned. = N+J: The solution corresponding to the Jth
344 right-hand side is not guaranteed. The solutions corresponding
345 to other right- hand sides K with K > J may not be guaranteed
346 as well, but only the first such right-hand side is reported.
347 If a small componentwise error is not requested (PARAMS(3) =
348 0.0) then the Jth right-hand side is the first with a normwise
349 error bound that is not guaranteed (the smallest J such that
350 ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) the Jth
351 right-hand side is the first with either a normwise or compo‐
352 nentwise error bound that is not guaranteed (the smallest J
353 such that either ERR_BNDS_NORM(J,1) = 0.0 or ERR_BNDS_COMP(J,1)
354 = 0.0). See the definition of ERR_BNDS_NORM(:,1) and
355 ERR_BNDS_COMP(:,1). To get information about all of the right-
356 hand sides check ERR_BNDS_NORM or ERR_BNDS_COMP.
357 LAPACK driver routine (versionNo3v.e2m)ber 2008 SGESVXX(1)