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