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