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