1*> \brief \b SLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLA_GBRFSX_EXTENDED + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_gbrfsx_extended.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_gbrfsx_extended.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_gbrfsx_extended.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 22* NRHS, AB, LDAB, AFB, LDAFB, IPIV, 23* COLEQU, C, B, LDB, Y, LDY, 24* BERR_OUT, N_NORMS, ERR_BNDS_NORM, 25* ERR_BNDS_COMP, RES, AYB, DY, 26* Y_TAIL, RCOND, ITHRESH, RTHRESH, 27* DZ_UB, IGNORE_CWISE, INFO ) 28* 29* .. Scalar Arguments .. 30* INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS, 31* $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH 32* LOGICAL COLEQU, IGNORE_CWISE 33* REAL RTHRESH, DZ_UB 34* .. 35* .. Array Arguments .. 36* INTEGER IPIV( * ) 37* REAL AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 38* $ Y( LDY, * ), RES(*), DY(*), Y_TAIL(*) 39* REAL C( * ), AYB(*), RCOND, BERR_OUT(*), 40* $ ERR_BNDS_NORM( NRHS, * ), 41* $ ERR_BNDS_COMP( NRHS, * ) 42* .. 43* 44* 45*> \par Purpose: 46* ============= 47*> 48*> \verbatim 49*> 50*> SLA_GBRFSX_EXTENDED improves the computed solution to a system of 51*> linear equations by performing extra-precise iterative refinement 52*> and provides error bounds and backward error estimates for the solution. 53*> This subroutine is called by SGBRFSX to perform iterative refinement. 54*> In addition to normwise error bound, the code provides maximum 55*> componentwise error bound if possible. See comments for ERR_BNDS_NORM 56*> and ERR_BNDS_COMP for details of the error bounds. Note that this 57*> subroutine is only resonsible for setting the second fields of 58*> ERR_BNDS_NORM and ERR_BNDS_COMP. 59*> \endverbatim 60* 61* Arguments: 62* ========== 63* 64*> \param[in] PREC_TYPE 65*> \verbatim 66*> PREC_TYPE is INTEGER 67*> Specifies the intermediate precision to be used in refinement. 68*> The value is defined by ILAPREC(P) where P is a CHARACTER and P 69*> = 'S': Single 70*> = 'D': Double 71*> = 'I': Indigenous 72*> = 'X' or 'E': Extra 73*> \endverbatim 74*> 75*> \param[in] TRANS_TYPE 76*> \verbatim 77*> TRANS_TYPE is INTEGER 78*> Specifies the transposition operation on A. 79*> The value is defined by ILATRANS(T) where T is a CHARACTER and T 80*> = 'N': No transpose 81*> = 'T': Transpose 82*> = 'C': Conjugate transpose 83*> \endverbatim 84*> 85*> \param[in] N 86*> \verbatim 87*> N is INTEGER 88*> The number of linear equations, i.e., the order of the 89*> matrix A. N >= 0. 90*> \endverbatim 91*> 92*> \param[in] KL 93*> \verbatim 94*> KL is INTEGER 95*> The number of subdiagonals within the band of A. KL >= 0. 96*> \endverbatim 97*> 98*> \param[in] KU 99*> \verbatim 100*> KU is INTEGER 101*> The number of superdiagonals within the band of A. KU >= 0 102*> \endverbatim 103*> 104*> \param[in] NRHS 105*> \verbatim 106*> NRHS is INTEGER 107*> The number of right-hand-sides, i.e., the number of columns of the 108*> matrix B. 109*> \endverbatim 110*> 111*> \param[in] AB 112*> \verbatim 113*> AB is REAL array, dimension (LDAB,N) 114*> On entry, the N-by-N matrix AB. 115*> \endverbatim 116*> 117*> \param[in] LDAB 118*> \verbatim 119*> LDAB is INTEGER 120*> The leading dimension of the array AB. LDAB >= max(1,N). 121*> \endverbatim 122*> 123*> \param[in] AFB 124*> \verbatim 125*> AFB is REAL array, dimension (LDAFB,N) 126*> The factors L and U from the factorization 127*> A = P*L*U as computed by SGBTRF. 128*> \endverbatim 129*> 130*> \param[in] LDAFB 131*> \verbatim 132*> LDAFB is INTEGER 133*> The leading dimension of the array AF. LDAFB >= max(1,N). 134*> \endverbatim 135*> 136*> \param[in] IPIV 137*> \verbatim 138*> IPIV is INTEGER array, dimension (N) 139*> The pivot indices from the factorization A = P*L*U 140*> as computed by SGBTRF; row i of the matrix was interchanged 141*> with row IPIV(i). 142*> \endverbatim 143*> 144*> \param[in] COLEQU 145*> \verbatim 146*> COLEQU is LOGICAL 147*> If .TRUE. then column equilibration was done to A before calling 148*> this routine. This is needed to compute the solution and error 149*> bounds correctly. 150*> \endverbatim 151*> 152*> \param[in] C 153*> \verbatim 154*> C is REAL array, dimension (N) 155*> The column scale factors for A. If COLEQU = .FALSE., C 156*> is not accessed. If C is input, each element of C should be a power 157*> of the radix to ensure a reliable solution and error estimates. 158*> Scaling by powers of the radix does not cause rounding errors unless 159*> the result underflows or overflows. Rounding errors during scaling 160*> lead to refining with a matrix that is not equivalent to the 161*> input matrix, producing error estimates that may not be 162*> reliable. 163*> \endverbatim 164*> 165*> \param[in] B 166*> \verbatim 167*> B is REAL array, dimension (LDB,NRHS) 168*> The right-hand-side matrix B. 169*> \endverbatim 170*> 171*> \param[in] LDB 172*> \verbatim 173*> LDB is INTEGER 174*> The leading dimension of the array B. LDB >= max(1,N). 175*> \endverbatim 176*> 177*> \param[in,out] Y 178*> \verbatim 179*> Y is REAL array, dimension (LDY,NRHS) 180*> On entry, the solution matrix X, as computed by SGBTRS. 181*> On exit, the improved solution matrix Y. 182*> \endverbatim 183*> 184*> \param[in] LDY 185*> \verbatim 186*> LDY is INTEGER 187*> The leading dimension of the array Y. LDY >= max(1,N). 188*> \endverbatim 189*> 190*> \param[out] BERR_OUT 191*> \verbatim 192*> BERR_OUT is REAL array, dimension (NRHS) 193*> On exit, BERR_OUT(j) contains the componentwise relative backward 194*> error for right-hand-side j from the formula 195*> max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 196*> where abs(Z) is the componentwise absolute value of the matrix 197*> or vector Z. This is computed by SLA_LIN_BERR. 198*> \endverbatim 199*> 200*> \param[in] N_NORMS 201*> \verbatim 202*> N_NORMS is INTEGER 203*> Determines which error bounds to return (see ERR_BNDS_NORM 204*> and ERR_BNDS_COMP). 205*> If N_NORMS >= 1 return normwise error bounds. 206*> If N_NORMS >= 2 return componentwise error bounds. 207*> \endverbatim 208*> 209*> \param[in,out] ERR_BNDS_NORM 210*> \verbatim 211*> ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS) 212*> For each right-hand side, this array contains information about 213*> various error bounds and condition numbers corresponding to the 214*> normwise relative error, which is defined as follows: 215*> 216*> Normwise relative error in the ith solution vector: 217*> max_j (abs(XTRUE(j,i) - X(j,i))) 218*> ------------------------------ 219*> max_j abs(X(j,i)) 220*> 221*> The array is indexed by the type of error information as described 222*> below. There currently are up to three pieces of information 223*> returned. 224*> 225*> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith 226*> right-hand side. 227*> 228*> The second index in ERR_BNDS_NORM(:,err) contains the following 229*> three fields: 230*> err = 1 "Trust/don't trust" boolean. Trust the answer if the 231*> reciprocal condition number is less than the threshold 232*> sqrt(n) * slamch('Epsilon'). 233*> 234*> err = 2 "Guaranteed" error bound: The estimated forward error, 235*> almost certainly within a factor of 10 of the true error 236*> so long as the next entry is greater than the threshold 237*> sqrt(n) * slamch('Epsilon'). This error bound should only 238*> be trusted if the previous boolean is true. 239*> 240*> err = 3 Reciprocal condition number: Estimated normwise 241*> reciprocal condition number. Compared with the threshold 242*> sqrt(n) * slamch('Epsilon') to determine if the error 243*> estimate is "guaranteed". These reciprocal condition 244*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 245*> appropriately scaled matrix Z. 246*> Let Z = S*A, where S scales each row by a power of the 247*> radix so all absolute row sums of Z are approximately 1. 248*> 249*> This subroutine is only responsible for setting the second field 250*> above. 251*> See Lapack Working Note 165 for further details and extra 252*> cautions. 253*> \endverbatim 254*> 255*> \param[in,out] ERR_BNDS_COMP 256*> \verbatim 257*> ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS) 258*> For each right-hand side, this array contains information about 259*> various error bounds and condition numbers corresponding to the 260*> componentwise relative error, which is defined as follows: 261*> 262*> Componentwise relative error in the ith solution vector: 263*> abs(XTRUE(j,i) - X(j,i)) 264*> max_j ---------------------- 265*> abs(X(j,i)) 266*> 267*> The array is indexed by the right-hand side i (on which the 268*> componentwise relative error depends), and the type of error 269*> information as described below. There currently are up to three 270*> pieces of information returned for each right-hand side. If 271*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then 272*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most 273*> the first (:,N_ERR_BNDS) entries are returned. 274*> 275*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith 276*> right-hand side. 277*> 278*> The second index in ERR_BNDS_COMP(:,err) contains the following 279*> three fields: 280*> err = 1 "Trust/don't trust" boolean. Trust the answer if the 281*> reciprocal condition number is less than the threshold 282*> sqrt(n) * slamch('Epsilon'). 283*> 284*> err = 2 "Guaranteed" error bound: The estimated forward error, 285*> almost certainly within a factor of 10 of the true error 286*> so long as the next entry is greater than the threshold 287*> sqrt(n) * slamch('Epsilon'). This error bound should only 288*> be trusted if the previous boolean is true. 289*> 290*> err = 3 Reciprocal condition number: Estimated componentwise 291*> reciprocal condition number. Compared with the threshold 292*> sqrt(n) * slamch('Epsilon') to determine if the error 293*> estimate is "guaranteed". These reciprocal condition 294*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 295*> appropriately scaled matrix Z. 296*> Let Z = S*(A*diag(x)), where x is the solution for the 297*> current right-hand side and S scales each row of 298*> A*diag(x) by a power of the radix so all absolute row 299*> sums of Z are approximately 1. 300*> 301*> This subroutine is only responsible for setting the second field 302*> above. 303*> See Lapack Working Note 165 for further details and extra 304*> cautions. 305*> \endverbatim 306*> 307*> \param[in] RES 308*> \verbatim 309*> RES is REAL array, dimension (N) 310*> Workspace to hold the intermediate residual. 311*> \endverbatim 312*> 313*> \param[in] AYB 314*> \verbatim 315*> AYB is REAL array, dimension (N) 316*> Workspace. This can be the same workspace passed for Y_TAIL. 317*> \endverbatim 318*> 319*> \param[in] DY 320*> \verbatim 321*> DY is REAL array, dimension (N) 322*> Workspace to hold the intermediate solution. 323*> \endverbatim 324*> 325*> \param[in] Y_TAIL 326*> \verbatim 327*> Y_TAIL is REAL array, dimension (N) 328*> Workspace to hold the trailing bits of the intermediate solution. 329*> \endverbatim 330*> 331*> \param[in] RCOND 332*> \verbatim 333*> RCOND is REAL 334*> Reciprocal scaled condition number. This is an estimate of the 335*> reciprocal Skeel condition number of the matrix A after 336*> equilibration (if done). If this is less than the machine 337*> precision (in particular, if it is zero), the matrix is singular 338*> to working precision. Note that the error may still be small even 339*> if this number is very small and the matrix appears ill- 340*> conditioned. 341*> \endverbatim 342*> 343*> \param[in] ITHRESH 344*> \verbatim 345*> ITHRESH is INTEGER 346*> The maximum number of residual computations allowed for 347*> refinement. The default is 10. For 'aggressive' set to 100 to 348*> permit convergence using approximate factorizations or 349*> factorizations other than LU. If the factorization uses a 350*> technique other than Gaussian elimination, the guarantees in 351*> ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy. 352*> \endverbatim 353*> 354*> \param[in] RTHRESH 355*> \verbatim 356*> RTHRESH is REAL 357*> Determines when to stop refinement if the error estimate stops 358*> decreasing. Refinement will stop when the next solution no longer 359*> satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is 360*> the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The 361*> default value is 0.5. For 'aggressive' set to 0.9 to permit 362*> convergence on extremely ill-conditioned matrices. See LAWN 165 363*> for more details. 364*> \endverbatim 365*> 366*> \param[in] DZ_UB 367*> \verbatim 368*> DZ_UB is REAL 369*> Determines when to start considering componentwise convergence. 370*> Componentwise convergence is only considered after each component 371*> of the solution Y is stable, which we define as the relative 372*> change in each component being less than DZ_UB. The default value 373*> is 0.25, requiring the first bit to be stable. See LAWN 165 for 374*> more details. 375*> \endverbatim 376*> 377*> \param[in] IGNORE_CWISE 378*> \verbatim 379*> IGNORE_CWISE is LOGICAL 380*> If .TRUE. then ignore componentwise convergence. Default value 381*> is .FALSE.. 382*> \endverbatim 383*> 384*> \param[out] INFO 385*> \verbatim 386*> INFO is INTEGER 387*> = 0: Successful exit. 388*> < 0: if INFO = -i, the ith argument to SGBTRS had an illegal 389*> value 390*> \endverbatim 391* 392* Authors: 393* ======== 394* 395*> \author Univ. of Tennessee 396*> \author Univ. of California Berkeley 397*> \author Univ. of Colorado Denver 398*> \author NAG Ltd. 399* 400*> \ingroup realGBcomputational 401* 402* ===================================================================== 403 SUBROUTINE SLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 404 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, 405 $ COLEQU, C, B, LDB, Y, LDY, 406 $ BERR_OUT, N_NORMS, ERR_BNDS_NORM, 407 $ ERR_BNDS_COMP, RES, AYB, DY, 408 $ Y_TAIL, RCOND, ITHRESH, RTHRESH, 409 $ DZ_UB, IGNORE_CWISE, INFO ) 410* 411* -- LAPACK computational routine -- 412* -- LAPACK is a software package provided by Univ. of Tennessee, -- 413* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 414* 415* .. Scalar Arguments .. 416 INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS, 417 $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH 418 LOGICAL COLEQU, IGNORE_CWISE 419 REAL RTHRESH, DZ_UB 420* .. 421* .. Array Arguments .. 422 INTEGER IPIV( * ) 423 REAL AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 424 $ Y( LDY, * ), RES(*), DY(*), Y_TAIL(*) 425 REAL C( * ), AYB(*), RCOND, BERR_OUT(*), 426 $ ERR_BNDS_NORM( NRHS, * ), 427 $ ERR_BNDS_COMP( NRHS, * ) 428* .. 429* 430* ===================================================================== 431* 432* .. Local Scalars .. 433 CHARACTER TRANS 434 INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE 435 REAL YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT, 436 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX, 437 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z, 438 $ EPS, HUGEVAL, INCR_THRESH 439 LOGICAL INCR_PREC 440* .. 441* .. Parameters .. 442 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE, 443 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL, 444 $ EXTRA_Y 445 PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1, 446 $ CONV_STATE = 2, NOPROG_STATE = 3 ) 447 PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1, 448 $ EXTRA_Y = 2 ) 449 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 450 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 451 INTEGER CMP_ERR_I, PIV_GROWTH_I 452 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 453 $ BERR_I = 3 ) 454 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 455 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8, 456 $ PIV_GROWTH_I = 9 ) 457 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 458 $ LA_LINRX_CWISE_I 459 PARAMETER ( LA_LINRX_ITREF_I = 1, 460 $ LA_LINRX_ITHRESH_I = 2 ) 461 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 462 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 463 $ LA_LINRX_RCOND_I 464 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 465 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 466* .. 467* .. External Subroutines .. 468 EXTERNAL SAXPY, SCOPY, SGBTRS, SGBMV, BLAS_SGBMV_X, 469 $ BLAS_SGBMV2_X, SLA_GBAMV, SLA_WWADDW, SLAMCH, 470 $ CHLA_TRANSTYPE, SLA_LIN_BERR 471 REAL SLAMCH 472 CHARACTER CHLA_TRANSTYPE 473* .. 474* .. Intrinsic Functions .. 475 INTRINSIC ABS, MAX, MIN 476* .. 477* .. Executable Statements .. 478* 479 IF (INFO.NE.0) RETURN 480 TRANS = CHLA_TRANSTYPE(TRANS_TYPE) 481 EPS = SLAMCH( 'Epsilon' ) 482 HUGEVAL = SLAMCH( 'Overflow' ) 483* Force HUGEVAL to Inf 484 HUGEVAL = HUGEVAL * HUGEVAL 485* Using HUGEVAL may lead to spurious underflows. 486 INCR_THRESH = REAL( N ) * EPS 487 M = KL+KU+1 488 489 DO J = 1, NRHS 490 Y_PREC_STATE = EXTRA_RESIDUAL 491 IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN 492 DO I = 1, N 493 Y_TAIL( I ) = 0.0 494 END DO 495 END IF 496 497 DXRAT = 0.0 498 DXRATMAX = 0.0 499 DZRAT = 0.0 500 DZRATMAX = 0.0 501 FINAL_DX_X = HUGEVAL 502 FINAL_DZ_Z = HUGEVAL 503 PREVNORMDX = HUGEVAL 504 PREV_DZ_Z = HUGEVAL 505 DZ_Z = HUGEVAL 506 DX_X = HUGEVAL 507 508 X_STATE = WORKING_STATE 509 Z_STATE = UNSTABLE_STATE 510 INCR_PREC = .FALSE. 511 512 DO CNT = 1, ITHRESH 513* 514* Compute residual RES = B_s - op(A_s) * Y, 515* op(A) = A, A**T, or A**H depending on TRANS (and type). 516* 517 CALL SCOPY( N, B( 1, J ), 1, RES, 1 ) 518 IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN 519 CALL SGBMV( TRANS, M, N, KL, KU, -1.0, AB, LDAB, 520 $ Y( 1, J ), 1, 1.0, RES, 1 ) 521 ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN 522 CALL BLAS_SGBMV_X( TRANS_TYPE, N, N, KL, KU, 523 $ -1.0, AB, LDAB, Y( 1, J ), 1, 1.0, RES, 1, 524 $ PREC_TYPE ) 525 ELSE 526 CALL BLAS_SGBMV2_X( TRANS_TYPE, N, N, KL, KU, -1.0, 527 $ AB, LDAB, Y( 1, J ), Y_TAIL, 1, 1.0, RES, 1, 528 $ PREC_TYPE ) 529 END IF 530 531! XXX: RES is no longer needed. 532 CALL SCOPY( N, RES, 1, DY, 1 ) 533 CALL SGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, DY, N, 534 $ INFO ) 535* 536* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. 537* 538 NORMX = 0.0 539 NORMY = 0.0 540 NORMDX = 0.0 541 DZ_Z = 0.0 542 YMIN = HUGEVAL 543 544 DO I = 1, N 545 YK = ABS( Y( I, J ) ) 546 DYK = ABS( DY( I ) ) 547 548 IF ( YK .NE. 0.0 ) THEN 549 DZ_Z = MAX( DZ_Z, DYK / YK ) 550 ELSE IF ( DYK .NE. 0.0 ) THEN 551 DZ_Z = HUGEVAL 552 END IF 553 554 YMIN = MIN( YMIN, YK ) 555 556 NORMY = MAX( NORMY, YK ) 557 558 IF ( COLEQU ) THEN 559 NORMX = MAX( NORMX, YK * C( I ) ) 560 NORMDX = MAX( NORMDX, DYK * C( I ) ) 561 ELSE 562 NORMX = NORMY 563 NORMDX = MAX( NORMDX, DYK ) 564 END IF 565 END DO 566 567 IF ( NORMX .NE. 0.0 ) THEN 568 DX_X = NORMDX / NORMX 569 ELSE IF ( NORMDX .EQ. 0.0 ) THEN 570 DX_X = 0.0 571 ELSE 572 DX_X = HUGEVAL 573 END IF 574 575 DXRAT = NORMDX / PREVNORMDX 576 DZRAT = DZ_Z / PREV_DZ_Z 577* 578* Check termination criteria. 579* 580 IF ( .NOT.IGNORE_CWISE 581 $ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY 582 $ .AND. Y_PREC_STATE .LT. EXTRA_Y ) 583 $ INCR_PREC = .TRUE. 584 585 IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH ) 586 $ X_STATE = WORKING_STATE 587 IF ( X_STATE .EQ. WORKING_STATE ) THEN 588 IF ( DX_X .LE. EPS ) THEN 589 X_STATE = CONV_STATE 590 ELSE IF ( DXRAT .GT. RTHRESH ) THEN 591 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 592 INCR_PREC = .TRUE. 593 ELSE 594 X_STATE = NOPROG_STATE 595 END IF 596 ELSE 597 IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT 598 END IF 599 IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X 600 END IF 601 602 IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB ) 603 $ Z_STATE = WORKING_STATE 604 IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH ) 605 $ Z_STATE = WORKING_STATE 606 IF ( Z_STATE .EQ. WORKING_STATE ) THEN 607 IF ( DZ_Z .LE. EPS ) THEN 608 Z_STATE = CONV_STATE 609 ELSE IF ( DZ_Z .GT. DZ_UB ) THEN 610 Z_STATE = UNSTABLE_STATE 611 DZRATMAX = 0.0 612 FINAL_DZ_Z = HUGEVAL 613 ELSE IF ( DZRAT .GT. RTHRESH ) THEN 614 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 615 INCR_PREC = .TRUE. 616 ELSE 617 Z_STATE = NOPROG_STATE 618 END IF 619 ELSE 620 IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT 621 END IF 622 IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 623 END IF 624* 625* Exit if both normwise and componentwise stopped working, 626* but if componentwise is unstable, let it go at least two 627* iterations. 628* 629 IF ( X_STATE.NE.WORKING_STATE ) THEN 630 IF ( IGNORE_CWISE ) GOTO 666 631 IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE ) 632 $ GOTO 666 633 IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666 634 END IF 635 636 IF ( INCR_PREC ) THEN 637 INCR_PREC = .FALSE. 638 Y_PREC_STATE = Y_PREC_STATE + 1 639 DO I = 1, N 640 Y_TAIL( I ) = 0.0 641 END DO 642 END IF 643 644 PREVNORMDX = NORMDX 645 PREV_DZ_Z = DZ_Z 646* 647* Update soluton. 648* 649 IF (Y_PREC_STATE .LT. EXTRA_Y) THEN 650 CALL SAXPY( N, 1.0, DY, 1, Y(1,J), 1 ) 651 ELSE 652 CALL SLA_WWADDW( N, Y(1,J), Y_TAIL, DY ) 653 END IF 654 655 END DO 656* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT. 657 666 CONTINUE 658* 659* Set final_* when cnt hits ithresh. 660* 661 IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X 662 IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 663* 664* Compute error bounds. 665* 666 IF ( N_NORMS .GE. 1 ) THEN 667 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 668 $ FINAL_DX_X / (1 - DXRATMAX) 669 END IF 670 IF (N_NORMS .GE. 2) THEN 671 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 672 $ FINAL_DZ_Z / (1 - DZRATMAX) 673 END IF 674* 675* Compute componentwise relative backward error from formula 676* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 677* where abs(Z) is the componentwise absolute value of the matrix 678* or vector Z. 679* 680* Compute residual RES = B_s - op(A_s) * Y, 681* op(A) = A, A**T, or A**H depending on TRANS (and type). 682* 683 CALL SCOPY( N, B( 1, J ), 1, RES, 1 ) 684 CALL SGBMV(TRANS, N, N, KL, KU, -1.0, AB, LDAB, Y(1,J), 685 $ 1, 1.0, RES, 1 ) 686 687 DO I = 1, N 688 AYB( I ) = ABS( B( I, J ) ) 689 END DO 690* 691* Compute abs(op(A_s))*abs(Y) + abs(B_s). 692* 693 CALL SLA_GBAMV( TRANS_TYPE, N, N, KL, KU, 1.0, 694 $ AB, LDAB, Y(1, J), 1, 1.0, AYB, 1 ) 695 696 CALL SLA_LIN_BERR( N, N, 1, RES, AYB, BERR_OUT( J ) ) 697* 698* End of loop for each RHS 699* 700 END DO 701* 702 RETURN 703* 704* End of SLA_GBRFSX_EXTENDED 705* 706 END 707