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*> \date April 2012 388* 389*> \ingroup doublePOcomputational 390* 391* ===================================================================== 392 SUBROUTINE DPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B, 393 $ LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, 394 $ ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, 395 $ WORK, IWORK, INFO ) 396* 397* -- LAPACK computational routine (version 3.7.0) -- 398* -- LAPACK is a software package provided by Univ. of Tennessee, -- 399* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 400* April 2012 401* 402* .. Scalar Arguments .. 403 CHARACTER UPLO, EQUED 404 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 405 $ N_ERR_BNDS 406 DOUBLE PRECISION RCOND 407* .. 408* .. Array Arguments .. 409 INTEGER IWORK( * ) 410 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 411 $ X( LDX, * ), WORK( * ) 412 DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ), 413 $ ERR_BNDS_NORM( NRHS, * ), 414 $ ERR_BNDS_COMP( NRHS, * ) 415* .. 416* 417* ================================================================== 418* 419* .. Parameters .. 420 DOUBLE PRECISION ZERO, ONE 421 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 422 DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT 423 DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT 424 DOUBLE PRECISION DZTHRESH_DEFAULT 425 PARAMETER ( ITREF_DEFAULT = 1.0D+0 ) 426 PARAMETER ( ITHRESH_DEFAULT = 10.0D+0 ) 427 PARAMETER ( COMPONENTWISE_DEFAULT = 1.0D+0 ) 428 PARAMETER ( RTHRESH_DEFAULT = 0.5D+0 ) 429 PARAMETER ( DZTHRESH_DEFAULT = 0.25D+0 ) 430 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 431 $ LA_LINRX_CWISE_I 432 PARAMETER ( LA_LINRX_ITREF_I = 1, 433 $ LA_LINRX_ITHRESH_I = 2 ) 434 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 435 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 436 $ LA_LINRX_RCOND_I 437 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 438 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 439* .. 440* .. Local Scalars .. 441 CHARACTER(1) NORM 442 LOGICAL RCEQU 443 INTEGER J, PREC_TYPE, REF_TYPE 444 INTEGER N_NORMS 445 DOUBLE PRECISION ANORM, RCOND_TMP 446 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG 447 LOGICAL IGNORE_CWISE 448 INTEGER ITHRESH 449 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH 450* .. 451* .. External Subroutines .. 452 EXTERNAL XERBLA, DPOCON, DLA_PORFSX_EXTENDED 453* .. 454* .. Intrinsic Functions .. 455 INTRINSIC MAX, SQRT 456* .. 457* .. External Functions .. 458 EXTERNAL LSAME, ILAPREC 459 EXTERNAL DLAMCH, DLANSY, DLA_PORCOND 460 DOUBLE PRECISION DLAMCH, DLANSY, DLA_PORCOND 461 LOGICAL LSAME 462 INTEGER ILAPREC 463* .. 464* .. Executable Statements .. 465* 466* Check the input parameters. 467* 468 INFO = 0 469 REF_TYPE = INT( ITREF_DEFAULT ) 470 IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN 471 IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN 472 PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT 473 ELSE 474 REF_TYPE = PARAMS( LA_LINRX_ITREF_I ) 475 END IF 476 END IF 477* 478* Set default parameters. 479* 480 ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' ) 481 ITHRESH = INT( ITHRESH_DEFAULT ) 482 RTHRESH = RTHRESH_DEFAULT 483 UNSTABLE_THRESH = DZTHRESH_DEFAULT 484 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0 485* 486 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN 487 IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN 488 PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH 489 ELSE 490 ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) ) 491 END IF 492 END IF 493 IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN 494 IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN 495 IF ( IGNORE_CWISE ) THEN 496 PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0 497 ELSE 498 PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0 499 END IF 500 ELSE 501 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0 502 END IF 503 END IF 504 IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN 505 N_NORMS = 0 506 ELSE IF ( IGNORE_CWISE ) THEN 507 N_NORMS = 1 508 ELSE 509 N_NORMS = 2 510 END IF 511* 512 RCEQU = LSAME( EQUED, 'Y' ) 513* 514* Test input parameters. 515* 516 IF (.NOT.LSAME(UPLO, 'U') .AND. .NOT.LSAME(UPLO, 'L')) THEN 517 INFO = -1 518 ELSE IF( .NOT.RCEQU .AND. .NOT.LSAME( EQUED, 'N' ) ) THEN 519 INFO = -2 520 ELSE IF( N.LT.0 ) THEN 521 INFO = -3 522 ELSE IF( NRHS.LT.0 ) THEN 523 INFO = -4 524 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 525 INFO = -6 526 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 527 INFO = -8 528 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 529 INFO = -11 530 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 531 INFO = -13 532 END IF 533 IF( INFO.NE.0 ) THEN 534 CALL XERBLA( 'DPORFSX', -INFO ) 535 RETURN 536 END IF 537* 538* Quick return if possible. 539* 540 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 541 RCOND = 1.0D+0 542 DO J = 1, NRHS 543 BERR( J ) = 0.0D+0 544 IF ( N_ERR_BNDS .GE. 1 ) THEN 545 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 546 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 547 END IF 548 IF ( N_ERR_BNDS .GE. 2 ) THEN 549 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0 550 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0 551 END IF 552 IF ( N_ERR_BNDS .GE. 3 ) THEN 553 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0 554 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0 555 END IF 556 END DO 557 RETURN 558 END IF 559* 560* Default to failure. 561* 562 RCOND = 0.0D+0 563 DO J = 1, NRHS 564 BERR( J ) = 1.0D+0 565 IF ( N_ERR_BNDS .GE. 1 ) THEN 566 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 567 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 568 END IF 569 IF ( N_ERR_BNDS .GE. 2 ) THEN 570 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 571 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 572 END IF 573 IF ( N_ERR_BNDS .GE. 3 ) THEN 574 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0 575 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0 576 END IF 577 END DO 578* 579* Compute the norm of A and the reciprocal of the condition 580* number of A. 581* 582 NORM = 'I' 583 ANORM = DLANSY( NORM, UPLO, N, A, LDA, WORK ) 584 CALL DPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, 585 $ IWORK, INFO ) 586* 587* Perform refinement on each right-hand side 588* 589 IF ( REF_TYPE .NE. 0 ) THEN 590 591 PREC_TYPE = ILAPREC( 'E' ) 592 593 CALL DLA_PORFSX_EXTENDED( PREC_TYPE, UPLO, N, 594 $ NRHS, A, LDA, AF, LDAF, RCEQU, S, B, 595 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP, 596 $ WORK( N+1 ), WORK( 1 ), WORK( 2*N+1 ), WORK( 1 ), RCOND, 597 $ ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 598 $ INFO ) 599 END IF 600 601 ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' ) 602 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN 603* 604* Compute scaled normwise condition number cond(A*C). 605* 606 IF ( RCEQU ) THEN 607 RCOND_TMP = DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 608 $ -1, S, INFO, WORK, IWORK ) 609 ELSE 610 RCOND_TMP = DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 611 $ 0, S, INFO, WORK, IWORK ) 612 END IF 613 DO J = 1, NRHS 614* 615* Cap the error at 1.0. 616* 617 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 618 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 619 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 620* 621* Threshold the error (see LAWN). 622* 623 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 624 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 625 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0 626 IF ( INFO .LE. N ) INFO = N + J 627 ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND ) 628 $ THEN 629 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND 630 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 631 END IF 632* 633* Save the condition number. 634* 635 IF (N_ERR_BNDS .GE. LA_LINRX_RCOND_I) THEN 636 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP 637 END IF 638 END DO 639 END IF 640 641 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2 ) THEN 642* 643* Compute componentwise condition number cond(A*diag(Y(:,J))) for 644* each right-hand side using the current solution as an estimate of 645* the true solution. If the componentwise error estimate is too 646* large, then the solution is a lousy estimate of truth and the 647* estimated RCOND may be too optimistic. To avoid misleading users, 648* the inverse condition number is set to 0.0 when the estimated 649* cwise error is at least CWISE_WRONG. 650* 651 CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) ) 652 DO J = 1, NRHS 653 IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG ) 654 $ THEN 655 RCOND_TMP = DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 1, 656 $ X( 1, J ), INFO, WORK, IWORK ) 657 ELSE 658 RCOND_TMP = 0.0D+0 659 END IF 660* 661* Cap the error at 1.0. 662* 663 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 664 $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 665 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 666* 667* Threshold the error (see LAWN). 668* 669 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 670 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 671 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0 672 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0 673 $ .AND. INFO.LT.N + J ) INFO = N + J 674 ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) 675 $ .LT. ERR_LBND ) THEN 676 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND 677 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 678 END IF 679* 680* Save the condition number. 681* 682 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 683 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP 684 END IF 685 686 END DO 687 END IF 688* 689 RETURN 690* 691* End of DPORFSX 692* 693 END 694