LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cgerfsx.f
Go to the documentation of this file.
1*> \brief \b CGERFSX
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGERFSX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgerfsx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgerfsx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgerfsx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGERFSX( 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, RWORK, 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( * )
34* COMPLEX 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, * ), RWORK( * )
39* ..
40*
41*
42*> \par Purpose:
43* =============
44*>
45*> \verbatim
46*>
47*> CGERFSX 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)
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 COMPLEX 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 COMPLEX array, dimension (LDAF,N)
124*> The factors L and U from the factorization A = P*L*U
125*> as computed by CGETRF.
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 CGETRF; 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 COMPLEX 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 COMPLEX array, dimension (LDX,NRHS)
186*> On entry, the solution matrix X, as computed by CGETRS.
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 COMPLEX array, dimension (2*N)
365*> \endverbatim
366*>
367*> \param[out] RWORK
368*> \verbatim
369*> RWORK is REAL array, dimension (2*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 gerfsx
408*
409* =====================================================================
410 SUBROUTINE cgerfsx( 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, RWORK, 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( * )
427 COMPLEX 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, * ), rwork( * )
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, transfer
473* ..
474* .. External Functions ..
475 EXTERNAL lsame, ilatrans, ilaprec
477 REAL SLAMCH, CLANGE, CLA_GERCOND_X, CLA_GERCOND_C
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( 'CGERFSX', -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 = clange( norm, n, n, a, lda, rwork )
609 CALL cgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, 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 cla_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, rwork, work(n+1),
622 $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
623 $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
624 $ info )
625 ELSE
626 CALL cla_gerfsx_extended( prec_type, trans_type, n,
627 $ nrhs, a, lda, af, ldaf, ipiv, rowequ, r, b,
628 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
629 $ err_bnds_comp, work, rwork, work(n+1),
630 $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
631 $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
632 $ info )
633 END IF
634 END IF
635
636 err_lbnd = max( 10.0, sqrt( real( n ) ) ) * slamch( 'Epsilon' )
637 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
638*
639* Compute scaled normwise condition number cond(A*C).
640*
641 IF ( colequ .AND. notran ) THEN
642 rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
643 $ c, .true., info, work, rwork )
644 ELSE IF ( rowequ .AND. .NOT. notran ) THEN
645 rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
646 $ r, .true., info, work, rwork )
647 ELSE
648 rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
649 $ c, .false., info, work, rwork )
650 END IF
651 DO j = 1, nrhs
652*
653* Cap the error at 1.0.
654*
655 IF ( n_err_bnds .GE. la_linrx_err_i
656 $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0 )
657 $ err_bnds_norm( j, la_linrx_err_i ) = 1.0
658*
659* Threshold the error (see LAWN).
660*
661 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
662 err_bnds_norm( j, la_linrx_err_i ) = 1.0
663 err_bnds_norm( j, la_linrx_trust_i ) = 0.0
664 IF ( info .LE. n ) info = n + j
665 ELSE IF (err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd)
666 $ THEN
667 err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
668 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
669 END IF
670*
671* Save the condition number.
672*
673 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
674 err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
675 END IF
676 END DO
677 END IF
678
679 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
680*
681* Compute componentwise condition number cond(A*diag(Y(:,J))) for
682* each right-hand side using the current solution as an estimate of
683* the true solution. If the componentwise error estimate is too
684* large, then the solution is a lousy estimate of truth and the
685* estimated RCOND may be too optimistic. To avoid misleading users,
686* the inverse condition number is set to 0.0 when the estimated
687* cwise error is at least CWISE_WRONG.
688*
689 cwise_wrong = sqrt( slamch( 'Epsilon' ) )
690 DO j = 1, nrhs
691 IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
692 $ THEN
693 rcond_tmp = cla_gercond_x( trans, n, a, lda, af, ldaf,
694 $ ipiv, x(1,j), info, work, rwork )
695 ELSE
696 rcond_tmp = 0.0
697 END IF
698*
699* Cap the error at 1.0.
700*
701 IF ( n_err_bnds .GE. la_linrx_err_i
702 $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0 )
703 $ err_bnds_comp( j, la_linrx_err_i ) = 1.0
704*
705* Threshold the error (see LAWN).
706*
707 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
708 err_bnds_comp( j, la_linrx_err_i ) = 1.0
709 err_bnds_comp( j, la_linrx_trust_i ) = 0.0
710 IF ( params( la_linrx_cwise_i ) .EQ. 1.0
711 $ .AND. info.LT.n + j ) info = n + j
712 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
713 $ .LT. err_lbnd ) THEN
714 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
715 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
716 END IF
717*
718* Save the condition number.
719*
720 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
721 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
722 END IF
723
724 END DO
725 END IF
726*
727 RETURN
728*
729* End of CGERFSX
730*
731 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgecon(norm, n, a, lda, anorm, rcond, work, rwork, info)
CGECON
Definition cgecon.f:132
subroutine cgerfsx(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, rwork, info)
CGERFSX
Definition cgerfsx.f:414
integer function ilaprec(prec)
ILAPREC
Definition ilaprec.f:58
integer function ilatrans(trans)
ILATRANS
Definition ilatrans.f:58
real function cla_gercond_c(trans, n, a, lda, af, ldaf, ipiv, c, capply, info, work, rwork)
CLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices.
real function cla_gercond_x(trans, n, a, lda, af, ldaf, ipiv, x, info, work, rwork)
CLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices.
subroutine cla_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)
CLA_GERFSX_EXTENDED
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48