1*> \brief \b CPORFSX 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CPORFSX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cporfsx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cporfsx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cporfsx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CPORFSX( 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* REAL RCOND 31* .. 32* .. Array Arguments .. 33* COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 34* $ X( LDX, * ), WORK( * ) 35* REAL 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*> CPORFSX 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 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 array, dimension (LDAF,N) 124*> The triangular factor U or L from the Cholesky factorization 125*> A = U**T*U or A = L*L**T, as computed by SPOTRF. 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 REAL 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 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 array, dimension (LDX,NRHS) 165*> On entry, the solution matrix X, as computed by SGETRS. 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 REAL 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 REAL 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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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.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 array, dimension (2*N) 344*> \endverbatim 345*> 346*> \param[out] RWORK 347*> \verbatim 348*> RWORK is REAL 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*> \date April 2012 387* 388*> \ingroup complexPOcomputational 389* 390* ===================================================================== 391 SUBROUTINE CPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B, 392 $ LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, 393 $ ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, 394 $ WORK, RWORK, INFO ) 395* 396* -- LAPACK computational routine (version 3.7.0) -- 397* -- LAPACK is a software package provided by Univ. of Tennessee, -- 398* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 399* April 2012 400* 401* .. Scalar Arguments .. 402 CHARACTER UPLO, EQUED 403 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 404 $ N_ERR_BNDS 405 REAL RCOND 406* .. 407* .. Array Arguments .. 408 COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 409 $ X( LDX, * ), WORK( * ) 410 REAL RWORK( * ), S( * ), PARAMS(*), BERR( * ), 411 $ ERR_BNDS_NORM( NRHS, * ), 412 $ ERR_BNDS_COMP( NRHS, * ) 413* .. 414* 415* ================================================================== 416* 417* .. Parameters .. 418 REAL ZERO, ONE 419 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 420 REAL ITREF_DEFAULT, ITHRESH_DEFAULT, 421 $ COMPONENTWISE_DEFAULT 422 REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT 423 PARAMETER ( ITREF_DEFAULT = 1.0 ) 424 PARAMETER ( ITHRESH_DEFAULT = 10.0 ) 425 PARAMETER ( COMPONENTWISE_DEFAULT = 1.0 ) 426 PARAMETER ( RTHRESH_DEFAULT = 0.5 ) 427 PARAMETER ( DZTHRESH_DEFAULT = 0.25 ) 428 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 429 $ LA_LINRX_CWISE_I 430 PARAMETER ( LA_LINRX_ITREF_I = 1, 431 $ LA_LINRX_ITHRESH_I = 2 ) 432 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 433 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 434 $ LA_LINRX_RCOND_I 435 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 436 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 437* .. 438* .. Local Scalars .. 439 CHARACTER(1) NORM 440 LOGICAL RCEQU 441 INTEGER J, PREC_TYPE, REF_TYPE 442 INTEGER N_NORMS 443 REAL ANORM, RCOND_TMP 444 REAL ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG 445 LOGICAL IGNORE_CWISE 446 INTEGER ITHRESH 447 REAL RTHRESH, UNSTABLE_THRESH 448* .. 449* .. External Subroutines .. 450 EXTERNAL XERBLA, CPOCON, CLA_PORFSX_EXTENDED 451* .. 452* .. Intrinsic Functions .. 453 INTRINSIC MAX, SQRT, TRANSFER 454* .. 455* .. External Functions .. 456 EXTERNAL LSAME, ILAPREC 457 EXTERNAL SLAMCH, CLANHE, CLA_PORCOND_X, CLA_PORCOND_C 458 REAL SLAMCH, CLANHE, CLA_PORCOND_X, CLA_PORCOND_C 459 LOGICAL LSAME 460 INTEGER ILAPREC 461* .. 462* .. Executable Statements .. 463* 464* Check the input parameters. 465* 466 INFO = 0 467 REF_TYPE = INT( ITREF_DEFAULT ) 468 IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN 469 IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0 ) THEN 470 PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT 471 ELSE 472 REF_TYPE = PARAMS( LA_LINRX_ITREF_I ) 473 END IF 474 END IF 475* 476* Set default parameters. 477* 478 ILLRCOND_THRESH = REAL( N ) * SLAMCH( 'Epsilon' ) 479 ITHRESH = INT( ITHRESH_DEFAULT ) 480 RTHRESH = RTHRESH_DEFAULT 481 UNSTABLE_THRESH = DZTHRESH_DEFAULT 482 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0 483* 484 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN 485 IF ( PARAMS(LA_LINRX_ITHRESH_I ).LT.0.0 ) THEN 486 PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH 487 ELSE 488 ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) ) 489 END IF 490 END IF 491 IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN 492 IF ( PARAMS(LA_LINRX_CWISE_I ).LT.0.0 ) THEN 493 IF ( IGNORE_CWISE ) THEN 494 PARAMS( LA_LINRX_CWISE_I ) = 0.0 495 ELSE 496 PARAMS( LA_LINRX_CWISE_I ) = 1.0 497 END IF 498 ELSE 499 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0 500 END IF 501 END IF 502 IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN 503 N_NORMS = 0 504 ELSE IF ( IGNORE_CWISE ) THEN 505 N_NORMS = 1 506 ELSE 507 N_NORMS = 2 508 END IF 509* 510 RCEQU = LSAME( EQUED, 'Y' ) 511* 512* Test input parameters. 513* 514 IF (.NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 515 INFO = -1 516 ELSE IF( .NOT.RCEQU .AND. .NOT.LSAME( EQUED, 'N' ) ) THEN 517 INFO = -2 518 ELSE IF( N.LT.0 ) THEN 519 INFO = -3 520 ELSE IF( NRHS.LT.0 ) THEN 521 INFO = -4 522 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 523 INFO = -6 524 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 525 INFO = -8 526 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 527 INFO = -11 528 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 529 INFO = -13 530 END IF 531 IF( INFO.NE.0 ) THEN 532 CALL XERBLA( 'CPORFSX', -INFO ) 533 RETURN 534 END IF 535* 536* Quick return if possible. 537* 538 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 539 RCOND = 1.0 540 DO J = 1, NRHS 541 BERR( J ) = 0.0 542 IF ( N_ERR_BNDS .GE. 1 ) THEN 543 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0 544 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0 545 END IF 546 IF ( N_ERR_BNDS .GE. 2 ) THEN 547 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0 548 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0 549 END IF 550 IF ( N_ERR_BNDS .GE. 3 ) THEN 551 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0 552 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0 553 END IF 554 END DO 555 RETURN 556 END IF 557* 558* Default to failure. 559* 560 RCOND = 0.0 561 DO J = 1, NRHS 562 BERR( J ) = 1.0 563 IF ( N_ERR_BNDS .GE. 1 ) THEN 564 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0 565 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0 566 END IF 567 IF ( N_ERR_BNDS .GE. 2 ) THEN 568 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0 569 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0 570 END IF 571 IF ( N_ERR_BNDS .GE. 3 ) THEN 572 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0 573 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0 574 END IF 575 END DO 576* 577* Compute the norm of A and the reciprocal of the condition 578* number of A. 579* 580 NORM = 'I' 581 ANORM = CLANHE( NORM, UPLO, N, A, LDA, RWORK ) 582 CALL CPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, 583 $ INFO ) 584* 585* Perform refinement on each right-hand side 586* 587 IF ( REF_TYPE .NE. 0 ) THEN 588 589 PREC_TYPE = ILAPREC( 'D' ) 590 591 CALL CLA_PORFSX_EXTENDED( PREC_TYPE, UPLO, N, 592 $ NRHS, A, LDA, AF, LDAF, RCEQU, S, B, 593 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP, 594 $ WORK, RWORK, WORK(N+1), 595 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), RCOND, 596 $ ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 597 $ INFO ) 598 END IF 599 600 ERR_LBND = MAX( 10.0, SQRT( REAL( N ) ) ) * SLAMCH( 'Epsilon' ) 601 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN 602* 603* Compute scaled normwise condition number cond(A*C). 604* 605 IF ( RCEQU ) THEN 606 RCOND_TMP = CLA_PORCOND_C( UPLO, N, A, LDA, AF, LDAF, 607 $ S, .TRUE., INFO, WORK, RWORK ) 608 ELSE 609 RCOND_TMP = CLA_PORCOND_C( UPLO, N, A, LDA, AF, LDAF, 610 $ S, .FALSE., INFO, WORK, RWORK ) 611 END IF 612 DO J = 1, NRHS 613* 614* Cap the error at 1.0. 615* 616 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 617 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0 ) 618 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0 619* 620* Threshold the error (see LAWN). 621* 622 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 623 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0 624 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0 625 IF ( INFO .LE. N ) INFO = N + J 626 ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND ) 627 $ THEN 628 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND 629 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0 630 END IF 631* 632* Save the condition number. 633* 634 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 635 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP 636 END IF 637 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( SLAMCH( 'Epsilon' ) ) 652 DO J = 1, NRHS 653 IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG ) 654 $ THEN 655 RCOND_TMP = CLA_PORCOND_X( UPLO, N, A, LDA, AF, LDAF, 656 $ X(1,J), INFO, WORK, RWORK ) 657 ELSE 658 RCOND_TMP = 0.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.0 ) 665 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.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.0 671 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0 672 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.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.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 CPORFSX 692* 693 END 694