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