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