1*> \brief \b CLA_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 CLA_GBRFSX_EXTENDED + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_gbrfsx_extended.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_gbrfsx_extended.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_gbrfsx_extended.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CLA_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* COMPLEX 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*> CLA_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 CGBRFSX 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 69*> P = 'S': Single 70*> = 'D': Double 71*> = 'I': Indigenous 72*> = 'X', '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 80*> T = '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 COMPLEX 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 COMPLEX array, dimension (LDAF,N) 126*> The factors L and U from the factorization 127*> A = P*L*U as computed by CGBTRF. 128*> \endverbatim 129*> 130*> \param[in] LDAFB 131*> \verbatim 132*> LDAFB is INTEGER 133*> The leading dimension of the array AF. LDAF >= 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 CGBTRF; 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 COMPLEX 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 COMPLEX array, dimension (LDY,NRHS) 180*> On entry, the solution matrix X, as computed by CGBTRS. 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 CLA_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 212*> (NRHS, N_ERR_BNDS) 213*> For each right-hand side, this array contains information about 214*> various error bounds and condition numbers corresponding to the 215*> normwise relative error, which is defined as follows: 216*> 217*> Normwise relative error in the ith solution vector: 218*> max_j (abs(XTRUE(j,i) - X(j,i))) 219*> ------------------------------ 220*> max_j abs(X(j,i)) 221*> 222*> The array is indexed by the type of error information as described 223*> below. There currently are up to three pieces of information 224*> returned. 225*> 226*> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith 227*> right-hand side. 228*> 229*> The second index in ERR_BNDS_NORM(:,err) contains the following 230*> three fields: 231*> err = 1 "Trust/don't trust" boolean. Trust the answer if the 232*> reciprocal condition number is less than the threshold 233*> sqrt(n) * slamch('Epsilon'). 234*> 235*> err = 2 "Guaranteed" error bound: The estimated forward error, 236*> almost certainly within a factor of 10 of the true error 237*> so long as the next entry is greater than the threshold 238*> sqrt(n) * slamch('Epsilon'). This error bound should only 239*> be trusted if the previous boolean is true. 240*> 241*> err = 3 Reciprocal condition number: Estimated normwise 242*> reciprocal condition number. Compared with the threshold 243*> sqrt(n) * slamch('Epsilon') to determine if the error 244*> estimate is "guaranteed". These reciprocal condition 245*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 246*> appropriately scaled matrix Z. 247*> Let Z = S*A, where S scales each row by a power of the 248*> radix so all absolute row sums of Z are approximately 1. 249*> 250*> This subroutine is only responsible for setting the second field 251*> above. 252*> See Lapack Working Note 165 for further details and extra 253*> cautions. 254*> \endverbatim 255*> 256*> \param[in,out] ERR_BNDS_COMP 257*> \verbatim 258*> ERR_BNDS_COMP is REAL array, dimension 259*> (NRHS, N_ERR_BNDS) 260*> For each right-hand side, this array contains information about 261*> various error bounds and condition numbers corresponding to the 262*> componentwise relative error, which is defined as follows: 263*> 264*> Componentwise relative error in the ith solution vector: 265*> abs(XTRUE(j,i) - X(j,i)) 266*> max_j ---------------------- 267*> abs(X(j,i)) 268*> 269*> The array is indexed by the right-hand side i (on which the 270*> componentwise relative error depends), and the type of error 271*> information as described below. There currently are up to three 272*> pieces of information returned for each right-hand side. If 273*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then 274*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most 275*> the first (:,N_ERR_BNDS) entries are returned. 276*> 277*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith 278*> right-hand side. 279*> 280*> The second index in ERR_BNDS_COMP(:,err) contains the following 281*> three fields: 282*> err = 1 "Trust/don't trust" boolean. Trust the answer if the 283*> reciprocal condition number is less than the threshold 284*> sqrt(n) * slamch('Epsilon'). 285*> 286*> err = 2 "Guaranteed" error bound: The estimated forward error, 287*> almost certainly within a factor of 10 of the true error 288*> so long as the next entry is greater than the threshold 289*> sqrt(n) * slamch('Epsilon'). This error bound should only 290*> be trusted if the previous boolean is true. 291*> 292*> err = 3 Reciprocal condition number: Estimated componentwise 293*> reciprocal condition number. Compared with the threshold 294*> sqrt(n) * slamch('Epsilon') to determine if the error 295*> estimate is "guaranteed". These reciprocal condition 296*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 297*> appropriately scaled matrix Z. 298*> Let Z = S*(A*diag(x)), where x is the solution for the 299*> current right-hand side and S scales each row of 300*> A*diag(x) by a power of the radix so all absolute row 301*> sums of Z are approximately 1. 302*> 303*> This subroutine is only responsible for setting the second field 304*> above. 305*> See Lapack Working Note 165 for further details and extra 306*> cautions. 307*> \endverbatim 308*> 309*> \param[in] RES 310*> \verbatim 311*> RES is COMPLEX array, dimension (N) 312*> Workspace to hold the intermediate residual. 313*> \endverbatim 314*> 315*> \param[in] AYB 316*> \verbatim 317*> AYB is REAL array, dimension (N) 318*> Workspace. 319*> \endverbatim 320*> 321*> \param[in] DY 322*> \verbatim 323*> DY is COMPLEX array, dimension (N) 324*> Workspace to hold the intermediate solution. 325*> \endverbatim 326*> 327*> \param[in] Y_TAIL 328*> \verbatim 329*> Y_TAIL is COMPLEX array, dimension (N) 330*> Workspace to hold the trailing bits of the intermediate solution. 331*> \endverbatim 332*> 333*> \param[in] RCOND 334*> \verbatim 335*> RCOND is REAL 336*> Reciprocal scaled condition number. This is an estimate of the 337*> reciprocal Skeel condition number of the matrix A after 338*> equilibration (if done). If this is less than the machine 339*> precision (in particular, if it is zero), the matrix is singular 340*> to working precision. Note that the error may still be small even 341*> if this number is very small and the matrix appears ill- 342*> conditioned. 343*> \endverbatim 344*> 345*> \param[in] ITHRESH 346*> \verbatim 347*> ITHRESH is INTEGER 348*> The maximum number of residual computations allowed for 349*> refinement. The default is 10. For 'aggressive' set to 100 to 350*> permit convergence using approximate factorizations or 351*> factorizations other than LU. If the factorization uses a 352*> technique other than Gaussian elimination, the guarantees in 353*> ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy. 354*> \endverbatim 355*> 356*> \param[in] RTHRESH 357*> \verbatim 358*> RTHRESH is REAL 359*> Determines when to stop refinement if the error estimate stops 360*> decreasing. Refinement will stop when the next solution no longer 361*> satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is 362*> the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The 363*> default value is 0.5. For 'aggressive' set to 0.9 to permit 364*> convergence on extremely ill-conditioned matrices. See LAWN 165 365*> for more details. 366*> \endverbatim 367*> 368*> \param[in] DZ_UB 369*> \verbatim 370*> DZ_UB is REAL 371*> Determines when to start considering componentwise convergence. 372*> Componentwise convergence is only considered after each component 373*> of the solution Y is stable, which we definte as the relative 374*> change in each component being less than DZ_UB. The default value 375*> is 0.25, requiring the first bit to be stable. See LAWN 165 for 376*> more details. 377*> \endverbatim 378*> 379*> \param[in] IGNORE_CWISE 380*> \verbatim 381*> IGNORE_CWISE is LOGICAL 382*> If .TRUE. then ignore componentwise convergence. Default value 383*> is .FALSE.. 384*> \endverbatim 385*> 386*> \param[out] INFO 387*> \verbatim 388*> INFO is INTEGER 389*> = 0: Successful exit. 390*> < 0: if INFO = -i, the ith argument to CGBTRS had an illegal 391*> value 392*> \endverbatim 393* 394* Authors: 395* ======== 396* 397*> \author Univ. of Tennessee 398*> \author Univ. of California Berkeley 399*> \author Univ. of Colorado Denver 400*> \author NAG Ltd. 401* 402*> \date September 2012 403* 404*> \ingroup complexGBcomputational 405* 406* ===================================================================== 407 SUBROUTINE CLA_GBRFSX_EXTENDED ( PREC_TYPE, TRANS_TYPE, N, KL, KU, 408 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, 409 $ COLEQU, C, B, LDB, Y, LDY, 410 $ BERR_OUT, N_NORMS, ERR_BNDS_NORM, 411 $ ERR_BNDS_COMP, RES, AYB, DY, 412 $ Y_TAIL, RCOND, ITHRESH, RTHRESH, 413 $ DZ_UB, IGNORE_CWISE, INFO ) 414* 415* -- LAPACK computational routine (version 3.4.2) -- 416* -- LAPACK is a software package provided by Univ. of Tennessee, -- 417* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 418* September 2012 419* 420* .. Scalar Arguments .. 421 INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS, 422 $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH 423 LOGICAL COLEQU, IGNORE_CWISE 424 REAL RTHRESH, DZ_UB 425* .. 426* .. Array Arguments .. 427 INTEGER IPIV( * ) 428 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 429 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * ) 430 REAL C( * ), AYB(*), RCOND, BERR_OUT( * ), 431 $ ERR_BNDS_NORM( NRHS, * ), 432 $ ERR_BNDS_COMP( NRHS, * ) 433* .. 434* 435* ===================================================================== 436* 437* .. Local Scalars .. 438 CHARACTER TRANS 439 INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE 440 REAL YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT, 441 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX, 442 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z, 443 $ EPS, HUGEVAL, INCR_THRESH 444 LOGICAL INCR_PREC 445 COMPLEX ZDUM 446* .. 447* .. Parameters .. 448 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE, 449 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL, 450 $ EXTRA_Y 451 PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1, 452 $ CONV_STATE = 2, NOPROG_STATE = 3 ) 453 PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1, 454 $ EXTRA_Y = 2 ) 455 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 456 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 457 INTEGER CMP_ERR_I, PIV_GROWTH_I 458 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 459 $ BERR_I = 3 ) 460 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 461 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8, 462 $ PIV_GROWTH_I = 9 ) 463 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 464 $ LA_LINRX_CWISE_I 465 PARAMETER ( LA_LINRX_ITREF_I = 1, 466 $ LA_LINRX_ITHRESH_I = 2 ) 467 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 468 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 469 $ LA_LINRX_RCOND_I 470 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 471 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 472* .. 473* .. External Subroutines .. 474 EXTERNAL CAXPY, CCOPY, CGBTRS, CGBMV, BLAS_CGBMV_X, 475 $ BLAS_CGBMV2_X, CLA_GBAMV, CLA_WWADDW, SLAMCH, 476 $ CHLA_TRANSTYPE, CLA_LIN_BERR 477 REAL SLAMCH 478 CHARACTER CHLA_TRANSTYPE 479* .. 480* .. Intrinsic Functions.. 481 INTRINSIC ABS, MAX, MIN 482* .. 483* .. Statement Functions .. 484 REAL CABS1 485* .. 486* .. Statement Function Definitions .. 487 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 488* .. 489* .. Executable Statements .. 490* 491 IF (INFO.NE.0) RETURN 492 TRANS = CHLA_TRANSTYPE(TRANS_TYPE) 493 EPS = SLAMCH( 'Epsilon' ) 494 HUGEVAL = SLAMCH( 'Overflow' ) 495* Force HUGEVAL to Inf 496 HUGEVAL = HUGEVAL * HUGEVAL 497* Using HUGEVAL may lead to spurious underflows. 498 INCR_THRESH = REAL( N ) * EPS 499 M = KL+KU+1 500 501 DO J = 1, NRHS 502 Y_PREC_STATE = EXTRA_RESIDUAL 503 IF ( Y_PREC_STATE .EQ. EXTRA_Y ) then 504 DO I = 1, N 505 Y_TAIL( I ) = 0.0 506 END DO 507 END IF 508 509 DXRAT = 0.0E+0 510 DXRATMAX = 0.0E+0 511 DZRAT = 0.0E+0 512 DZRATMAX = 0.0E+0 513 FINAL_DX_X = HUGEVAL 514 FINAL_DZ_Z = HUGEVAL 515 PREVNORMDX = HUGEVAL 516 PREV_DZ_Z = HUGEVAL 517 DZ_Z = HUGEVAL 518 DX_X = HUGEVAL 519 520 X_STATE = WORKING_STATE 521 Z_STATE = UNSTABLE_STATE 522 INCR_PREC = .FALSE. 523 524 DO CNT = 1, ITHRESH 525* 526* Compute residual RES = B_s - op(A_s) * Y, 527* op(A) = A, A**T, or A**H depending on TRANS (and type). 528* 529 CALL CCOPY( N, B( 1, J ), 1, RES, 1 ) 530 IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN 531 CALL CGBMV( TRANS, M, N, KL, KU, (-1.0E+0,0.0E+0), AB, 532 $ LDAB, Y( 1, J ), 1, (1.0E+0,0.0E+0), RES, 1 ) 533 ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN 534 CALL BLAS_CGBMV_X( TRANS_TYPE, N, N, KL, KU, 535 $ (-1.0E+0,0.0E+0), AB, LDAB, Y( 1, J ), 1, 536 $ (1.0E+0,0.0E+0), RES, 1, PREC_TYPE ) 537 ELSE 538 CALL BLAS_CGBMV2_X( TRANS_TYPE, N, N, KL, KU, 539 $ (-1.0E+0,0.0E+0), AB, LDAB, Y( 1, J ), Y_TAIL, 1, 540 $ (1.0E+0,0.0E+0), RES, 1, PREC_TYPE ) 541 END IF 542 543! XXX: RES is no longer needed. 544 CALL CCOPY( N, RES, 1, DY, 1 ) 545 CALL CGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, DY, N, 546 $ INFO ) 547* 548* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. 549* 550 NORMX = 0.0E+0 551 NORMY = 0.0E+0 552 NORMDX = 0.0E+0 553 DZ_Z = 0.0E+0 554 YMIN = HUGEVAL 555 556 DO I = 1, N 557 YK = CABS1( Y( I, J ) ) 558 DYK = CABS1( DY( I ) ) 559 560 IF (YK .NE. 0.0) THEN 561 DZ_Z = MAX( DZ_Z, DYK / YK ) 562 ELSE IF ( DYK .NE. 0.0 ) THEN 563 DZ_Z = HUGEVAL 564 END IF 565 566 YMIN = MIN( YMIN, YK ) 567 568 NORMY = MAX( NORMY, YK ) 569 570 IF ( COLEQU ) THEN 571 NORMX = MAX( NORMX, YK * C( I ) ) 572 NORMDX = MAX(NORMDX, DYK * C(I)) 573 ELSE 574 NORMX = NORMY 575 NORMDX = MAX( NORMDX, DYK ) 576 END IF 577 END DO 578 579 IF ( NORMX .NE. 0.0 ) THEN 580 DX_X = NORMDX / NORMX 581 ELSE IF ( NORMDX .EQ. 0.0 ) THEN 582 DX_X = 0.0 583 ELSE 584 DX_X = HUGEVAL 585 END IF 586 587 DXRAT = NORMDX / PREVNORMDX 588 DZRAT = DZ_Z / PREV_DZ_Z 589* 590* Check termination criteria. 591* 592 IF (.NOT.IGNORE_CWISE 593 $ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY 594 $ .AND. Y_PREC_STATE .LT. EXTRA_Y ) 595 $ INCR_PREC = .TRUE. 596 597 IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH ) 598 $ X_STATE = WORKING_STATE 599 IF ( X_STATE .EQ. WORKING_STATE ) THEN 600 IF ( DX_X .LE. EPS ) THEN 601 X_STATE = CONV_STATE 602 ELSE IF ( DXRAT .GT. RTHRESH ) THEN 603 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 604 INCR_PREC = .TRUE. 605 ELSE 606 X_STATE = NOPROG_STATE 607 END IF 608 ELSE 609 IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT 610 END IF 611 IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X 612 END IF 613 614 IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB ) 615 $ Z_STATE = WORKING_STATE 616 IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH ) 617 $ Z_STATE = WORKING_STATE 618 IF ( Z_STATE .EQ. WORKING_STATE ) THEN 619 IF ( DZ_Z .LE. EPS ) THEN 620 Z_STATE = CONV_STATE 621 ELSE IF ( DZ_Z .GT. DZ_UB ) THEN 622 Z_STATE = UNSTABLE_STATE 623 DZRATMAX = 0.0 624 FINAL_DZ_Z = HUGEVAL 625 ELSE IF ( DZRAT .GT. RTHRESH ) THEN 626 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 627 INCR_PREC = .TRUE. 628 ELSE 629 Z_STATE = NOPROG_STATE 630 END IF 631 ELSE 632 IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT 633 END IF 634 IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 635 END IF 636* 637* Exit if both normwise and componentwise stopped working, 638* but if componentwise is unstable, let it go at least two 639* iterations. 640* 641 IF ( X_STATE.NE.WORKING_STATE ) THEN 642 IF ( IGNORE_CWISE ) GOTO 666 643 IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE ) 644 $ GOTO 666 645 IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666 646 END IF 647 648 IF ( INCR_PREC ) THEN 649 INCR_PREC = .FALSE. 650 Y_PREC_STATE = Y_PREC_STATE + 1 651 DO I = 1, N 652 Y_TAIL( I ) = 0.0 653 END DO 654 END IF 655 656 PREVNORMDX = NORMDX 657 PREV_DZ_Z = DZ_Z 658* 659* Update soluton. 660* 661 IF ( Y_PREC_STATE .LT. EXTRA_Y ) THEN 662 CALL CAXPY( N, (1.0E+0,0.0E+0), DY, 1, Y(1,J), 1 ) 663 ELSE 664 CALL CLA_WWADDW( N, Y(1,J), Y_TAIL, DY ) 665 END IF 666 667 END DO 668* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT. 669 666 CONTINUE 670* 671* Set final_* when cnt hits ithresh. 672* 673 IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X 674 IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 675* 676* Compute error bounds. 677* 678 IF ( N_NORMS .GE. 1 ) THEN 679 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 680 $ FINAL_DX_X / (1 - DXRATMAX) 681 END IF 682 IF ( N_NORMS .GE. 2 ) THEN 683 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 684 $ FINAL_DZ_Z / (1 - DZRATMAX) 685 END IF 686* 687* Compute componentwise relative backward error from formula 688* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 689* where abs(Z) is the componentwise absolute value of the matrix 690* or vector Z. 691* 692* Compute residual RES = B_s - op(A_s) * Y, 693* op(A) = A, A**T, or A**H depending on TRANS (and type). 694* 695 CALL CCOPY( N, B( 1, J ), 1, RES, 1 ) 696 CALL CGBMV( TRANS, N, N, KL, KU, (-1.0E+0,0.0E+0), AB, LDAB, 697 $ Y(1,J), 1, (1.0E+0,0.0E+0), RES, 1 ) 698 699 DO I = 1, N 700 AYB( I ) = CABS1( B( I, J ) ) 701 END DO 702* 703* Compute abs(op(A_s))*abs(Y) + abs(B_s). 704* 705 CALL CLA_GBAMV( TRANS_TYPE, N, N, KL, KU, 1.0E+0, 706 $ AB, LDAB, Y(1, J), 1, 1.0E+0, AYB, 1 ) 707 708 CALL CLA_LIN_BERR( N, N, 1, RES, AYB, BERR_OUT( J ) ) 709* 710* End of loop for each RHS. 711* 712 END DO 713* 714 RETURN 715 END 716