1*> \brief \b DGBRFSX 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DGBRFSX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbrfsx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbrfsx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbrfsx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGBRFSX( 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, IWORK, 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( * ), IWORK( * ) 35* DOUBLE PRECISION 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, * ) 40* .. 41* 42* 43*> \par Purpose: 44* ============= 45*> 46*> \verbatim 47*> 48*> DGBRFSX 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 = 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (4*N) 390*> \endverbatim 391*> 392*> \param[out] IWORK 393*> \verbatim 394*> IWORK is INTEGER array, dimension (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 doubleGBcomputational 435* 436* ===================================================================== 437 SUBROUTINE DGBRFSX( 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, IWORK, 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( * ), IWORK( * ) 456 DOUBLE PRECISION 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, * ) 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 489 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE 490 INTEGER N_NORMS 491 DOUBLE PRECISION ANORM, RCOND_TMP 492 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG 493 LOGICAL IGNORE_CWISE 494 INTEGER ITHRESH 495 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH 496* .. 497* .. External Subroutines .. 498 EXTERNAL XERBLA, DGBCON 499 EXTERNAL DLA_GBRFSX_EXTENDED 500* .. 501* .. Intrinsic Functions .. 502 INTRINSIC MAX, SQRT 503* .. 504* .. External Functions .. 505 EXTERNAL LSAME, ILATRANS, ILAPREC 506 EXTERNAL DLAMCH, DLANGB, DLA_GBRCOND 507 DOUBLE PRECISION DLAMCH, DLANGB, DLA_GBRCOND 508 LOGICAL LSAME 509 INTEGER ILATRANS, ILAPREC 510* .. 511* .. Executable Statements .. 512* 513* Check the input parameters. 514* 515 INFO = 0 516 TRANS_TYPE = ILATRANS( TRANS ) 517 REF_TYPE = INT( ITREF_DEFAULT ) 518 IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN 519 IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN 520 PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT 521 ELSE 522 REF_TYPE = PARAMS( LA_LINRX_ITREF_I ) 523 END IF 524 END IF 525* 526* Set default parameters. 527* 528 ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' ) 529 ITHRESH = INT( ITHRESH_DEFAULT ) 530 RTHRESH = RTHRESH_DEFAULT 531 UNSTABLE_THRESH = DZTHRESH_DEFAULT 532 IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0 533* 534 IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN 535 IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN 536 PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH 537 ELSE 538 ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) ) 539 END IF 540 END IF 541 IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN 542 IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN 543 IF ( IGNORE_CWISE ) THEN 544 PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0 545 ELSE 546 PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0 547 END IF 548 ELSE 549 IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0 550 END IF 551 END IF 552 IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN 553 N_NORMS = 0 554 ELSE IF ( IGNORE_CWISE ) THEN 555 N_NORMS = 1 556 ELSE 557 N_NORMS = 2 558 END IF 559* 560 NOTRAN = LSAME( TRANS, 'N' ) 561 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' ) 562 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' ) 563* 564* Test input parameters. 565* 566 IF( TRANS_TYPE.EQ.-1 ) THEN 567 INFO = -1 568 ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND. 569 $ .NOT.LSAME( EQUED, 'N' ) ) THEN 570 INFO = -2 571 ELSE IF( N.LT.0 ) THEN 572 INFO = -3 573 ELSE IF( KL.LT.0 ) THEN 574 INFO = -4 575 ELSE IF( KU.LT.0 ) THEN 576 INFO = -5 577 ELSE IF( NRHS.LT.0 ) THEN 578 INFO = -6 579 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 580 INFO = -8 581 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 582 INFO = -10 583 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 584 INFO = -13 585 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 586 INFO = -15 587 END IF 588 IF( INFO.NE.0 ) THEN 589 CALL XERBLA( 'DGBRFSX', -INFO ) 590 RETURN 591 END IF 592* 593* Quick return if possible. 594* 595 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 596 RCOND = 1.0D+0 597 DO J = 1, NRHS 598 BERR( J ) = 0.0D+0 599 IF ( N_ERR_BNDS .GE. 1 ) THEN 600 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 601 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 602 END IF 603 IF ( N_ERR_BNDS .GE. 2 ) THEN 604 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0 605 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0 606 END IF 607 IF ( N_ERR_BNDS .GE. 3 ) THEN 608 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0 609 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0 610 END IF 611 END DO 612 RETURN 613 END IF 614* 615* Default to failure. 616* 617 RCOND = 0.0D+0 618 DO J = 1, NRHS 619 BERR( J ) = 1.0D+0 620 IF ( N_ERR_BNDS .GE. 1 ) THEN 621 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 622 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 623 END IF 624 IF ( N_ERR_BNDS .GE. 2 ) THEN 625 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 626 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 627 END IF 628 IF ( N_ERR_BNDS .GE. 3 ) THEN 629 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0 630 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0 631 END IF 632 END DO 633* 634* Compute the norm of A and the reciprocal of the condition 635* number of A. 636* 637 IF( NOTRAN ) THEN 638 NORM = 'I' 639 ELSE 640 NORM = '1' 641 END IF 642 ANORM = DLANGB( NORM, N, KL, KU, AB, LDAB, WORK ) 643 CALL DGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND, 644 $ WORK, IWORK, INFO ) 645* 646* Perform refinement on each right-hand side 647* 648 IF ( REF_TYPE .NE. 0 .AND. INFO .EQ. 0 ) THEN 649 650 PREC_TYPE = ILAPREC( 'E' ) 651 652 IF ( NOTRAN ) THEN 653 CALL DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 654 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B, 655 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, 656 $ ERR_BNDS_COMP, WORK( N+1 ), WORK( 1 ), WORK( 2*N+1 ), 657 $ WORK( 1 ), RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, 658 $ IGNORE_CWISE, INFO ) 659 ELSE 660 CALL DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 661 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, ROWEQU, R, B, 662 $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, 663 $ ERR_BNDS_COMP, WORK( N+1 ), WORK( 1 ), WORK( 2*N+1 ), 664 $ WORK( 1 ), RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, 665 $ IGNORE_CWISE, INFO ) 666 END IF 667 END IF 668 669 ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' ) 670 IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN 671* 672* Compute scaled normwise condition number cond(A*C). 673* 674 IF ( COLEQU .AND. NOTRAN ) THEN 675 RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, 676 $ LDAFB, IPIV, -1, C, INFO, WORK, IWORK ) 677 ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN 678 RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, 679 $ LDAFB, IPIV, -1, R, INFO, WORK, IWORK ) 680 ELSE 681 RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, 682 $ LDAFB, IPIV, 0, R, INFO, WORK, IWORK ) 683 END IF 684 DO J = 1, NRHS 685* 686* Cap the error at 1.0. 687* 688 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 689 $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 690 $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 691* 692* Threshold the error (see LAWN). 693* 694 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 695 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0 696 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0 697 IF ( INFO .LE. N ) INFO = N + J 698 ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND ) 699 $ THEN 700 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND 701 ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0 702 END IF 703* 704* Save the condition number. 705* 706 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 707 ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP 708 END IF 709 710 END DO 711 END IF 712 713 IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN 714* 715* Compute componentwise condition number cond(A*diag(Y(:,J))) for 716* each right-hand side using the current solution as an estimate of 717* the true solution. If the componentwise error estimate is too 718* large, then the solution is a lousy estimate of truth and the 719* estimated RCOND may be too optimistic. To avoid misleading users, 720* the inverse condition number is set to 0.0 when the estimated 721* cwise error is at least CWISE_WRONG. 722* 723 CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) ) 724 DO J = 1, NRHS 725 IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG ) 726 $ THEN 727 RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB, 728 $ LDAFB, IPIV, 1, X( 1, J ), INFO, WORK, IWORK ) 729 ELSE 730 RCOND_TMP = 0.0D+0 731 END IF 732* 733* Cap the error at 1.0. 734* 735 IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I 736 $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 ) 737 $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 738* 739* Threshold the error (see LAWN). 740* 741 IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN 742 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0 743 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0 744 IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0 745 $ .AND. INFO.LT.N + J ) INFO = N + J 746 ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) 747 $ .LT. ERR_LBND ) THEN 748 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND 749 ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0 750 END IF 751* 752* Save the condition number. 753* 754 IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN 755 ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP 756 END IF 757 758 END DO 759 END IF 760* 761 RETURN 762* 763* End of DGBRFSX 764* 765 END 766