LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dgesvxx.f
Go to the documentation of this file.
1 *> \brief <b> DGESVXX computes the solution to system of linear equations A * X = B for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DGESVXX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvxx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvxx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvxx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
22 * EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW,
23 * BERR, N_ERR_BNDS, ERR_BNDS_NORM,
24 * ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK,
25 * INFO )
26 *
27 * .. Scalar Arguments ..
28 * CHARACTER EQUED, FACT, TRANS
29 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
30 * $ N_ERR_BNDS
31 * DOUBLE PRECISION RCOND, RPVGRW
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IPIV( * ), IWORK( * )
35 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
36 * $ X( LDX , * ),WORK( * )
37 * DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
38 * $ ERR_BNDS_NORM( NRHS, * ),
39 * $ ERR_BNDS_COMP( NRHS, * )
40 * ..
41 *
42 *
43 *> \par Purpose:
44 * =============
45 *>
46 *> \verbatim
47 *>
48 *> DGESVXX uses the LU factorization to compute the solution to a
49 *> double precision system of linear equations A * X = B, where A is an
50 *> N-by-N matrix and X and B are N-by-NRHS matrices.
51 *>
52 *> If requested, both normwise and maximum componentwise error bounds
53 *> are returned. DGESVXX will return a solution with a tiny
54 *> guaranteed error (O(eps) where eps is the working machine
55 *> precision) unless the matrix is very ill-conditioned, in which
56 *> case a warning is returned. Relevant condition numbers also are
57 *> calculated and returned.
58 *>
59 *> DGESVXX accepts user-provided factorizations and equilibration
60 *> factors; see the definitions of the FACT and EQUED options.
61 *> Solving with refinement and using a factorization from a previous
62 *> DGESVXX call will also produce a solution with either O(eps)
63 *> errors or warnings, but we cannot make that claim for general
64 *> user-provided factorizations and equilibration factors if they
65 *> differ from what DGESVXX would itself produce.
66 *> \endverbatim
67 *
68 *> \par Description:
69 * =================
70 *>
71 *> \verbatim
72 *>
73 *> The following steps are performed:
74 *>
75 *> 1. If FACT = 'E', double precision scaling factors are computed to equilibrate
76 *> the system:
77 *>
78 *> TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
79 *> TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
80 *> TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
81 *>
82 *> Whether or not the system will be equilibrated depends on the
83 *> scaling of the matrix A, but if equilibration is used, A is
84 *> overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
85 *> or diag(C)*B (if TRANS = 'T' or 'C').
86 *>
87 *> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor
88 *> the matrix A (after equilibration if FACT = 'E') as
89 *>
90 *> A = P * L * U,
91 *>
92 *> where P is a permutation matrix, L is a unit lower triangular
93 *> matrix, and U is upper triangular.
94 *>
95 *> 3. If some U(i,i)=0, so that U is exactly singular, then the
96 *> routine returns with INFO = i. Otherwise, the factored form of A
97 *> is used to estimate the condition number of the matrix A (see
98 *> argument RCOND). If the reciprocal of the condition number is less
99 *> than machine precision, the routine still goes on to solve for X
100 *> and compute error bounds as described below.
101 *>
102 *> 4. The system of equations is solved for X using the factored form
103 *> of A.
104 *>
105 *> 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
106 *> the routine will use iterative refinement to try to get a small
107 *> error and error bounds. Refinement calculates the residual to at
108 *> least twice the working precision.
109 *>
110 *> 6. If equilibration was used, the matrix X is premultiplied by
111 *> diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
112 *> that it solves the original system before equilibration.
113 *> \endverbatim
114 *
115 * Arguments:
116 * ==========
117 *
118 *> \verbatim
119 *> Some optional parameters are bundled in the PARAMS array. These
120 *> settings determine how refinement is performed, but often the
121 *> defaults are acceptable. If the defaults are acceptable, users
122 *> can pass NPARAMS = 0 which prevents the source code from accessing
123 *> the PARAMS argument.
124 *> \endverbatim
125 *>
126 *> \param[in] FACT
127 *> \verbatim
128 *> FACT is CHARACTER*1
129 *> Specifies whether or not the factored form of the matrix A is
130 *> supplied on entry, and if not, whether the matrix A should be
131 *> equilibrated before it is factored.
132 *> = 'F': On entry, AF and IPIV contain the factored form of A.
133 *> If EQUED is not 'N', the matrix A has been
134 *> equilibrated with scaling factors given by R and C.
135 *> A, AF, and IPIV are not modified.
136 *> = 'N': The matrix A will be copied to AF and factored.
137 *> = 'E': The matrix A will be equilibrated if necessary, then
138 *> copied to AF and factored.
139 *> \endverbatim
140 *>
141 *> \param[in] TRANS
142 *> \verbatim
143 *> TRANS is CHARACTER*1
144 *> Specifies the form of the system of equations:
145 *> = 'N': A * X = B (No transpose)
146 *> = 'T': A**T * X = B (Transpose)
147 *> = 'C': A**H * X = B (Conjugate Transpose = Transpose)
148 *> \endverbatim
149 *>
150 *> \param[in] N
151 *> \verbatim
152 *> N is INTEGER
153 *> The number of linear equations, i.e., the order of the
154 *> matrix A. N >= 0.
155 *> \endverbatim
156 *>
157 *> \param[in] NRHS
158 *> \verbatim
159 *> NRHS is INTEGER
160 *> The number of right hand sides, i.e., the number of columns
161 *> of the matrices B and X. NRHS >= 0.
162 *> \endverbatim
163 *>
164 *> \param[in,out] A
165 *> \verbatim
166 *> A is DOUBLE PRECISION array, dimension (LDA,N)
167 *> On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is
168 *> not 'N', then A must have been equilibrated by the scaling
169 *> factors in R and/or C. A is not modified if FACT = 'F' or
170 *> 'N', or if FACT = 'E' and EQUED = 'N' on exit.
171 *>
172 *> On exit, if EQUED .ne. 'N', A is scaled as follows:
173 *> EQUED = 'R': A := diag(R) * A
174 *> EQUED = 'C': A := A * diag(C)
175 *> EQUED = 'B': A := diag(R) * A * diag(C).
176 *> \endverbatim
177 *>
178 *> \param[in] LDA
179 *> \verbatim
180 *> LDA is INTEGER
181 *> The leading dimension of the array A. LDA >= max(1,N).
182 *> \endverbatim
183 *>
184 *> \param[in,out] AF
185 *> \verbatim
186 *> AF is DOUBLE PRECISION array, dimension (LDAF,N)
187 *> If FACT = 'F', then AF is an input argument and on entry
188 *> contains the factors L and U from the factorization
189 *> A = P*L*U as computed by DGETRF. If EQUED .ne. 'N', then
190 *> AF is the factored form of the equilibrated matrix A.
191 *>
192 *> If FACT = 'N', then AF is an output argument and on exit
193 *> returns the factors L and U from the factorization A = P*L*U
194 *> of the original matrix A.
195 *>
196 *> If FACT = 'E', then AF is an output argument and on exit
197 *> returns the factors L and U from the factorization A = P*L*U
198 *> of the equilibrated matrix A (see the description of A for
199 *> the form of the equilibrated matrix).
200 *> \endverbatim
201 *>
202 *> \param[in] LDAF
203 *> \verbatim
204 *> LDAF is INTEGER
205 *> The leading dimension of the array AF. LDAF >= max(1,N).
206 *> \endverbatim
207 *>
208 *> \param[in,out] IPIV
209 *> \verbatim
210 *> IPIV is INTEGER array, dimension (N)
211 *> If FACT = 'F', then IPIV is an input argument and on entry
212 *> contains the pivot indices from the factorization A = P*L*U
213 *> as computed by DGETRF; row i of the matrix was interchanged
214 *> with row IPIV(i).
215 *>
216 *> If FACT = 'N', then IPIV is an output argument and on exit
217 *> contains the pivot indices from the factorization A = P*L*U
218 *> of the original matrix A.
219 *>
220 *> If FACT = 'E', then IPIV is an output argument and on exit
221 *> contains the pivot indices from the factorization A = P*L*U
222 *> of the equilibrated matrix A.
223 *> \endverbatim
224 *>
225 *> \param[in,out] EQUED
226 *> \verbatim
227 *> EQUED is CHARACTER*1
228 *> Specifies the form of equilibration that was done.
229 *> = 'N': No equilibration (always true if FACT = 'N').
230 *> = 'R': Row equilibration, i.e., A has been premultiplied by
231 *> diag(R).
232 *> = 'C': Column equilibration, i.e., A has been postmultiplied
233 *> by diag(C).
234 *> = 'B': Both row and column equilibration, i.e., A has been
235 *> replaced by diag(R) * A * diag(C).
236 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
237 *> output argument.
238 *> \endverbatim
239 *>
240 *> \param[in,out] R
241 *> \verbatim
242 *> R is DOUBLE PRECISION array, dimension (N)
243 *> The row scale factors for A. If EQUED = 'R' or 'B', A is
244 *> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
245 *> is not accessed. R is an input argument if FACT = 'F';
246 *> otherwise, R is an output argument. If FACT = 'F' and
247 *> EQUED = 'R' or 'B', each element of R must be positive.
248 *> If R is output, each element of R is a power of the radix.
249 *> If R is input, each element of R should be a power of the radix
250 *> to ensure a reliable solution and error estimates. Scaling by
251 *> powers of the radix does not cause rounding errors unless the
252 *> result underflows or overflows. Rounding errors during scaling
253 *> lead to refining with a matrix that is not equivalent to the
254 *> input matrix, producing error estimates that may not be
255 *> reliable.
256 *> \endverbatim
257 *>
258 *> \param[in,out] C
259 *> \verbatim
260 *> C is DOUBLE PRECISION array, dimension (N)
261 *> The column scale factors for A. If EQUED = 'C' or 'B', A is
262 *> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
263 *> is not accessed. C is an input argument if FACT = 'F';
264 *> otherwise, C is an output argument. If FACT = 'F' and
265 *> EQUED = 'C' or 'B', each element of C must be positive.
266 *> If C is output, each element of C is a power of the radix.
267 *> If C is input, each element of C should be a power of the radix
268 *> to ensure a reliable solution and error estimates. Scaling by
269 *> powers of the radix does not cause rounding errors unless the
270 *> result underflows or overflows. Rounding errors during scaling
271 *> lead to refining with a matrix that is not equivalent to the
272 *> input matrix, producing error estimates that may not be
273 *> reliable.
274 *> \endverbatim
275 *>
276 *> \param[in,out] B
277 *> \verbatim
278 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
279 *> On entry, the N-by-NRHS right hand side matrix B.
280 *> On exit,
281 *> if EQUED = 'N', B is not modified;
282 *> if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
283 *> diag(R)*B;
284 *> if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
285 *> overwritten by diag(C)*B.
286 *> \endverbatim
287 *>
288 *> \param[in] LDB
289 *> \verbatim
290 *> LDB is INTEGER
291 *> The leading dimension of the array B. LDB >= max(1,N).
292 *> \endverbatim
293 *>
294 *> \param[out] X
295 *> \verbatim
296 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
297 *> If INFO = 0, the N-by-NRHS solution matrix X to the original
298 *> system of equations. Note that A and B are modified on exit
299 *> if EQUED .ne. 'N', and the solution to the equilibrated system is
300 *> inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
301 *> inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
302 *> \endverbatim
303 *>
304 *> \param[in] LDX
305 *> \verbatim
306 *> LDX is INTEGER
307 *> The leading dimension of the array X. LDX >= max(1,N).
308 *> \endverbatim
309 *>
310 *> \param[out] RCOND
311 *> \verbatim
312 *> RCOND is DOUBLE PRECISION
313 *> Reciprocal scaled condition number. This is an estimate of the
314 *> reciprocal Skeel condition number of the matrix A after
315 *> equilibration (if done). If this is less than the machine
316 *> precision (in particular, if it is zero), the matrix is singular
317 *> to working precision. Note that the error may still be small even
318 *> if this number is very small and the matrix appears ill-
319 *> conditioned.
320 *> \endverbatim
321 *>
322 *> \param[out] RPVGRW
323 *> \verbatim
324 *> RPVGRW is DOUBLE PRECISION
325 *> Reciprocal pivot growth. On exit, this contains the reciprocal
326 *> pivot growth factor norm(A)/norm(U). The "max absolute element"
327 *> norm is used. If this is much less than 1, then the stability of
328 *> the LU factorization of the (equilibrated) matrix A could be poor.
329 *> This also means that the solution X, estimated condition numbers,
330 *> and error bounds could be unreliable. If factorization fails with
331 *> 0<INFO<=N, then this contains the reciprocal pivot growth factor
332 *> for the leading INFO columns of A. In DGESVX, this quantity is
333 *> returned in WORK(1).
334 *> \endverbatim
335 *>
336 *> \param[out] BERR
337 *> \verbatim
338 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
339 *> Componentwise relative backward error. This is the
340 *> componentwise relative backward error of each solution vector X(j)
341 *> (i.e., the smallest relative change in any element of A or B that
342 *> makes X(j) an exact solution).
343 *> \endverbatim
344 *>
345 *> \param[in] N_ERR_BNDS
346 *> \verbatim
347 *> N_ERR_BNDS is INTEGER
348 *> Number of error bounds to return for each right hand side
349 *> and each type (normwise or componentwise). See ERR_BNDS_NORM and
350 *> ERR_BNDS_COMP below.
351 *> \endverbatim
352 *>
353 *> \param[out] ERR_BNDS_NORM
354 *> \verbatim
355 *> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
356 *> For each right-hand side, this array contains information about
357 *> various error bounds and condition numbers corresponding to the
358 *> normwise relative error, which is defined as follows:
359 *>
360 *> Normwise relative error in the ith solution vector:
361 *> max_j (abs(XTRUE(j,i) - X(j,i)))
362 *> ------------------------------
363 *> max_j abs(X(j,i))
364 *>
365 *> The array is indexed by the type of error information as described
366 *> below. There currently are up to three pieces of information
367 *> returned.
368 *>
369 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
370 *> right-hand side.
371 *>
372 *> The second index in ERR_BNDS_NORM(:,err) contains the following
373 *> three fields:
374 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
375 *> reciprocal condition number is less than the threshold
376 *> sqrt(n) * dlamch('Epsilon').
377 *>
378 *> err = 2 "Guaranteed" error bound: The estimated forward error,
379 *> almost certainly within a factor of 10 of the true error
380 *> so long as the next entry is greater than the threshold
381 *> sqrt(n) * dlamch('Epsilon'). This error bound should only
382 *> be trusted if the previous boolean is true.
383 *>
384 *> err = 3 Reciprocal condition number: Estimated normwise
385 *> reciprocal condition number. Compared with the threshold
386 *> sqrt(n) * dlamch('Epsilon') to determine if the error
387 *> estimate is "guaranteed". These reciprocal condition
388 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
389 *> appropriately scaled matrix Z.
390 *> Let Z = S*A, where S scales each row by a power of the
391 *> radix so all absolute row sums of Z are approximately 1.
392 *>
393 *> See Lapack Working Note 165 for further details and extra
394 *> cautions.
395 *> \endverbatim
396 *>
397 *> \param[out] ERR_BNDS_COMP
398 *> \verbatim
399 *> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
400 *> For each right-hand side, this array contains information about
401 *> various error bounds and condition numbers corresponding to the
402 *> componentwise relative error, which is defined as follows:
403 *>
404 *> Componentwise relative error in the ith solution vector:
405 *> abs(XTRUE(j,i) - X(j,i))
406 *> max_j ----------------------
407 *> abs(X(j,i))
408 *>
409 *> The array is indexed by the right-hand side i (on which the
410 *> componentwise relative error depends), and the type of error
411 *> information as described below. There currently are up to three
412 *> pieces of information returned for each right-hand side. If
413 *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
414 *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
415 *> the first (:,N_ERR_BNDS) entries are returned.
416 *>
417 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
418 *> right-hand side.
419 *>
420 *> The second index in ERR_BNDS_COMP(:,err) contains the following
421 *> three fields:
422 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
423 *> reciprocal condition number is less than the threshold
424 *> sqrt(n) * dlamch('Epsilon').
425 *>
426 *> err = 2 "Guaranteed" error bound: The estimated forward error,
427 *> almost certainly within a factor of 10 of the true error
428 *> so long as the next entry is greater than the threshold
429 *> sqrt(n) * dlamch('Epsilon'). This error bound should only
430 *> be trusted if the previous boolean is true.
431 *>
432 *> err = 3 Reciprocal condition number: Estimated componentwise
433 *> reciprocal condition number. Compared with the threshold
434 *> sqrt(n) * dlamch('Epsilon') to determine if the error
435 *> estimate is "guaranteed". These reciprocal condition
436 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
437 *> appropriately scaled matrix Z.
438 *> Let Z = S*(A*diag(x)), where x is the solution for the
439 *> current right-hand side and S scales each row of
440 *> A*diag(x) by a power of the radix so all absolute row
441 *> sums of Z are approximately 1.
442 *>
443 *> See Lapack Working Note 165 for further details and extra
444 *> cautions.
445 *> \endverbatim
446 *>
447 *> \param[in] NPARAMS
448 *> \verbatim
449 *> NPARAMS is INTEGER
450 *> Specifies the number of parameters set in PARAMS. If .LE. 0, the
451 *> PARAMS array is never referenced and default values are used.
452 *> \endverbatim
453 *>
454 *> \param[in,out] PARAMS
455 *> \verbatim
456 *> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
457 *> Specifies algorithm parameters. If an entry is .LT. 0.0, then
458 *> that entry will be filled with default value used for that
459 *> parameter. Only positions up to NPARAMS are accessed; defaults
460 *> are used for higher-numbered parameters.
461 *>
462 *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
463 *> refinement or not.
464 *> Default: 1.0D+0
465 *> = 0.0 : No refinement is performed, and no error bounds are
466 *> computed.
467 *> = 1.0 : Use the extra-precise refinement algorithm.
468 *> (other values are reserved for future use)
469 *>
470 *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
471 *> computations allowed for refinement.
472 *> Default: 10
473 *> Aggressive: Set to 100 to permit convergence using approximate
474 *> factorizations or factorizations other than LU. If
475 *> the factorization uses a technique other than
476 *> Gaussian elimination, the guarantees in
477 *> err_bnds_norm and err_bnds_comp may no longer be
478 *> trustworthy.
479 *>
480 *> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
481 *> will attempt to find a solution with small componentwise
482 *> relative error in the double-precision algorithm. Positive
483 *> is true, 0.0 is false.
484 *> Default: 1.0 (attempt componentwise convergence)
485 *> \endverbatim
486 *>
487 *> \param[out] WORK
488 *> \verbatim
489 *> WORK is DOUBLE PRECISION array, dimension (4*N)
490 *> \endverbatim
491 *>
492 *> \param[out] IWORK
493 *> \verbatim
494 *> IWORK is INTEGER array, dimension (N)
495 *> \endverbatim
496 *>
497 *> \param[out] INFO
498 *> \verbatim
499 *> INFO is INTEGER
500 *> = 0: Successful exit. The solution to every right-hand side is
501 *> guaranteed.
502 *> < 0: If INFO = -i, the i-th argument had an illegal value
503 *> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
504 *> has been completed, but the factor U is exactly singular, so
505 *> the solution and error bounds could not be computed. RCOND = 0
506 *> is returned.
507 *> = N+J: The solution corresponding to the Jth right-hand side is
508 *> not guaranteed. The solutions corresponding to other right-
509 *> hand sides K with K > J may not be guaranteed as well, but
510 *> only the first such right-hand side is reported. If a small
511 *> componentwise error is not requested (PARAMS(3) = 0.0) then
512 *> the Jth right-hand side is the first with a normwise error
513 *> bound that is not guaranteed (the smallest J such
514 *> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
515 *> the Jth right-hand side is the first with either a normwise or
516 *> componentwise error bound that is not guaranteed (the smallest
517 *> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
518 *> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
519 *> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
520 *> about all of the right-hand sides check ERR_BNDS_NORM or
521 *> ERR_BNDS_COMP.
522 *> \endverbatim
523 *
524 * Authors:
525 * ========
526 *
527 *> \author Univ. of Tennessee
528 *> \author Univ. of California Berkeley
529 *> \author Univ. of Colorado Denver
530 *> \author NAG Ltd.
531 *
532 *> \date April 2012
533 *
534 *> \ingroup doubleGEsolve
535 *
536 * =====================================================================
537  SUBROUTINE dgesvxx( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
538  $ equed, r, c, b, ldb, x, ldx, rcond, rpvgrw,
539  $ berr, n_err_bnds, err_bnds_norm,
540  $ err_bnds_comp, nparams, params, work, iwork,
541  $ info )
542 *
543 * -- LAPACK driver routine (version 3.4.1) --
544 * -- LAPACK is a software package provided by Univ. of Tennessee, --
545 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
546 * April 2012
547 *
548 * .. Scalar Arguments ..
549  CHARACTER equed, fact, trans
550  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs, nparams,
551  $ n_err_bnds
552  DOUBLE PRECISION rcond, rpvgrw
553 * ..
554 * .. Array Arguments ..
555  INTEGER ipiv( * ), iwork( * )
556  DOUBLE PRECISION a( lda, * ), af( ldaf, * ), b( ldb, * ),
557  $ x( ldx , * ),work( * )
558  DOUBLE PRECISION r( * ), c( * ), params( * ), berr( * ),
559  $ err_bnds_norm( nrhs, * ),
560  $ err_bnds_comp( nrhs, * )
561 * ..
562 *
563 * =====================================================================
564 *
565 * .. Parameters ..
566  DOUBLE PRECISION zero, one
567  parameter( zero = 0.0d+0, one = 1.0d+0 )
568  INTEGER final_nrm_err_i, final_cmp_err_i, berr_i
569  INTEGER rcond_i, nrm_rcond_i, nrm_err_i, cmp_rcond_i
570  INTEGER cmp_err_i, piv_growth_i
571  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
572  $ berr_i = 3 )
573  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
574  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
575  $ piv_growth_i = 9 )
576 * ..
577 * .. Local Scalars ..
578  LOGICAL colequ, equil, nofact, notran, rowequ
579  INTEGER infequ, j
580  DOUBLE PRECISION amax, bignum, colcnd, rcmax, rcmin, rowcnd,
581  $ smlnum
582 * ..
583 * .. External Functions ..
584  EXTERNAL lsame, dlamch, dla_gerpvgrw
585  LOGICAL lsame
586  DOUBLE PRECISION dlamch, dla_gerpvgrw
587 * ..
588 * .. External Subroutines ..
589  EXTERNAL dgeequb, dgetrf, dgetrs, dlacpy, dlaqge,
591 * ..
592 * .. Intrinsic Functions ..
593  INTRINSIC max, min
594 * ..
595 * .. Executable Statements ..
596 *
597  info = 0
598  nofact = lsame( fact, 'N' )
599  equil = lsame( fact, 'E' )
600  notran = lsame( trans, 'N' )
601  smlnum = dlamch( 'Safe minimum' )
602  bignum = one / smlnum
603  IF( nofact .OR. equil ) THEN
604  equed = 'N'
605  rowequ = .false.
606  colequ = .false.
607  ELSE
608  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
609  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
610  END IF
611 *
612 * Default is failure. If an input parameter is wrong or
613 * factorization fails, make everything look horrible. Only the
614 * pivot growth is set here, the rest is initialized in DGERFSX.
615 *
616  rpvgrw = zero
617 *
618 * Test the input parameters. PARAMS is not tested until DGERFSX.
619 *
620  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
621  $ lsame( fact, 'F' ) ) THEN
622  info = -1
623  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
624  $ lsame( trans, 'C' ) ) THEN
625  info = -2
626  ELSE IF( n.LT.0 ) THEN
627  info = -3
628  ELSE IF( nrhs.LT.0 ) THEN
629  info = -4
630  ELSE IF( lda.LT.max( 1, n ) ) THEN
631  info = -6
632  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
633  info = -8
634  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
635  $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
636  info = -10
637  ELSE
638  IF( rowequ ) THEN
639  rcmin = bignum
640  rcmax = zero
641  DO 10 j = 1, n
642  rcmin = min( rcmin, r( j ) )
643  rcmax = max( rcmax, r( j ) )
644  10 continue
645  IF( rcmin.LE.zero ) THEN
646  info = -11
647  ELSE IF( n.GT.0 ) THEN
648  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
649  ELSE
650  rowcnd = one
651  END IF
652  END IF
653  IF( colequ .AND. info.EQ.0 ) THEN
654  rcmin = bignum
655  rcmax = zero
656  DO 20 j = 1, n
657  rcmin = min( rcmin, c( j ) )
658  rcmax = max( rcmax, c( j ) )
659  20 continue
660  IF( rcmin.LE.zero ) THEN
661  info = -12
662  ELSE IF( n.GT.0 ) THEN
663  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
664  ELSE
665  colcnd = one
666  END IF
667  END IF
668  IF( info.EQ.0 ) THEN
669  IF( ldb.LT.max( 1, n ) ) THEN
670  info = -14
671  ELSE IF( ldx.LT.max( 1, n ) ) THEN
672  info = -16
673  END IF
674  END IF
675  END IF
676 *
677  IF( info.NE.0 ) THEN
678  CALL xerbla( 'DGESVXX', -info )
679  return
680  END IF
681 *
682  IF( equil ) THEN
683 *
684 * Compute row and column scalings to equilibrate the matrix A.
685 *
686  CALL dgeequb( n, n, a, lda, r, c, rowcnd, colcnd, amax,
687  $ infequ )
688  IF( infequ.EQ.0 ) THEN
689 *
690 * Equilibrate the matrix.
691 *
692  CALL dlaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
693  $ equed )
694  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
695  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
696  END IF
697 *
698 * If the scaling factors are not applied, set them to 1.0.
699 *
700  IF ( .NOT.rowequ ) THEN
701  DO j = 1, n
702  r( j ) = 1.0d+0
703  END DO
704  END IF
705  IF ( .NOT.colequ ) THEN
706  DO j = 1, n
707  c( j ) = 1.0d+0
708  END DO
709  END IF
710  END IF
711 *
712 * Scale the right-hand side.
713 *
714  IF( notran ) THEN
715  IF( rowequ ) CALL dlascl2( n, nrhs, r, b, ldb )
716  ELSE
717  IF( colequ ) CALL dlascl2( n, nrhs, c, b, ldb )
718  END IF
719 *
720  IF( nofact .OR. equil ) THEN
721 *
722 * Compute the LU factorization of A.
723 *
724  CALL dlacpy( 'Full', n, n, a, lda, af, ldaf )
725  CALL dgetrf( n, n, af, ldaf, ipiv, info )
726 *
727 * Return if INFO is non-zero.
728 *
729  IF( info.GT.0 ) THEN
730 *
731 * Pivot in column INFO is exactly 0
732 * Compute the reciprocal pivot growth factor of the
733 * leading rank-deficient INFO columns of A.
734 *
735  rpvgrw = dla_gerpvgrw( n, info, a, lda, af, ldaf )
736  return
737  END IF
738  END IF
739 *
740 * Compute the reciprocal pivot growth factor RPVGRW.
741 *
742  rpvgrw = dla_gerpvgrw( n, n, a, lda, af, ldaf )
743 *
744 * Compute the solution matrix X.
745 *
746  CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
747  CALL dgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
748 *
749 * Use iterative refinement to improve the computed solution and
750 * compute error bounds and backward error estimates for it.
751 *
752  CALL dgerfsx( trans, equed, n, nrhs, a, lda, af, ldaf,
753  $ ipiv, r, c, b, ldb, x, ldx, rcond, berr,
754  $ n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params,
755  $ work, iwork, info )
756 *
757 * Scale solutions.
758 *
759  IF ( colequ .AND. notran ) THEN
760  CALL dlascl2( n, nrhs, c, x, ldx )
761  ELSE IF ( rowequ .AND. .NOT.notran ) THEN
762  CALL dlascl2( n, nrhs, r, x, ldx )
763  END IF
764 *
765  return
766 *
767 * End of DGESVXX
768 
769  END