Annotation of rpl/lapack/lapack/dla_syrfsx_extended.f, revision 1.6

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

CVSweb interface <joel.bertrand@systella.fr>