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 definte 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*> \date June 2017 388* 389*> \ingroup complex16GEcomputational 390* 391* ===================================================================== 392 SUBROUTINE ZLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, NRHS, A, 393 $ LDA, AF, LDAF, IPIV, COLEQU, C, B, 394 $ LDB, Y, LDY, BERR_OUT, N_NORMS, 395 $ ERRS_N, ERRS_C, RES, AYB, DY, 396 $ Y_TAIL, RCOND, ITHRESH, RTHRESH, 397 $ DZ_UB, IGNORE_CWISE, INFO ) 398* 399* -- LAPACK computational routine (version 3.7.1) -- 400* -- LAPACK is a software package provided by Univ. of Tennessee, -- 401* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 402* June 2017 403* 404* .. Scalar Arguments .. 405 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE, 406 $ TRANS_TYPE, N_NORMS 407 LOGICAL COLEQU, IGNORE_CWISE 408 INTEGER ITHRESH 409 DOUBLE PRECISION RTHRESH, DZ_UB 410* .. 411* .. Array Arguments 412 INTEGER IPIV( * ) 413 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 414 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * ) 415 DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ), 416 $ ERRS_N( NRHS, * ), ERRS_C( NRHS, * ) 417* .. 418* 419* ===================================================================== 420* 421* .. Local Scalars .. 422 CHARACTER TRANS 423 INTEGER CNT, I, J, X_STATE, Z_STATE, Y_PREC_STATE 424 DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT, 425 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX, 426 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z, 427 $ EPS, HUGEVAL, INCR_THRESH 428 LOGICAL INCR_PREC 429 COMPLEX*16 ZDUM 430* .. 431* .. Parameters .. 432 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE, 433 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL, 434 $ EXTRA_Y 435 PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1, 436 $ CONV_STATE = 2, 437 $ NOPROG_STATE = 3 ) 438 PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1, 439 $ EXTRA_Y = 2 ) 440 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 441 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 442 INTEGER CMP_ERR_I, PIV_GROWTH_I 443 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 444 $ BERR_I = 3 ) 445 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 446 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8, 447 $ PIV_GROWTH_I = 9 ) 448 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 449 $ LA_LINRX_CWISE_I 450 PARAMETER ( LA_LINRX_ITREF_I = 1, 451 $ LA_LINRX_ITHRESH_I = 2 ) 452 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 453 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 454 $ LA_LINRX_RCOND_I 455 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 456 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 457* .. 458* .. External Subroutines .. 459 EXTERNAL ZAXPY, ZCOPY, ZGETRS, ZGEMV, BLAS_ZGEMV_X, 460 $ BLAS_ZGEMV2_X, ZLA_GEAMV, ZLA_WWADDW, DLAMCH, 461 $ CHLA_TRANSTYPE, ZLA_LIN_BERR 462 DOUBLE PRECISION DLAMCH 463 CHARACTER CHLA_TRANSTYPE 464* .. 465* .. Intrinsic Functions .. 466 INTRINSIC ABS, MAX, MIN 467* .. 468* .. Statement Functions .. 469 DOUBLE PRECISION CABS1 470* .. 471* .. Statement Function Definitions .. 472 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 473* .. 474* .. Executable Statements .. 475* 476 IF ( INFO.NE.0 ) RETURN 477 TRANS = CHLA_TRANSTYPE(TRANS_TYPE) 478 EPS = DLAMCH( 'Epsilon' ) 479 HUGEVAL = DLAMCH( 'Overflow' ) 480* Force HUGEVAL to Inf 481 HUGEVAL = HUGEVAL * HUGEVAL 482* Using HUGEVAL may lead to spurious underflows. 483 INCR_THRESH = DBLE( N ) * EPS 484* 485 DO J = 1, NRHS 486 Y_PREC_STATE = EXTRA_RESIDUAL 487 IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN 488 DO I = 1, N 489 Y_TAIL( I ) = 0.0D+0 490 END DO 491 END IF 492 493 DXRAT = 0.0D+0 494 DXRATMAX = 0.0D+0 495 DZRAT = 0.0D+0 496 DZRATMAX = 0.0D+0 497 FINAL_DX_X = HUGEVAL 498 FINAL_DZ_Z = HUGEVAL 499 PREVNORMDX = HUGEVAL 500 PREV_DZ_Z = HUGEVAL 501 DZ_Z = HUGEVAL 502 DX_X = HUGEVAL 503 504 X_STATE = WORKING_STATE 505 Z_STATE = UNSTABLE_STATE 506 INCR_PREC = .FALSE. 507 508 DO CNT = 1, ITHRESH 509* 510* Compute residual RES = B_s - op(A_s) * Y, 511* op(A) = A, A**T, or A**H depending on TRANS (and type). 512* 513 CALL ZCOPY( N, B( 1, J ), 1, RES, 1 ) 514 IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN 515 CALL ZGEMV( TRANS, N, N, (-1.0D+0,0.0D+0), A, LDA, 516 $ Y( 1, J ), 1, (1.0D+0,0.0D+0), RES, 1) 517 ELSE IF (Y_PREC_STATE .EQ. EXTRA_RESIDUAL) THEN 518 CALL BLAS_ZGEMV_X( TRANS_TYPE, N, N, (-1.0D+0,0.0D+0), A, 519 $ LDA, Y( 1, J ), 1, (1.0D+0,0.0D+0), 520 $ RES, 1, PREC_TYPE ) 521 ELSE 522 CALL BLAS_ZGEMV2_X( TRANS_TYPE, N, N, (-1.0D+0,0.0D+0), 523 $ A, LDA, Y(1, J), Y_TAIL, 1, (1.0D+0,0.0D+0), RES, 1, 524 $ PREC_TYPE) 525 END IF 526 527! XXX: RES is no longer needed. 528 CALL ZCOPY( N, RES, 1, DY, 1 ) 529 CALL ZGETRS( TRANS, N, 1, AF, LDAF, IPIV, DY, N, INFO ) 530* 531* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. 532* 533 NORMX = 0.0D+0 534 NORMY = 0.0D+0 535 NORMDX = 0.0D+0 536 DZ_Z = 0.0D+0 537 YMIN = HUGEVAL 538* 539 DO I = 1, N 540 YK = CABS1( Y( I, J ) ) 541 DYK = CABS1( DY( I ) ) 542 543 IF ( YK .NE. 0.0D+0 ) THEN 544 DZ_Z = MAX( DZ_Z, DYK / YK ) 545 ELSE IF ( DYK .NE. 0.0D+0 ) THEN 546 DZ_Z = HUGEVAL 547 END IF 548 549 YMIN = MIN( YMIN, YK ) 550 551 NORMY = MAX( NORMY, YK ) 552 553 IF ( COLEQU ) THEN 554 NORMX = MAX( NORMX, YK * C( I ) ) 555 NORMDX = MAX( NORMDX, DYK * C( I ) ) 556 ELSE 557 NORMX = NORMY 558 NORMDX = MAX(NORMDX, DYK) 559 END IF 560 END DO 561 562 IF ( NORMX .NE. 0.0D+0 ) THEN 563 DX_X = NORMDX / NORMX 564 ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN 565 DX_X = 0.0D+0 566 ELSE 567 DX_X = HUGEVAL 568 END IF 569 570 DXRAT = NORMDX / PREVNORMDX 571 DZRAT = DZ_Z / PREV_DZ_Z 572* 573* Check termination criteria 574* 575 IF (.NOT.IGNORE_CWISE 576 $ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY 577 $ .AND. Y_PREC_STATE .LT. EXTRA_Y ) 578 $ INCR_PREC = .TRUE. 579 580 IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH ) 581 $ X_STATE = WORKING_STATE 582 IF ( X_STATE .EQ. WORKING_STATE ) THEN 583 IF (DX_X .LE. EPS) THEN 584 X_STATE = CONV_STATE 585 ELSE IF ( DXRAT .GT. RTHRESH ) THEN 586 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 587 INCR_PREC = .TRUE. 588 ELSE 589 X_STATE = NOPROG_STATE 590 END IF 591 ELSE 592 IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT 593 END IF 594 IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X 595 END IF 596 597 IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB ) 598 $ Z_STATE = WORKING_STATE 599 IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH ) 600 $ Z_STATE = WORKING_STATE 601 IF ( Z_STATE .EQ. WORKING_STATE ) THEN 602 IF ( DZ_Z .LE. EPS ) THEN 603 Z_STATE = CONV_STATE 604 ELSE IF ( DZ_Z .GT. DZ_UB ) THEN 605 Z_STATE = UNSTABLE_STATE 606 DZRATMAX = 0.0D+0 607 FINAL_DZ_Z = HUGEVAL 608 ELSE IF ( DZRAT .GT. RTHRESH ) THEN 609 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 610 INCR_PREC = .TRUE. 611 ELSE 612 Z_STATE = NOPROG_STATE 613 END IF 614 ELSE 615 IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT 616 END IF 617 IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 618 END IF 619* 620* Exit if both normwise and componentwise stopped working, 621* but if componentwise is unstable, let it go at least two 622* iterations. 623* 624 IF ( X_STATE.NE.WORKING_STATE ) THEN 625 IF ( IGNORE_CWISE ) GOTO 666 626 IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE ) 627 $ GOTO 666 628 IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666 629 END IF 630 631 IF ( INCR_PREC ) THEN 632 INCR_PREC = .FALSE. 633 Y_PREC_STATE = Y_PREC_STATE + 1 634 DO I = 1, N 635 Y_TAIL( I ) = 0.0D+0 636 END DO 637 END IF 638 639 PREVNORMDX = NORMDX 640 PREV_DZ_Z = DZ_Z 641* 642* Update soluton. 643* 644 IF ( Y_PREC_STATE .LT. EXTRA_Y ) THEN 645 CALL ZAXPY( N, (1.0D+0,0.0D+0), DY, 1, Y(1,J), 1 ) 646 ELSE 647 CALL ZLA_WWADDW( N, Y( 1, J ), Y_TAIL, DY ) 648 END IF 649 650 END DO 651* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT. 652 666 CONTINUE 653* 654* Set final_* when cnt hits ithresh 655* 656 IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X 657 IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 658* 659* Compute error bounds 660* 661 IF (N_NORMS .GE. 1) THEN 662 ERRS_N( J, LA_LINRX_ERR_I ) = FINAL_DX_X / (1 - DXRATMAX) 663 664 END IF 665 IF ( N_NORMS .GE. 2 ) THEN 666 ERRS_C( J, LA_LINRX_ERR_I ) = FINAL_DZ_Z / (1 - DZRATMAX) 667 END IF 668* 669* Compute componentwise relative backward error from formula 670* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 671* where abs(Z) is the componentwise absolute value of the matrix 672* or vector Z. 673* 674* Compute residual RES = B_s - op(A_s) * Y, 675* op(A) = A, A**T, or A**H depending on TRANS (and type). 676* 677 CALL ZCOPY( N, B( 1, J ), 1, RES, 1 ) 678 CALL ZGEMV( TRANS, N, N, (-1.0D+0,0.0D+0), A, LDA, Y(1,J), 1, 679 $ (1.0D+0,0.0D+0), RES, 1 ) 680 681 DO I = 1, N 682 AYB( I ) = CABS1( B( I, J ) ) 683 END DO 684* 685* Compute abs(op(A_s))*abs(Y) + abs(B_s). 686* 687 CALL ZLA_GEAMV ( TRANS_TYPE, N, N, 1.0D+0, 688 $ A, LDA, Y(1, J), 1, 1.0D+0, AYB, 1 ) 689 690 CALL ZLA_LIN_BERR ( N, N, 1, RES, AYB, BERR_OUT( J ) ) 691* 692* End of loop for each RHS. 693* 694 END DO 695* 696 RETURN 697 END 698