1*> \brief \b CGBRFSX 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CGBRFSX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbrfsx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbrfsx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbrfsx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CGBRFSX( 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* REAL RCOND 32* .. 33* .. Array Arguments .. 34* INTEGER IPIV( * ) 35* COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 36* $ X( LDX , * ),WORK( * ) 37* REAL 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*> CGBRFSX 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 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 array, dimension (LDAFB,N) 140*> Details of the LU factorization of the band matrix A, as 141*> computed by CGBTRF. 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 CGETRF; 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 REAL 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 REAL 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 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 array, dimension (LDX,NRHS) 211*> On entry, the solution matrix X, as computed by CGETRS. 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 REAL 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 REAL 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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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.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 array, dimension (2*N) 390*> \endverbatim 391*> 392*> \param[out] RWORK 393*> \verbatim 394*> RWORK is REAL 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*> \ingroup complexGBcomputational 433* 434* ===================================================================== 435 SUBROUTINE CGBRFSX( TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB, 436 $ LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND, 437 $ BERR, N_ERR_BNDS, ERR_BNDS_NORM, 438 $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, 439 $ INFO ) 440* 441* -- LAPACK computational routine -- 442* -- LAPACK is a software package provided by Univ. of Tennessee, -- 443* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 444* 445* .. Scalar Arguments .. 446 CHARACTER TRANS, EQUED 447 INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS, 448 $ NPARAMS, N_ERR_BNDS 449 REAL RCOND 450* .. 451* .. Array Arguments .. 452 INTEGER IPIV( * ) 453 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 454 $ X( LDX , * ),WORK( * ) 455 REAL R( * ), C( * ), PARAMS( * ), BERR( * ), 456 $ ERR_BNDS_NORM( NRHS, * ), 457 $ ERR_BNDS_COMP( NRHS, * ), RWORK( * ) 458* .. 459* 460* ================================================================== 461* 462* .. Parameters .. 463 REAL ZERO, ONE 464 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 465 REAL ITREF_DEFAULT, ITHRESH_DEFAULT, 466 $ COMPONENTWISE_DEFAULT 467 REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT 468 PARAMETER ( ITREF_DEFAULT = 1.0 ) 469 PARAMETER ( ITHRESH_DEFAULT = 10.0 ) 470 PARAMETER ( COMPONENTWISE_DEFAULT = 1.0 ) 471 PARAMETER ( RTHRESH_DEFAULT = 0.5 ) 472 PARAMETER ( DZTHRESH_DEFAULT = 0.25 ) 473 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 474 $ LA_LINRX_CWISE_I 475 PARAMETER ( LA_LINRX_ITREF_I = 1, 476 $ LA_LINRX_ITHRESH_I = 2 ) 477 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 478 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 479 $ LA_LINRX_RCOND_I 480 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 481 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 482* .. 483* .. Local Scalars .. 484 CHARACTER(1) NORM 485 LOGICAL ROWEQU, COLEQU, NOTRAN, IGNORE_CWISE 486 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE, N_NORMS, 487 $ ITHRESH 488 REAL ANORM, RCOND_TMP, ILLRCOND_THRESH, ERR_LBND, 489 $ CWISE_WRONG, RTHRESH, UNSTABLE_THRESH 490* .. 491* .. External Subroutines .. 492 EXTERNAL XERBLA, CGBCON, CLA_GBRFSX_EXTENDED 493* .. 494* .. Intrinsic Functions .. 495 INTRINSIC MAX, SQRT, TRANSFER 496* .. 497* .. External Functions .. 498 EXTERNAL LSAME, ILATRANS, ILAPREC 499 EXTERNAL SLAMCH, CLANGB, CLA_GBRCOND_X, CLA_GBRCOND_C 500 REAL SLAMCH, CLANGB, CLA_GBRCOND_X, CLA_GBRCOND_C 501 LOGICAL LSAME 502 INTEGER ILATRANS, ILAPREC 503* .. 504* .. Executable Statements .. 505* 506* Check the input parameters. 507* 508 INFO = 0 509 TRANS_TYPE = ILATRANS( TRANS ) 510 REF_TYPE = INT( ITREF_DEFAULT ) 511 IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN 512 IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0 ) THEN 513 PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT 514 ELSE 515 REF_TYPE = PARAMS( LA_LINRX_ITREF_I ) 516 END IF 517 END IF 518* 519* Set default parameters. 520* 521 ILLRCOND_THRESH = REAL( N ) * SLAMCH( 'Epsilon' ) 522 ITHRESH = INT( ITHRESH_DEFAULT ) 523 RTHRESH = RTHRESH_DEFAULT 524 UNSTABLE_THRESH = DZTHRESH_DEFAULT 525 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0 526* 527 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN 528 IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0 ) THEN 529 PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH 530 ELSE 531 ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) ) 532 END IF 533 END IF 534 IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN 535 IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0 ) THEN 536 IF ( IGNORE_CWISE ) THEN 537 PARAMS( LA_LINRX_CWISE_I ) = 0.0 538 ELSE 539 PARAMS( LA_LINRX_CWISE_I ) = 1.0 540 END IF 541 ELSE 542 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0 543 END IF 544 END IF 545 IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN 546 N_NORMS = 0 547 ELSE IF ( IGNORE_CWISE ) THEN 548 N_NORMS = 1 549 ELSE 550 N_NORMS = 2 551 END IF 552* 553 NOTRAN = LSAME( TRANS, 'N' ) 554 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' ) 555 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' ) 556* 557* Test input parameters. 558* 559 IF( TRANS_TYPE.EQ.-1 ) THEN 560 INFO = -1 561 ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND. 562 $ .NOT.LSAME( EQUED, 'N' ) ) THEN 563 INFO = -2 564 ELSE IF( N.LT.0 ) THEN 565 INFO = -3 566 ELSE IF( KL.LT.0 ) THEN 567 INFO = -4 568 ELSE IF( KU.LT.0 ) THEN 569 INFO = -5 570 ELSE IF( NRHS.LT.0 ) THEN 571 INFO = -6 572 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 573 INFO = -8 574 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 575 INFO = -10 576 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 577 INFO = -13 578 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 579 INFO = -15 580 END IF 581 IF( INFO.NE.0 ) THEN 582 CALL XERBLA( 'CGBRFSX', -INFO ) 583 RETURN 584 END IF 585* 586* Quick return if possible. 587* 588 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 589 RCOND = 1.0 590 DO J = 1, NRHS 591 BERR( J ) = 0.0 592 IF ( N_ERR_BNDS .GE. 1 ) THEN 593 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0 594 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0 595 END IF 596 IF ( N_ERR_BNDS .GE. 2 ) THEN 597 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0 598 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0 599 END IF 600 IF ( N_ERR_BNDS .GE. 3 ) THEN 601 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0 602 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0 603 END IF 604 END DO 605 RETURN 606 END IF 607* 608* Default to failure. 609* 610 RCOND = 0.0 611 DO J = 1, NRHS 612 BERR( J ) = 1.0 613 IF ( N_ERR_BNDS .GE. 1 ) THEN 614 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0 615 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0 616 END IF 617 IF ( N_ERR_BNDS .GE. 2 ) THEN 618 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0 619 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0 620 END IF 621 IF ( N_ERR_BNDS .GE. 3 ) THEN 622 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0 623 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0 624 END IF 625 END DO 626* 627* Compute the norm of A and the reciprocal of the condition 628* number of A. 629* 630 IF( NOTRAN ) THEN 631 NORM = 'I' 632 ELSE 633 NORM = '1' 634 END IF 635 ANORM = CLANGB( NORM, N, KL, KU, AB, LDAB, RWORK ) 636 CALL CGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND, 637 $ WORK, RWORK, INFO ) 638* 639* Perform refinement on each right-hand side 640* 641 IF ( REF_TYPE .NE. 0 .AND. INFO .EQ. 0 ) THEN 642 643 PREC_TYPE = ILAPREC( 'D' ) 644 645 IF ( NOTRAN ) THEN 646 CALL CLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 647 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B, 648 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, 649 $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1), 650 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), 651 $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 652 $ INFO ) 653 ELSE 654 CALL CLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 655 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, ROWEQU, R, B, 656 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, 657 $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1), 658 $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), 659 $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE, 660 $ INFO ) 661 END IF 662 END IF 663 664 ERR_LBND = MAX( 10.0, SQRT( REAL( N ) ) ) * SLAMCH( 'Epsilon' ) 665 IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1) THEN 666* 667* Compute scaled normwise condition number cond(A*C). 668* 669 IF ( COLEQU .AND. NOTRAN ) THEN 670 RCOND_TMP = CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB, 671 $ LDAFB, IPIV, C, .TRUE., INFO, WORK, RWORK ) 672 ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN 673 RCOND_TMP = CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB, 674 $ LDAFB, IPIV, R, .TRUE., INFO, WORK, RWORK ) 675 ELSE 676 RCOND_TMP = CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB, 677 $ LDAFB, IPIV, C, .FALSE., INFO, WORK, RWORK ) 678 END IF 679 DO J = 1, NRHS 680* 681* Cap the error at 1.0. 682* 683 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 684 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0) 685 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0 686* 687* Threshold the error (see LAWN). 688* 689 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 690 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0 691 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0 692 IF ( INFO .LE. N ) INFO = N + J 693 ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND ) 694 $ THEN 695 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND 696 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0 697 END IF 698* 699* Save the condition number. 700* 701 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 702 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP 703 END IF 704 705 END DO 706 END IF 707 708 IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN 709* 710* Compute componentwise condition number cond(A*diag(Y(:,J))) for 711* each right-hand side using the current solution as an estimate of 712* the true solution. If the componentwise error estimate is too 713* large, then the solution is a lousy estimate of truth and the 714* estimated RCOND may be too optimistic. To avoid misleading users, 715* the inverse condition number is set to 0.0 when the estimated 716* cwise error is at least CWISE_WRONG. 717* 718 CWISE_WRONG = SQRT( SLAMCH( 'Epsilon' ) ) 719 DO J = 1, NRHS 720 IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG ) 721 $ THEN 722 RCOND_TMP = CLA_GBRCOND_X( TRANS, N, KL, KU, AB, LDAB, 723 $ AFB, LDAFB, IPIV, X( 1, J ), INFO, WORK, RWORK ) 724 ELSE 725 RCOND_TMP = 0.0 726 END IF 727* 728* Cap the error at 1.0. 729* 730 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 731 $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0 ) 732 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0 733* 734* Threshold the error (see LAWN). 735* 736 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 737 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0 738 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0 739 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0 740 $ .AND. INFO.LT.N + J ) INFO = N + J 741 ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) 742 $ .LT. ERR_LBND ) THEN 743 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND 744 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0 745 END IF 746* 747* Save the condition number. 748* 749 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 750 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP 751 END IF 752 753 END DO 754 END IF 755* 756 RETURN 757* 758* End of CGBRFSX 759* 760 END 761