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