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