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