Annotation of rpl/lapack/lapack/zla_porfsx_extended.f, revision 1.1

1.1     ! bertrand    1:       SUBROUTINE ZLA_PORFSX_EXTENDED( PREC_TYPE, UPLO, N, NRHS, A, LDA,
        !             2:      $                                AF, LDAF, COLEQU, C, B, LDB, Y,
        !             3:      $                                LDY, BERR_OUT, N_NORMS,
        !             4:      $                                ERR_BNDS_NORM, ERR_BNDS_COMP, RES,
        !             5:      $                                AYB, DY, Y_TAIL, RCOND, ITHRESH,
        !             6:      $                                RTHRESH, DZ_UB, IGNORE_CWISE,
        !             7:      $                                INFO )
        !             8: *
        !             9: *     -- LAPACK routine (version 3.2.2)                                 --
        !            10: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
        !            11: *     -- Jason Riedy of Univ. of California Berkeley.                 --
        !            12: *     -- June 2010                                                    --
        !            13: *
        !            14: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
        !            15: *     -- Univ. of California Berkeley and NAG Ltd.                    --
        !            16: *
        !            17:       IMPLICIT NONE
        !            18: *     ..
        !            19: *     .. Scalar Arguments ..
        !            20:       INTEGER            INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
        !            21:      $                   N_NORMS, ITHRESH
        !            22:       CHARACTER          UPLO
        !            23:       LOGICAL            COLEQU, IGNORE_CWISE
        !            24:       DOUBLE PRECISION   RTHRESH, DZ_UB
        !            25: *     ..
        !            26: *     .. Array Arguments ..
        !            27:       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
        !            28:      $                   Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
        !            29:       DOUBLE PRECISION   C( * ), AYB( * ), RCOND, BERR_OUT( * ),
        !            30:      $                   ERR_BNDS_NORM( NRHS, * ),
        !            31:      $                   ERR_BNDS_COMP( NRHS, * )
        !            32: *     ..
        !            33: *
        !            34: *  Purpose
        !            35: *  =======
        !            36: *
        !            37: *  ZLA_PORFSX_EXTENDED improves the computed solution to a system of
        !            38: *  linear equations by performing extra-precise iterative refinement
        !            39: *  and provides error bounds and backward error estimates for the solution.
        !            40: *  This subroutine is called by ZPORFSX to perform iterative refinement.
        !            41: *  In addition to normwise error bound, the code provides maximum
        !            42: *  componentwise error bound if possible. See comments for ERR_BNDS_NORM
        !            43: *  and ERR_BNDS_COMP for details of the error bounds. Note that this
        !            44: *  subroutine is only resonsible for setting the second fields of
        !            45: *  ERR_BNDS_NORM and ERR_BNDS_COMP.
        !            46: *
        !            47: *  Arguments
        !            48: *  =========
        !            49: *
        !            50: *     PREC_TYPE      (input) INTEGER
        !            51: *     Specifies the intermediate precision to be used in refinement.
        !            52: *     The value is defined by ILAPREC(P) where P is a CHARACTER and
        !            53: *     P    = 'S':  Single
        !            54: *          = 'D':  Double
        !            55: *          = 'I':  Indigenous
        !            56: *          = 'X', 'E':  Extra
        !            57: *
        !            58: *     UPLO    (input) CHARACTER*1
        !            59: *       = 'U':  Upper triangle of A is stored;
        !            60: *       = 'L':  Lower triangle of A is stored.
        !            61: *
        !            62: *     N              (input) INTEGER
        !            63: *     The number of linear equations, i.e., the order of the
        !            64: *     matrix A.  N >= 0.
        !            65: *
        !            66: *     NRHS           (input) INTEGER
        !            67: *     The number of right-hand-sides, i.e., the number of columns of the
        !            68: *     matrix B.
        !            69: *
        !            70: *     A              (input) COMPLEX*16 array, dimension (LDA,N)
        !            71: *     On entry, the N-by-N matrix A.
        !            72: *
        !            73: *     LDA            (input) INTEGER
        !            74: *     The leading dimension of the array A.  LDA >= max(1,N).
        !            75: *
        !            76: *     AF             (input) COMPLEX*16 array, dimension (LDAF,N)
        !            77: *     The triangular factor U or L from the Cholesky factorization
        !            78: *     A = U**T*U or A = L*L**T, as computed by ZPOTRF.
        !            79: *
        !            80: *     LDAF           (input) INTEGER
        !            81: *     The leading dimension of the array AF.  LDAF >= max(1,N).
        !            82: *
        !            83: *     COLEQU         (input) LOGICAL
        !            84: *     If .TRUE. then column equilibration was done to A before calling
        !            85: *     this routine. This is needed to compute the solution and error
        !            86: *     bounds correctly.
        !            87: *
        !            88: *     C              (input) DOUBLE PRECISION array, dimension (N)
        !            89: *     The column scale factors for A. If COLEQU = .FALSE., C
        !            90: *     is not accessed. If C is input, each element of C should be a power
        !            91: *     of the radix to ensure a reliable solution and error estimates.
        !            92: *     Scaling by powers of the radix does not cause rounding errors unless
        !            93: *     the result underflows or overflows. Rounding errors during scaling
        !            94: *     lead to refining with a matrix that is not equivalent to the
        !            95: *     input matrix, producing error estimates that may not be
        !            96: *     reliable.
        !            97: *
        !            98: *     B              (input) COMPLEX*16 array, dimension (LDB,NRHS)
        !            99: *     The right-hand-side matrix B.
        !           100: *
        !           101: *     LDB            (input) INTEGER
        !           102: *     The leading dimension of the array B.  LDB >= max(1,N).
        !           103: *
        !           104: *     Y              (input/output) COMPLEX*16 array, dimension
        !           105: *                    (LDY,NRHS)
        !           106: *     On entry, the solution matrix X, as computed by ZPOTRS.
        !           107: *     On exit, the improved solution matrix Y.
        !           108: *
        !           109: *     LDY            (input) INTEGER
        !           110: *     The leading dimension of the array Y.  LDY >= max(1,N).
        !           111: *
        !           112: *     BERR_OUT       (output) DOUBLE PRECISION array, dimension (NRHS)
        !           113: *     On exit, BERR_OUT(j) contains the componentwise relative backward
        !           114: *     error for right-hand-side j from the formula
        !           115: *         max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
        !           116: *     where abs(Z) is the componentwise absolute value of the matrix
        !           117: *     or vector Z. This is computed by ZLA_LIN_BERR.
        !           118: *
        !           119: *     N_NORMS        (input) INTEGER
        !           120: *     Determines which error bounds to return (see ERR_BNDS_NORM
        !           121: *     and ERR_BNDS_COMP).
        !           122: *     If N_NORMS >= 1 return normwise error bounds.
        !           123: *     If N_NORMS >= 2 return componentwise error bounds.
        !           124: *
        !           125: *     ERR_BNDS_NORM  (input/output) DOUBLE PRECISION array, dimension
        !           126: *                    (NRHS, N_ERR_BNDS)
        !           127: *     For each right-hand side, this array contains information about
        !           128: *     various error bounds and condition numbers corresponding to the
        !           129: *     normwise relative error, which is defined as follows:
        !           130: *
        !           131: *     Normwise relative error in the ith solution vector:
        !           132: *             max_j (abs(XTRUE(j,i) - X(j,i)))
        !           133: *            ------------------------------
        !           134: *                  max_j abs(X(j,i))
        !           135: *
        !           136: *     The array is indexed by the type of error information as described
        !           137: *     below. There currently are up to three pieces of information
        !           138: *     returned.
        !           139: *
        !           140: *     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
        !           141: *     right-hand side.
        !           142: *
        !           143: *     The second index in ERR_BNDS_NORM(:,err) contains the following
        !           144: *     three fields:
        !           145: *     err = 1 "Trust/don't trust" boolean. Trust the answer if the
        !           146: *              reciprocal condition number is less than the threshold
        !           147: *              sqrt(n) * slamch('Epsilon').
        !           148: *
        !           149: *     err = 2 "Guaranteed" error bound: The estimated forward error,
        !           150: *              almost certainly within a factor of 10 of the true error
        !           151: *              so long as the next entry is greater than the threshold
        !           152: *              sqrt(n) * slamch('Epsilon'). This error bound should only
        !           153: *              be trusted if the previous boolean is true.
        !           154: *
        !           155: *     err = 3  Reciprocal condition number: Estimated normwise
        !           156: *              reciprocal condition number.  Compared with the threshold
        !           157: *              sqrt(n) * slamch('Epsilon') to determine if the error
        !           158: *              estimate is "guaranteed". These reciprocal condition
        !           159: *              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
        !           160: *              appropriately scaled matrix Z.
        !           161: *              Let Z = S*A, where S scales each row by a power of the
        !           162: *              radix so all absolute row sums of Z are approximately 1.
        !           163: *
        !           164: *     This subroutine is only responsible for setting the second field
        !           165: *     above.
        !           166: *     See Lapack Working Note 165 for further details and extra
        !           167: *     cautions.
        !           168: *
        !           169: *     ERR_BNDS_COMP  (input/output) DOUBLE PRECISION array, dimension
        !           170: *                    (NRHS, N_ERR_BNDS)
        !           171: *     For each right-hand side, this array contains information about
        !           172: *     various error bounds and condition numbers corresponding to the
        !           173: *     componentwise relative error, which is defined as follows:
        !           174: *
        !           175: *     Componentwise relative error in the ith solution vector:
        !           176: *                    abs(XTRUE(j,i) - X(j,i))
        !           177: *             max_j ----------------------
        !           178: *                         abs(X(j,i))
        !           179: *
        !           180: *     The array is indexed by the right-hand side i (on which the
        !           181: *     componentwise relative error depends), and the type of error
        !           182: *     information as described below. There currently are up to three
        !           183: *     pieces of information returned for each right-hand side. If
        !           184: *     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
        !           185: *     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
        !           186: *     the first (:,N_ERR_BNDS) entries are returned.
        !           187: *
        !           188: *     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
        !           189: *     right-hand side.
        !           190: *
        !           191: *     The second index in ERR_BNDS_COMP(:,err) contains the following
        !           192: *     three fields:
        !           193: *     err = 1 "Trust/don't trust" boolean. Trust the answer if the
        !           194: *              reciprocal condition number is less than the threshold
        !           195: *              sqrt(n) * slamch('Epsilon').
        !           196: *
        !           197: *     err = 2 "Guaranteed" error bound: The estimated forward error,
        !           198: *              almost certainly within a factor of 10 of the true error
        !           199: *              so long as the next entry is greater than the threshold
        !           200: *              sqrt(n) * slamch('Epsilon'). This error bound should only
        !           201: *              be trusted if the previous boolean is true.
        !           202: *
        !           203: *     err = 3  Reciprocal condition number: Estimated componentwise
        !           204: *              reciprocal condition number.  Compared with the threshold
        !           205: *              sqrt(n) * slamch('Epsilon') to determine if the error
        !           206: *              estimate is "guaranteed". These reciprocal condition
        !           207: *              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
        !           208: *              appropriately scaled matrix Z.
        !           209: *              Let Z = S*(A*diag(x)), where x is the solution for the
        !           210: *              current right-hand side and S scales each row of
        !           211: *              A*diag(x) by a power of the radix so all absolute row
        !           212: *              sums of Z are approximately 1.
        !           213: *
        !           214: *     This subroutine is only responsible for setting the second field
        !           215: *     above.
        !           216: *     See Lapack Working Note 165 for further details and extra
        !           217: *     cautions.
        !           218: *
        !           219: *     RES            (input) COMPLEX*16 array, dimension (N)
        !           220: *     Workspace to hold the intermediate residual.
        !           221: *
        !           222: *     AYB            (input) DOUBLE PRECISION array, dimension (N)
        !           223: *     Workspace.
        !           224: *
        !           225: *     DY             (input) COMPLEX*16 PRECISION array, dimension (N)
        !           226: *     Workspace to hold the intermediate solution.
        !           227: *
        !           228: *     Y_TAIL         (input) COMPLEX*16 array, dimension (N)
        !           229: *     Workspace to hold the trailing bits of the intermediate solution.
        !           230: *
        !           231: *     RCOND          (input) DOUBLE PRECISION
        !           232: *     Reciprocal scaled condition number.  This is an estimate of the
        !           233: *     reciprocal Skeel condition number of the matrix A after
        !           234: *     equilibration (if done).  If this is less than the machine
        !           235: *     precision (in particular, if it is zero), the matrix is singular
        !           236: *     to working precision.  Note that the error may still be small even
        !           237: *     if this number is very small and the matrix appears ill-
        !           238: *     conditioned.
        !           239: *
        !           240: *     ITHRESH        (input) INTEGER
        !           241: *     The maximum number of residual computations allowed for
        !           242: *     refinement. The default is 10. For 'aggressive' set to 100 to
        !           243: *     permit convergence using approximate factorizations or
        !           244: *     factorizations other than LU. If the factorization uses a
        !           245: *     technique other than Gaussian elimination, the guarantees in
        !           246: *     ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
        !           247: *
        !           248: *     RTHRESH        (input) DOUBLE PRECISION
        !           249: *     Determines when to stop refinement if the error estimate stops
        !           250: *     decreasing. Refinement will stop when the next solution no longer
        !           251: *     satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
        !           252: *     the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
        !           253: *     default value is 0.5. For 'aggressive' set to 0.9 to permit
        !           254: *     convergence on extremely ill-conditioned matrices. See LAWN 165
        !           255: *     for more details.
        !           256: *
        !           257: *     DZ_UB          (input) DOUBLE PRECISION
        !           258: *     Determines when to start considering componentwise convergence.
        !           259: *     Componentwise convergence is only considered after each component
        !           260: *     of the solution Y is stable, which we definte as the relative
        !           261: *     change in each component being less than DZ_UB. The default value
        !           262: *     is 0.25, requiring the first bit to be stable. See LAWN 165 for
        !           263: *     more details.
        !           264: *
        !           265: *     IGNORE_CWISE   (input) LOGICAL
        !           266: *     If .TRUE. then ignore componentwise convergence. Default value
        !           267: *     is .FALSE..
        !           268: *
        !           269: *     INFO           (output) INTEGER
        !           270: *       = 0:  Successful exit.
        !           271: *       < 0:  if INFO = -i, the ith argument to ZPOTRS had an illegal
        !           272: *             value
        !           273: *
        !           274: *  =====================================================================
        !           275: *
        !           276: *     .. Local Scalars ..
        !           277:       INTEGER            UPLO2, CNT, I, J, X_STATE, Z_STATE,
        !           278:      $                   Y_PREC_STATE
        !           279:       DOUBLE PRECISION   YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
        !           280:      $                   DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
        !           281:      $                   DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
        !           282:      $                   EPS, HUGEVAL, INCR_THRESH
        !           283:       LOGICAL            INCR_PREC
        !           284:       COMPLEX*16         ZDUM
        !           285: *     ..
        !           286: *     .. Parameters ..
        !           287:       INTEGER            UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
        !           288:      $                   NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
        !           289:      $                   EXTRA_Y
        !           290:       PARAMETER          ( UNSTABLE_STATE = 0, WORKING_STATE = 1,
        !           291:      $                   CONV_STATE = 2, NOPROG_STATE = 3 )
        !           292:       PARAMETER          ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1,
        !           293:      $                   EXTRA_Y = 2 )
        !           294:       INTEGER            FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
        !           295:       INTEGER            RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
        !           296:       INTEGER            CMP_ERR_I, PIV_GROWTH_I
        !           297:       PARAMETER          ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
        !           298:      $                   BERR_I = 3 )
        !           299:       PARAMETER          ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
        !           300:       PARAMETER          ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
        !           301:      $                   PIV_GROWTH_I = 9 )
        !           302:       INTEGER            LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
        !           303:      $                   LA_LINRX_CWISE_I
        !           304:       PARAMETER          ( LA_LINRX_ITREF_I = 1,
        !           305:      $                   LA_LINRX_ITHRESH_I = 2 )
        !           306:       PARAMETER          ( LA_LINRX_CWISE_I = 3 )
        !           307:       INTEGER            LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
        !           308:      $                   LA_LINRX_RCOND_I
        !           309:       PARAMETER          ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
        !           310:       PARAMETER          ( LA_LINRX_RCOND_I = 3 )
        !           311: *     ..
        !           312: *     .. External Functions ..
        !           313:       LOGICAL            LSAME
        !           314:       EXTERNAL           ILAUPLO
        !           315:       INTEGER            ILAUPLO
        !           316: *     ..
        !           317: *     .. External Subroutines ..
        !           318:       EXTERNAL           ZAXPY, ZCOPY, ZPOTRS, ZHEMV, BLAS_ZHEMV_X,
        !           319:      $                   BLAS_ZHEMV2_X, ZLA_HEAMV, ZLA_WWADDW,
        !           320:      $                   ZLA_LIN_BERR, DLAMCH
        !           321:       DOUBLE PRECISION   DLAMCH
        !           322: *     ..
        !           323: *     .. Intrinsic Functions ..
        !           324:       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
        !           325: *     ..
        !           326: *     .. Statement Functions ..
        !           327:       DOUBLE PRECISION   CABS1
        !           328: *     ..
        !           329: *     .. Statement Function Definitions ..
        !           330:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
        !           331: *     ..
        !           332: *     .. Executable Statements ..
        !           333: *
        !           334:       IF (INFO.NE.0) RETURN
        !           335:       EPS = DLAMCH( 'Epsilon' )
        !           336:       HUGEVAL = DLAMCH( 'Overflow' )
        !           337: *     Force HUGEVAL to Inf
        !           338:       HUGEVAL = HUGEVAL * HUGEVAL
        !           339: *     Using HUGEVAL may lead to spurious underflows.
        !           340:       INCR_THRESH = DBLE(N) * EPS
        !           341: 
        !           342:       IF (LSAME (UPLO, 'L')) THEN
        !           343:          UPLO2 = ILAUPLO( 'L' )
        !           344:       ELSE
        !           345:          UPLO2 = ILAUPLO( 'U' )
        !           346:       ENDIF
        !           347: 
        !           348:       DO J = 1, NRHS
        !           349:          Y_PREC_STATE = EXTRA_RESIDUAL
        !           350:          IF (Y_PREC_STATE .EQ. EXTRA_Y) THEN
        !           351:             DO I = 1, N
        !           352:                Y_TAIL( I ) = 0.0D+0
        !           353:             END DO
        !           354:          END IF
        !           355: 
        !           356:          DXRAT = 0.0D+0
        !           357:          DXRATMAX = 0.0D+0
        !           358:          DZRAT = 0.0D+0
        !           359:          DZRATMAX = 0.0D+0
        !           360:          FINAL_DX_X = HUGEVAL
        !           361:          FINAL_DZ_Z = HUGEVAL
        !           362:          PREVNORMDX = HUGEVAL
        !           363:          PREV_DZ_Z = HUGEVAL
        !           364:          DZ_Z = HUGEVAL
        !           365:          DX_X = HUGEVAL
        !           366: 
        !           367:          X_STATE = WORKING_STATE
        !           368:          Z_STATE = UNSTABLE_STATE
        !           369:          INCR_PREC = .FALSE.
        !           370: 
        !           371:          DO CNT = 1, ITHRESH
        !           372: *
        !           373: *         Compute residual RES = B_s - op(A_s) * Y,
        !           374: *             op(A) = A, A**T, or A**H depending on TRANS (and type).
        !           375: *
        !           376:             CALL ZCOPY( N, B( 1, J ), 1, RES, 1 )
        !           377:             IF (Y_PREC_STATE .EQ. BASE_RESIDUAL) THEN
        !           378:                CALL ZHEMV(UPLO, N, DCMPLX(-1.0D+0), A, LDA, Y(1,J), 1,
        !           379:      $              DCMPLX(1.0D+0), RES, 1)
        !           380:             ELSE IF (Y_PREC_STATE .EQ. EXTRA_RESIDUAL) THEN
        !           381:                CALL BLAS_ZHEMV_X(UPLO2, N, DCMPLX(-1.0D+0), A, LDA,
        !           382:      $              Y( 1, J ), 1, DCMPLX(1.0D+0), RES, 1, PREC_TYPE)
        !           383:             ELSE
        !           384:                CALL BLAS_ZHEMV2_X(UPLO2, N, DCMPLX(-1.0D+0), A, LDA,
        !           385:      $              Y(1, J), Y_TAIL, 1, DCMPLX(1.0D+0), RES, 1,
        !           386:      $     PREC_TYPE)
        !           387:             END IF
        !           388: 
        !           389: !         XXX: RES is no longer needed.
        !           390:             CALL ZCOPY( N, RES, 1, DY, 1 )
        !           391:             CALL ZPOTRS( UPLO, N, 1, AF, LDAF, DY, N, INFO)
        !           392: *
        !           393: *         Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
        !           394: *
        !           395:             NORMX = 0.0D+0
        !           396:             NORMY = 0.0D+0
        !           397:             NORMDX = 0.0D+0
        !           398:             DZ_Z = 0.0D+0
        !           399:             YMIN = HUGEVAL
        !           400: 
        !           401:             DO I = 1, N
        !           402:                YK = CABS1(Y(I, J))
        !           403:                DYK = CABS1(DY(I))
        !           404: 
        !           405:                IF (YK .NE. 0.0D+0) THEN
        !           406:                   DZ_Z = MAX( DZ_Z, DYK / YK )
        !           407:                ELSE IF (DYK .NE. 0.0D+0) THEN
        !           408:                   DZ_Z = HUGEVAL
        !           409:                END IF
        !           410: 
        !           411:                YMIN = MIN( YMIN, YK )
        !           412: 
        !           413:                NORMY = MAX( NORMY, YK )
        !           414: 
        !           415:                IF ( COLEQU ) THEN
        !           416:                   NORMX = MAX(NORMX, YK * C(I))
        !           417:                   NORMDX = MAX(NORMDX, DYK * C(I))
        !           418:                ELSE
        !           419:                   NORMX = NORMY
        !           420:                   NORMDX = MAX(NORMDX, DYK)
        !           421:                END IF
        !           422:             END DO
        !           423: 
        !           424:             IF (NORMX .NE. 0.0D+0) THEN
        !           425:                DX_X = NORMDX / NORMX
        !           426:             ELSE IF (NORMDX .EQ. 0.0D+0) THEN
        !           427:                DX_X = 0.0D+0
        !           428:             ELSE
        !           429:                DX_X = HUGEVAL
        !           430:             END IF
        !           431: 
        !           432:             DXRAT = NORMDX / PREVNORMDX
        !           433:             DZRAT = DZ_Z / PREV_DZ_Z
        !           434: *
        !           435: *         Check termination criteria.
        !           436: *
        !           437:             IF (YMIN*RCOND .LT. INCR_THRESH*NORMY
        !           438:      $           .AND. Y_PREC_STATE .LT. EXTRA_Y)
        !           439:      $           INCR_PREC = .TRUE.
        !           440: 
        !           441:             IF (X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH)
        !           442:      $           X_STATE = WORKING_STATE
        !           443:             IF (X_STATE .EQ. WORKING_STATE) THEN
        !           444:                IF (DX_X .LE. EPS) THEN
        !           445:                   X_STATE = CONV_STATE
        !           446:                ELSE IF (DXRAT .GT. RTHRESH) THEN
        !           447:                   IF (Y_PREC_STATE .NE. EXTRA_Y) THEN
        !           448:                      INCR_PREC = .TRUE.
        !           449:                   ELSE
        !           450:                      X_STATE = NOPROG_STATE
        !           451:                   END IF
        !           452:                ELSE
        !           453:                   IF (DXRAT .GT. DXRATMAX) DXRATMAX = DXRAT
        !           454:                END IF
        !           455:                IF (X_STATE .GT. WORKING_STATE) FINAL_DX_X = DX_X
        !           456:             END IF
        !           457: 
        !           458:             IF (Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB)
        !           459:      $           Z_STATE = WORKING_STATE
        !           460:             IF (Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH)
        !           461:      $           Z_STATE = WORKING_STATE
        !           462:             IF (Z_STATE .EQ. WORKING_STATE) THEN
        !           463:                IF (DZ_Z .LE. EPS) THEN
        !           464:                   Z_STATE = CONV_STATE
        !           465:                ELSE IF (DZ_Z .GT. DZ_UB) THEN
        !           466:                   Z_STATE = UNSTABLE_STATE
        !           467:                   DZRATMAX = 0.0D+0
        !           468:                   FINAL_DZ_Z = HUGEVAL
        !           469:                ELSE IF (DZRAT .GT. RTHRESH) THEN
        !           470:                   IF (Y_PREC_STATE .NE. EXTRA_Y) THEN
        !           471:                      INCR_PREC = .TRUE.
        !           472:                   ELSE
        !           473:                      Z_STATE = NOPROG_STATE
        !           474:                   END IF
        !           475:                ELSE
        !           476:                   IF (DZRAT .GT. DZRATMAX) DZRATMAX = DZRAT
        !           477:                END IF
        !           478:                IF (Z_STATE .GT. WORKING_STATE) FINAL_DZ_Z = DZ_Z
        !           479:             END IF
        !           480: 
        !           481:             IF ( X_STATE.NE.WORKING_STATE.AND.
        !           482:      $           (IGNORE_CWISE.OR.Z_STATE.NE.WORKING_STATE) )
        !           483:      $           GOTO 666
        !           484: 
        !           485:             IF (INCR_PREC) THEN
        !           486:                INCR_PREC = .FALSE.
        !           487:                Y_PREC_STATE = Y_PREC_STATE + 1
        !           488:                DO I = 1, N
        !           489:                   Y_TAIL( I ) = 0.0D+0
        !           490:                END DO
        !           491:             END IF
        !           492: 
        !           493:             PREVNORMDX = NORMDX
        !           494:             PREV_DZ_Z = DZ_Z
        !           495: *
        !           496: *           Update soluton.
        !           497: *
        !           498:             IF (Y_PREC_STATE .LT. EXTRA_Y) THEN
        !           499:                CALL ZAXPY( N, DCMPLX(1.0D+0), DY, 1, Y(1,J), 1 )
        !           500:             ELSE
        !           501:                CALL ZLA_WWADDW(N, Y(1,J), Y_TAIL, DY)
        !           502:             END IF
        !           503: 
        !           504:          END DO
        !           505: *        Target of "IF (Z_STOP .AND. X_STOP)".  Sun's f77 won't EXIT.
        !           506:  666     CONTINUE
        !           507: *
        !           508: *     Set final_* when cnt hits ithresh.
        !           509: *
        !           510:          IF (X_STATE .EQ. WORKING_STATE) FINAL_DX_X = DX_X
        !           511:          IF (Z_STATE .EQ. WORKING_STATE) FINAL_DZ_Z = DZ_Z
        !           512: *
        !           513: *     Compute error bounds.
        !           514: *
        !           515:          IF (N_NORMS .GE. 1) THEN
        !           516:             ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) =
        !           517:      $           FINAL_DX_X / (1 - DXRATMAX)
        !           518:          END IF
        !           519:          IF (N_NORMS .GE. 2) THEN
        !           520:             ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) =
        !           521:      $           FINAL_DZ_Z / (1 - DZRATMAX)
        !           522:          END IF
        !           523: *
        !           524: *     Compute componentwise relative backward error from formula
        !           525: *         max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
        !           526: *     where abs(Z) is the componentwise absolute value of the matrix
        !           527: *     or vector Z.
        !           528: *
        !           529: *        Compute residual RES = B_s - op(A_s) * Y,
        !           530: *            op(A) = A, A**T, or A**H depending on TRANS (and type).
        !           531: *
        !           532:          CALL ZCOPY( N, B( 1, J ), 1, RES, 1 )
        !           533:          CALL ZHEMV(UPLO, N, DCMPLX(-1.0D+0), A, LDA, Y(1,J), 1,
        !           534:      $        DCMPLX(1.0D+0), RES, 1)
        !           535: 
        !           536:          DO I = 1, N
        !           537:             AYB( I ) = CABS1( B( I, J ) )
        !           538:          END DO
        !           539: *
        !           540: *     Compute abs(op(A_s))*abs(Y) + abs(B_s).
        !           541: *
        !           542:          CALL ZLA_HEAMV (UPLO2, N, 1.0D+0,
        !           543:      $        A, LDA, Y(1, J), 1, 1.0D+0, AYB, 1)
        !           544: 
        !           545:          CALL ZLA_LIN_BERR (N, N, 1, RES, AYB, BERR_OUT(J))
        !           546: *
        !           547: *     End of loop for each RHS.
        !           548: *
        !           549:       END DO
        !           550: *
        !           551:       RETURN
        !           552:       END

CVSweb interface <joel.bertrand@systella.fr>