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