LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
sgerfsx.f
Go to the documentation of this file.
1*> \brief \b SGERFSX
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgerfsx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgerfsx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgerfsx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
22* R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
23* ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
24* WORK, IWORK, INFO )
25*
26* .. Scalar Arguments ..
27* CHARACTER TRANS, EQUED
28* INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
29* \$ N_ERR_BNDS
30* REAL RCOND
31* ..
32* .. Array Arguments ..
33* INTEGER IPIV( * ), IWORK( * )
34* REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
35* \$ X( LDX , * ), WORK( * )
36* REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
37* \$ ERR_BNDS_NORM( NRHS, * ),
38* \$ ERR_BNDS_COMP( NRHS, * )
39* ..
40*
41*
42*> \par Purpose:
43* =============
44*>
45*> \verbatim
46*>
47*> SGERFSX improves the computed solution to a system of linear
48*> equations and provides error bounds and backward error estimates
49*> for the solution. In addition to normwise error bound, the code
50*> provides maximum componentwise error bound if possible. See
51*> comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
52*> error bounds.
53*>
54*> The original system of linear equations may have been equilibrated
55*> before calling this routine, as described by arguments EQUED, R
56*> and C below. In this case, the solution and error bounds returned
57*> are for the original unequilibrated system.
58*> \endverbatim
59*
60* Arguments:
61* ==========
62*
63*> \verbatim
64*> Some optional parameters are bundled in the PARAMS array. These
65*> settings determine how refinement is performed, but often the
66*> defaults are acceptable. If the defaults are acceptable, users
67*> can pass NPARAMS = 0 which prevents the source code from accessing
68*> the PARAMS argument.
69*> \endverbatim
70*>
71*> \param[in] TRANS
72*> \verbatim
73*> TRANS is CHARACTER*1
74*> Specifies the form of the system of equations:
75*> = 'N': A * X = B (No transpose)
76*> = 'T': A**T * X = B (Transpose)
77*> = 'C': A**H * X = B (Conjugate transpose = Transpose)
78*> \endverbatim
79*>
80*> \param[in] EQUED
81*> \verbatim
82*> EQUED is CHARACTER*1
83*> Specifies the form of equilibration that was done to A
84*> before calling this routine. This is needed to compute
85*> the solution and error bounds correctly.
86*> = 'N': No equilibration
87*> = 'R': Row equilibration, i.e., A has been premultiplied by
88*> diag(R).
89*> = 'C': Column equilibration, i.e., A has been postmultiplied
90*> by diag(C).
91*> = 'B': Both row and column equilibration, i.e., A has been
92*> replaced by diag(R) * A * diag(C).
93*> The right hand side B has been changed accordingly.
94*> \endverbatim
95*>
96*> \param[in] N
97*> \verbatim
98*> N is INTEGER
99*> The order of the matrix A. N >= 0.
100*> \endverbatim
101*>
102*> \param[in] NRHS
103*> \verbatim
104*> NRHS is INTEGER
105*> The number of right hand sides, i.e., the number of columns
106*> of the matrices B and X. NRHS >= 0.
107*> \endverbatim
108*>
109*> \param[in] A
110*> \verbatim
111*> A is REAL array, dimension (LDA,N)
112*> The original N-by-N matrix A.
113*> \endverbatim
114*>
115*> \param[in] LDA
116*> \verbatim
117*> LDA is INTEGER
118*> The leading dimension of the array A. LDA >= max(1,N).
119*> \endverbatim
120*>
121*> \param[in] AF
122*> \verbatim
123*> AF is REAL array, dimension (LDAF,N)
124*> The factors L and U from the factorization A = P*L*U
125*> as computed by SGETRF.
126*> \endverbatim
127*>
128*> \param[in] LDAF
129*> \verbatim
130*> LDAF is INTEGER
131*> The leading dimension of the array AF. LDAF >= max(1,N).
132*> \endverbatim
133*>
134*> \param[in] IPIV
135*> \verbatim
136*> IPIV is INTEGER array, dimension (N)
137*> The pivot indices from SGETRF; for 1<=i<=N, row i of the
138*> matrix was interchanged with row IPIV(i).
139*> \endverbatim
140*>
141*> \param[in] R
142*> \verbatim
143*> R is REAL array, dimension (N)
144*> The row scale factors for A. If EQUED = 'R' or 'B', A is
145*> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
146*> is not accessed.
147*> If R is accessed, each element of R should be a power of the radix
148*> to ensure a reliable solution and error estimates. Scaling by
149*> powers of the radix does not cause rounding errors unless the
150*> result underflows or overflows. Rounding errors during scaling
151*> lead to refining with a matrix that is not equivalent to the
152*> input matrix, producing error estimates that may not be
153*> reliable.
154*> \endverbatim
155*>
156*> \param[in] C
157*> \verbatim
158*> C is REAL array, dimension (N)
159*> The column scale factors for A. If EQUED = 'C' or 'B', A is
160*> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
161*> is not accessed.
162*> If C is accessed, each element of C should be a power of the radix
163*> to ensure a reliable solution and error estimates. Scaling by
164*> powers of the radix does not cause rounding errors unless the
165*> result underflows or overflows. Rounding errors during scaling
166*> lead to refining with a matrix that is not equivalent to the
167*> input matrix, producing error estimates that may not be
168*> reliable.
169*> \endverbatim
170*>
171*> \param[in] B
172*> \verbatim
173*> B is REAL array, dimension (LDB,NRHS)
174*> The right hand side matrix B.
175*> \endverbatim
176*>
177*> \param[in] LDB
178*> \verbatim
179*> LDB is INTEGER
180*> The leading dimension of the array B. LDB >= max(1,N).
181*> \endverbatim
182*>
183*> \param[in,out] X
184*> \verbatim
185*> X is REAL array, dimension (LDX,NRHS)
186*> On entry, the solution matrix X, as computed by SGETRS.
187*> On exit, the improved solution matrix X.
188*> \endverbatim
189*>
190*> \param[in] LDX
191*> \verbatim
192*> LDX is INTEGER
193*> The leading dimension of the array X. LDX >= max(1,N).
194*> \endverbatim
195*>
196*> \param[out] RCOND
197*> \verbatim
198*> RCOND is REAL
199*> Reciprocal scaled condition number. This is an estimate of the
200*> reciprocal Skeel condition number of the matrix A after
201*> equilibration (if done). If this is less than the machine
202*> precision (in particular, if it is zero), the matrix is singular
203*> to working precision. Note that the error may still be small even
204*> if this number is very small and the matrix appears ill-
205*> conditioned.
206*> \endverbatim
207*>
208*> \param[out] BERR
209*> \verbatim
210*> BERR is REAL array, dimension (NRHS)
211*> Componentwise relative backward error. This is the
212*> componentwise relative backward error of each solution vector X(j)
213*> (i.e., the smallest relative change in any element of A or B that
214*> makes X(j) an exact solution).
215*> \endverbatim
216*>
217*> \param[in] N_ERR_BNDS
218*> \verbatim
219*> N_ERR_BNDS is INTEGER
220*> Number of error bounds to return for each right hand side
221*> and each type (normwise or componentwise). See ERR_BNDS_NORM and
222*> ERR_BNDS_COMP below.
223*> \endverbatim
224*>
225*> \param[out] ERR_BNDS_NORM
226*> \verbatim
227*> ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
228*> For each right-hand side, this array contains information about
229*> various error bounds and condition numbers corresponding to the
230*> normwise relative error, which is defined as follows:
231*>
232*> Normwise relative error in the ith solution vector:
233*> max_j (abs(XTRUE(j,i) - X(j,i)))
234*> ------------------------------
235*> max_j abs(X(j,i))
236*>
237*> The array is indexed by the type of error information as described
238*> below. There currently are up to three pieces of information
239*> returned.
240*>
241*> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
242*> right-hand side.
243*>
244*> The second index in ERR_BNDS_NORM(:,err) contains the following
245*> three fields:
246*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
247*> reciprocal condition number is less than the threshold
248*> sqrt(n) * slamch('Epsilon').
249*>
250*> err = 2 "Guaranteed" error bound: The estimated forward error,
251*> almost certainly within a factor of 10 of the true error
252*> so long as the next entry is greater than the threshold
253*> sqrt(n) * slamch('Epsilon'). This error bound should only
254*> be trusted if the previous boolean is true.
255*>
256*> err = 3 Reciprocal condition number: Estimated normwise
257*> reciprocal condition number. Compared with the threshold
258*> sqrt(n) * slamch('Epsilon') to determine if the error
259*> estimate is "guaranteed". These reciprocal condition
260*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
261*> appropriately scaled matrix Z.
262*> Let Z = S*A, where S scales each row by a power of the
263*> radix so all absolute row sums of Z are approximately 1.
264*>
265*> See Lapack Working Note 165 for further details and extra
266*> cautions.
267*> \endverbatim
268*>
269*> \param[out] ERR_BNDS_COMP
270*> \verbatim
271*> ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
272*> For each right-hand side, this array contains information about
273*> various error bounds and condition numbers corresponding to the
274*> componentwise relative error, which is defined as follows:
275*>
276*> Componentwise relative error in the ith solution vector:
277*> abs(XTRUE(j,i) - X(j,i))
278*> max_j ----------------------
279*> abs(X(j,i))
280*>
281*> The array is indexed by the right-hand side i (on which the
282*> componentwise relative error depends), and the type of error
283*> information as described below. There currently are up to three
284*> pieces of information returned for each right-hand side. If
285*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
286*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
287*> the first (:,N_ERR_BNDS) entries are returned.
288*>
289*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
290*> right-hand side.
291*>
292*> The second index in ERR_BNDS_COMP(:,err) contains the following
293*> three fields:
294*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
295*> reciprocal condition number is less than the threshold
296*> sqrt(n) * slamch('Epsilon').
297*>
298*> err = 2 "Guaranteed" error bound: The estimated forward error,
299*> almost certainly within a factor of 10 of the true error
300*> so long as the next entry is greater than the threshold
301*> sqrt(n) * slamch('Epsilon'). This error bound should only
302*> be trusted if the previous boolean is true.
303*>
304*> err = 3 Reciprocal condition number: Estimated componentwise
305*> reciprocal condition number. Compared with the threshold
306*> sqrt(n) * slamch('Epsilon') to determine if the error
307*> estimate is "guaranteed". These reciprocal condition
308*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
309*> appropriately scaled matrix Z.
310*> Let Z = S*(A*diag(x)), where x is the solution for the
311*> current right-hand side and S scales each row of
312*> A*diag(x) by a power of the radix so all absolute row
313*> sums of Z are approximately 1.
314*>
315*> See Lapack Working Note 165 for further details and extra
316*> cautions.
317*> \endverbatim
318*>
319*> \param[in] NPARAMS
320*> \verbatim
321*> NPARAMS is INTEGER
322*> Specifies the number of parameters set in PARAMS. If <= 0, the
323*> PARAMS array is never referenced and default values are used.
324*> \endverbatim
325*>
326*> \param[in,out] PARAMS
327*> \verbatim
328*> PARAMS is REAL array, dimension NPARAMS
329*> Specifies algorithm parameters. If an entry is < 0.0, then
330*> that entry will be filled with default value used for that
331*> parameter. Only positions up to NPARAMS are accessed; defaults
332*> are used for higher-numbered parameters.
333*>
334*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
335*> refinement or not.
336*> Default: 1.0
337*> = 0.0: No refinement is performed, and no error bounds are
338*> computed.
339*> = 1.0: Use the double-precision refinement algorithm,
340*> possibly with doubled-single computations if the
341*> compilation environment does not support DOUBLE
342*> PRECISION.
343*> (other values are reserved for future use)
344*>
345*> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
346*> computations allowed for refinement.
347*> Default: 10
348*> Aggressive: Set to 100 to permit convergence using approximate
349*> factorizations or factorizations other than LU. If
350*> the factorization uses a technique other than
351*> Gaussian elimination, the guarantees in
352*> err_bnds_norm and err_bnds_comp may no longer be
353*> trustworthy.
354*>
355*> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
356*> will attempt to find a solution with small componentwise
357*> relative error in the double-precision algorithm. Positive
358*> is true, 0.0 is false.
359*> Default: 1.0 (attempt componentwise convergence)
360*> \endverbatim
361*>
362*> \param[out] WORK
363*> \verbatim
364*> WORK is REAL array, dimension (4*N)
365*> \endverbatim
366*>
367*> \param[out] IWORK
368*> \verbatim
369*> IWORK is INTEGER array, dimension (N)
370*> \endverbatim
371*>
372*> \param[out] INFO
373*> \verbatim
374*> INFO is INTEGER
375*> = 0: Successful exit. The solution to every right-hand side is
376*> guaranteed.
377*> < 0: If INFO = -i, the i-th argument had an illegal value
378*> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
379*> has been completed, but the factor U is exactly singular, so
380*> the solution and error bounds could not be computed. RCOND = 0
381*> is returned.
382*> = N+J: The solution corresponding to the Jth right-hand side is
383*> not guaranteed. The solutions corresponding to other right-
384*> hand sides K with K > J may not be guaranteed as well, but
385*> only the first such right-hand side is reported. If a small
386*> componentwise error is not requested (PARAMS(3) = 0.0) then
387*> the Jth right-hand side is the first with a normwise error
388*> bound that is not guaranteed (the smallest J such
389*> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
390*> the Jth right-hand side is the first with either a normwise or
391*> componentwise error bound that is not guaranteed (the smallest
392*> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
393*> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
394*> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
395*> about all of the right-hand sides check ERR_BNDS_NORM or
396*> ERR_BNDS_COMP.
397*> \endverbatim
398*
399* Authors:
400* ========
401*
402*> \author Univ. of Tennessee
403*> \author Univ. of California Berkeley
404*> \author Univ. of Colorado Denver
405*> \author NAG Ltd.
406*
407*> \ingroup realGEcomputational
408*
409* =====================================================================
410 SUBROUTINE sgerfsx( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
411 \$ R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
412 \$ ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
413 \$ WORK, IWORK, INFO )
414*
415* -- LAPACK computational routine --
416* -- LAPACK is a software package provided by Univ. of Tennessee, --
417* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
418*
419* .. Scalar Arguments ..
420 CHARACTER TRANS, EQUED
421 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
422 \$ N_ERR_BNDS
423 REAL RCOND
424* ..
425* .. Array Arguments ..
426 INTEGER IPIV( * ), IWORK( * )
427 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
428 \$ X( LDX , * ), WORK( * )
429 REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
430 \$ err_bnds_norm( nrhs, * ),
431 \$ err_bnds_comp( nrhs, * )
432* ..
433*
434* ==================================================================
435*
436* .. Parameters ..
437 REAL ZERO, ONE
438 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
439 REAL ITREF_DEFAULT, ITHRESH_DEFAULT,
440 \$ componentwise_default
441 REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT
442 parameter( itref_default = 1.0 )
443 parameter( ithresh_default = 10.0 )
444 parameter( componentwise_default = 1.0 )
445 parameter( rthresh_default = 0.5 )
446 parameter( dzthresh_default = 0.25 )
447 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
448 \$ la_linrx_cwise_i
449 parameter( la_linrx_itref_i = 1,
450 \$ la_linrx_ithresh_i = 2 )
451 parameter( la_linrx_cwise_i = 3 )
452 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
453 \$ LA_LINRX_RCOND_I
454 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
455 parameter( la_linrx_rcond_i = 3 )
456* ..
457* .. Local Scalars ..
458 CHARACTER(1) NORM
459 LOGICAL ROWEQU, COLEQU, NOTRAN
460 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
461 INTEGER N_NORMS
462 REAL ANORM, RCOND_TMP
463 REAL ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
464 LOGICAL IGNORE_CWISE
465 INTEGER ITHRESH
466 REAL RTHRESH, UNSTABLE_THRESH
467* ..
468* .. External Subroutines ..
470* ..
471* .. Intrinsic Functions ..
472 INTRINSIC max, sqrt
473* ..
474* .. External Functions ..
475 EXTERNAL lsame, ilatrans, ilaprec
476 EXTERNAL slamch, slange, sla_gercond
477 REAL SLAMCH, SLANGE, SLA_GERCOND
478 LOGICAL LSAME
479 INTEGER ILATRANS, ILAPREC
480* ..
481* .. Executable Statements ..
482*
483* Check the input parameters.
484*
485 info = 0
486 trans_type = ilatrans( trans )
487 ref_type = int( itref_default )
488 IF ( nparams .GE. la_linrx_itref_i ) THEN
489 IF ( params( la_linrx_itref_i ) .LT. 0.0 ) THEN
490 params( la_linrx_itref_i ) = itref_default
491 ELSE
492 ref_type = params( la_linrx_itref_i )
493 END IF
494 END IF
495*
496* Set default parameters.
497*
498 illrcond_thresh = real( n ) * slamch( 'Epsilon' )
499 ithresh = int( ithresh_default )
500 rthresh = rthresh_default
501 unstable_thresh = dzthresh_default
502 ignore_cwise = componentwise_default .EQ. 0.0
503*
504 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
505 IF ( params( la_linrx_ithresh_i ).LT.0.0 ) THEN
506 params( la_linrx_ithresh_i ) = ithresh
507 ELSE
508 ithresh = int( params( la_linrx_ithresh_i ) )
509 END IF
510 END IF
511 IF ( nparams.GE.la_linrx_cwise_i ) THEN
512 IF ( params( la_linrx_cwise_i ).LT.0.0 ) THEN
513 IF ( ignore_cwise ) THEN
514 params( la_linrx_cwise_i ) = 0.0
515 ELSE
516 params( la_linrx_cwise_i ) = 1.0
517 END IF
518 ELSE
519 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0
520 END IF
521 END IF
522 IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
523 n_norms = 0
524 ELSE IF ( ignore_cwise ) THEN
525 n_norms = 1
526 ELSE
527 n_norms = 2
528 END IF
529*
530 notran = lsame( trans, 'N' )
531 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
532 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
533*
534* Test input parameters.
535*
536 IF( trans_type.EQ.-1 ) THEN
537 info = -1
538 ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
539 \$ .NOT.lsame( equed, 'N' ) ) THEN
540 info = -2
541 ELSE IF( n.LT.0 ) THEN
542 info = -3
543 ELSE IF( nrhs.LT.0 ) THEN
544 info = -4
545 ELSE IF( lda.LT.max( 1, n ) ) THEN
546 info = -6
547 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
548 info = -8
549 ELSE IF( ldb.LT.max( 1, n ) ) THEN
550 info = -13
551 ELSE IF( ldx.LT.max( 1, n ) ) THEN
552 info = -15
553 END IF
554 IF( info.NE.0 ) THEN
555 CALL xerbla( 'SGERFSX', -info )
556 RETURN
557 END IF
558*
559* Quick return if possible.
560*
561 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
562 rcond = 1.0
563 DO j = 1, nrhs
564 berr( j ) = 0.0
565 IF ( n_err_bnds .GE. 1 ) THEN
566 err_bnds_norm( j, la_linrx_trust_i) = 1.0
567 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
568 END IF
569 IF ( n_err_bnds .GE. 2 ) THEN
570 err_bnds_norm( j, la_linrx_err_i) = 0.0
571 err_bnds_comp( j, la_linrx_err_i ) = 0.0
572 END IF
573 IF ( n_err_bnds .GE. 3 ) THEN
574 err_bnds_norm( j, la_linrx_rcond_i) = 1.0
575 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
576 END IF
577 END DO
578 RETURN
579 END IF
580*
581* Default to failure.
582*
583 rcond = 0.0
584 DO j = 1, nrhs
585 berr( j ) = 1.0
586 IF ( n_err_bnds .GE. 1 ) THEN
587 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
588 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
589 END IF
590 IF ( n_err_bnds .GE. 2 ) THEN
591 err_bnds_norm( j, la_linrx_err_i ) = 1.0
592 err_bnds_comp( j, la_linrx_err_i ) = 1.0
593 END IF
594 IF ( n_err_bnds .GE. 3 ) THEN
595 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
596 err_bnds_comp( j, la_linrx_rcond_i ) = 0.0
597 END IF
598 END DO
599*
600* Compute the norm of A and the reciprocal of the condition
601* number of A.
602*
603 IF( notran ) THEN
604 norm = 'I'
605 ELSE
606 norm = '1'
607 END IF
608 anorm = slange( norm, n, n, a, lda, work )
609 CALL sgecon( norm, n, af, ldaf, anorm, rcond, work, iwork, info )
610*
611* Perform refinement on each right-hand side
612*
613 IF ( ref_type .NE. 0 ) THEN
614
615 prec_type = ilaprec( 'D' )
616
617 IF ( notran ) THEN
618 CALL sla_gerfsx_extended( prec_type, trans_type, n,
619 \$ nrhs, a, lda, af, ldaf, ipiv, colequ, c, b,
620 \$ ldb, x, ldx, berr, n_norms, err_bnds_norm,
621 \$ err_bnds_comp, work(n+1), work(1), work(2*n+1),
622 \$ work(1), rcond, ithresh, rthresh, unstable_thresh,
623 \$ ignore_cwise, info )
624 ELSE
625 CALL sla_gerfsx_extended( prec_type, trans_type, n,
626 \$ nrhs, a, lda, af, ldaf, ipiv, rowequ, r, b,
627 \$ ldb, x, ldx, berr, n_norms, err_bnds_norm,
628 \$ err_bnds_comp, work(n+1), work(1), work(2*n+1),
629 \$ work(1), rcond, ithresh, rthresh, unstable_thresh,
630 \$ ignore_cwise, info )
631 END IF
632 END IF
633
634 err_lbnd = max( 10.0, sqrt( real( n ) ) ) * slamch( 'Epsilon' )
635 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
636*
637* Compute scaled normwise condition number cond(A*C).
638*
639 IF ( colequ .AND. notran ) THEN
640 rcond_tmp = sla_gercond( trans, n, a, lda, af, ldaf, ipiv,
641 \$ -1, c, info, work, iwork )
642 ELSE IF ( rowequ .AND. .NOT. notran ) THEN
643 rcond_tmp = sla_gercond( trans, n, a, lda, af, ldaf, ipiv,
644 \$ -1, r, info, work, iwork )
645 ELSE
646 rcond_tmp = sla_gercond( trans, n, a, lda, af, ldaf, ipiv,
647 \$ 0, r, info, work, iwork )
648 END IF
649 DO j = 1, nrhs
650*
651* Cap the error at 1.0.
652*
653 IF ( n_err_bnds .GE. la_linrx_err_i
654 \$ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0 )
655 \$ err_bnds_norm( j, la_linrx_err_i ) = 1.0
656*
657* Threshold the error (see LAWN).
658*
659 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
660 err_bnds_norm( j, la_linrx_err_i ) = 1.0
661 err_bnds_norm( j, la_linrx_trust_i ) = 0.0
662 IF ( info .LE. n ) info = n + j
663 ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
664 \$ THEN
665 err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
666 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
667 END IF
668*
669* Save the condition number.
670*
671 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
672 err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
673 END IF
674 END DO
675 END IF
676
677 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
678*
679* Compute componentwise condition number cond(A*diag(Y(:,J))) for
680* each right-hand side using the current solution as an estimate of
681* the true solution. If the componentwise error estimate is too
682* large, then the solution is a lousy estimate of truth and the
683* estimated RCOND may be too optimistic. To avoid misleading users,
684* the inverse condition number is set to 0.0 when the estimated
685* cwise error is at least CWISE_WRONG.
686*
687 cwise_wrong = sqrt( slamch( 'Epsilon' ) )
688 DO j = 1, nrhs
689 IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
690 \$ THEN
691 rcond_tmp = sla_gercond( trans, n, a, lda, af, ldaf,
692 \$ ipiv, 1, x(1,j), info, work, iwork )
693 ELSE
694 rcond_tmp = 0.0
695 END IF
696*
697* Cap the error at 1.0.
698*
699 IF ( n_err_bnds .GE. la_linrx_err_i
700 \$ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0 )
701 \$ err_bnds_comp( j, la_linrx_err_i ) = 1.0
702*
703* Threshold the error (see LAWN).
704*
705 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
706 err_bnds_comp( j, la_linrx_err_i ) = 1.0
707 err_bnds_comp( j, la_linrx_trust_i ) = 0.0
708 IF ( params( la_linrx_cwise_i ) .EQ. 1.0
709 \$ .AND. info.LT.n + j ) info = n + j
710 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
711 \$ .LT. err_lbnd ) THEN
712 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
713 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
714 END IF
715*
716* Save the condition number.
717*
718 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
719 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
720 END IF
721 END DO
722 END IF
723*
724 RETURN
725*
726* End of SGERFSX
727*
728 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:58
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:58
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
real function sla_gercond(TRANS, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
SLA_GERCOND estimates the Skeel condition number for a general matrix.
Definition: sla_gercond.f:150
subroutine sla_gerfsx_extended(PREC_TYPE, TRANS_TYPE, N, NRHS, A, LDA, AF, LDAF, IPIV, COLEQU, C, B, LDB, Y, LDY, BERR_OUT, N_NORMS, ERRS_N, ERRS_C, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)
SLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general matric...
subroutine sgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SGECON
Definition: sgecon.f:124
subroutine sgerfsx(TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)
SGERFSX
Definition: sgerfsx.f:414
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68