LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zsyrfsx.f
Go to the documentation of this file.
00001 *> \brief \b ZSYRFSX
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZSYRFSX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsyrfsx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsyrfsx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsyrfsx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
00022 *                           S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
00023 *                           ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
00024 *                           WORK, RWORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       CHARACTER          UPLO, EQUED
00028 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
00029 *      $                   N_ERR_BNDS
00030 *       DOUBLE PRECISION   RCOND
00031 *       ..
00032 *       .. Array Arguments ..
00033 *       INTEGER            IPIV( * )
00034 *       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00035 *      $                   X( LDX, * ), WORK( * )
00036 *       DOUBLE PRECISION   S( * ), PARAMS( * ), BERR( * ), RWORK( * ),
00037 *      $                   ERR_BNDS_NORM( NRHS, * ),
00038 *      $                   ERR_BNDS_COMP( NRHS, * )
00039 *       ..
00040 *  
00041 *
00042 *> \par Purpose:
00043 *  =============
00044 *>
00045 *> \verbatim
00046 *>
00047 *>    ZSYRFSX improves the computed solution to a system of linear
00048 *>    equations when the coefficient matrix is symmetric indefinite, and
00049 *>    provides error bounds and backward error estimates for the
00050 *>    solution.  In addition to normwise error bound, the code provides
00051 *>    maximum componentwise error bound if possible.  See comments for
00052 *>    ERR_BNDS_NORM and ERR_BNDS_COMP for details of the error bounds.
00053 *>
00054 *>    The original system of linear equations may have been equilibrated
00055 *>    before calling this routine, as described by arguments EQUED and S
00056 *>    below. In this case, the solution and error bounds returned are
00057 *>    for the original unequilibrated system.
00058 *> \endverbatim
00059 *
00060 *  Arguments:
00061 *  ==========
00062 *
00063 *> \verbatim
00064 *>     Some optional parameters are bundled in the PARAMS array.  These
00065 *>     settings determine how refinement is performed, but often the
00066 *>     defaults are acceptable.  If the defaults are acceptable, users
00067 *>     can pass NPARAMS = 0 which prevents the source code from accessing
00068 *>     the PARAMS argument.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] UPLO
00072 *> \verbatim
00073 *>          UPLO is CHARACTER*1
00074 *>       = 'U':  Upper triangle of A is stored;
00075 *>       = 'L':  Lower triangle of A is stored.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] EQUED
00079 *> \verbatim
00080 *>          EQUED is CHARACTER*1
00081 *>     Specifies the form of equilibration that was done to A
00082 *>     before calling this routine. This is needed to compute
00083 *>     the solution and error bounds correctly.
00084 *>       = 'N':  No equilibration
00085 *>       = 'Y':  Both row and column equilibration, i.e., A has been
00086 *>               replaced by diag(S) * A * diag(S).
00087 *>               The right hand side B has been changed accordingly.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] N
00091 *> \verbatim
00092 *>          N is INTEGER
00093 *>     The order of the matrix A.  N >= 0.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] NRHS
00097 *> \verbatim
00098 *>          NRHS is INTEGER
00099 *>     The number of right hand sides, i.e., the number of columns
00100 *>     of the matrices B and X.  NRHS >= 0.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] A
00104 *> \verbatim
00105 *>          A is COMPLEX*16 array, dimension (LDA,N)
00106 *>     The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
00107 *>     upper triangular part of A contains the upper triangular
00108 *>     part of the matrix A, and the strictly lower triangular
00109 *>     part of A is not referenced.  If UPLO = 'L', the leading
00110 *>     N-by-N lower triangular part of A contains the lower
00111 *>     triangular part of the matrix A, and the strictly upper
00112 *>     triangular part of A is not referenced.
00113 *> \endverbatim
00114 *>
00115 *> \param[in] LDA
00116 *> \verbatim
00117 *>          LDA is INTEGER
00118 *>     The leading dimension of the array A.  LDA >= max(1,N).
00119 *> \endverbatim
00120 *>
00121 *> \param[in] AF
00122 *> \verbatim
00123 *>          AF is COMPLEX*16 array, dimension (LDAF,N)
00124 *>     The factored form of the matrix A.  AF contains the block
00125 *>     diagonal matrix D and the multipliers used to obtain the
00126 *>     factor U or L from the factorization A = U*D*U**T or A =
00127 *>     L*D*L**T as computed by DSYTRF.
00128 *> \endverbatim
00129 *>
00130 *> \param[in] LDAF
00131 *> \verbatim
00132 *>          LDAF is INTEGER
00133 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
00134 *> \endverbatim
00135 *>
00136 *> \param[in] IPIV
00137 *> \verbatim
00138 *>          IPIV is INTEGER array, dimension (N)
00139 *>     Details of the interchanges and the block structure of D
00140 *>     as determined by DSYTRF.
00141 *> \endverbatim
00142 *>
00143 *> \param[in,out] S
00144 *> \verbatim
00145 *>          S is DOUBLE PRECISION array, dimension (N)
00146 *>     The scale factors for A.  If EQUED = 'Y', A is multiplied on
00147 *>     the left and right by diag(S).  S is an input argument if FACT =
00148 *>     'F'; otherwise, S is an output argument.  If FACT = 'F' and EQUED
00149 *>     = 'Y', each element of S must be positive.  If S is output, each
00150 *>     element of S is a power of the radix. If S is input, each element
00151 *>     of S should be a power of the radix to ensure a reliable solution
00152 *>     and error estimates. Scaling by powers of the radix does not cause
00153 *>     rounding errors unless the result underflows or overflows.
00154 *>     Rounding errors during scaling lead to refining with a matrix that
00155 *>     is not equivalent to the input matrix, producing error estimates
00156 *>     that may not be reliable.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] B
00160 *> \verbatim
00161 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
00162 *>     The right hand side matrix B.
00163 *> \endverbatim
00164 *>
00165 *> \param[in] LDB
00166 *> \verbatim
00167 *>          LDB is INTEGER
00168 *>     The leading dimension of the array B.  LDB >= max(1,N).
00169 *> \endverbatim
00170 *>
00171 *> \param[in,out] X
00172 *> \verbatim
00173 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
00174 *>     On entry, the solution matrix X, as computed by DGETRS.
00175 *>     On exit, the improved solution matrix X.
00176 *> \endverbatim
00177 *>
00178 *> \param[in] LDX
00179 *> \verbatim
00180 *>          LDX is INTEGER
00181 *>     The leading dimension of the array X.  LDX >= max(1,N).
00182 *> \endverbatim
00183 *>
00184 *> \param[out] RCOND
00185 *> \verbatim
00186 *>          RCOND is DOUBLE PRECISION
00187 *>     Reciprocal scaled condition number.  This is an estimate of the
00188 *>     reciprocal Skeel condition number of the matrix A after
00189 *>     equilibration (if done).  If this is less than the machine
00190 *>     precision (in particular, if it is zero), the matrix is singular
00191 *>     to working precision.  Note that the error may still be small even
00192 *>     if this number is very small and the matrix appears ill-
00193 *>     conditioned.
00194 *> \endverbatim
00195 *>
00196 *> \param[out] BERR
00197 *> \verbatim
00198 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
00199 *>     Componentwise relative backward error.  This is the
00200 *>     componentwise relative backward error of each solution vector X(j)
00201 *>     (i.e., the smallest relative change in any element of A or B that
00202 *>     makes X(j) an exact solution).
00203 *> \endverbatim
00204 *>
00205 *> \param[in] N_ERR_BNDS
00206 *> \verbatim
00207 *>          N_ERR_BNDS is INTEGER
00208 *>     Number of error bounds to return for each right hand side
00209 *>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
00210 *>     ERR_BNDS_COMP below.
00211 *> \endverbatim
00212 *>
00213 *> \param[out] ERR_BNDS_NORM
00214 *> \verbatim
00215 *>          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
00216 *>     For each right-hand side, this array contains information about
00217 *>     various error bounds and condition numbers corresponding to the
00218 *>     normwise relative error, which is defined as follows:
00219 *>
00220 *>     Normwise relative error in the ith solution vector:
00221 *>             max_j (abs(XTRUE(j,i) - X(j,i)))
00222 *>            ------------------------------
00223 *>                  max_j abs(X(j,i))
00224 *>
00225 *>     The array is indexed by the type of error information as described
00226 *>     below. There currently are up to three pieces of information
00227 *>     returned.
00228 *>
00229 *>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
00230 *>     right-hand side.
00231 *>
00232 *>     The second index in ERR_BNDS_NORM(:,err) contains the following
00233 *>     three fields:
00234 *>     err = 1 "Trust/don't trust" boolean. Trust the answer if the
00235 *>              reciprocal condition number is less than the threshold
00236 *>              sqrt(n) * dlamch('Epsilon').
00237 *>
00238 *>     err = 2 "Guaranteed" error bound: The estimated forward error,
00239 *>              almost certainly within a factor of 10 of the true error
00240 *>              so long as the next entry is greater than the threshold
00241 *>              sqrt(n) * dlamch('Epsilon'). This error bound should only
00242 *>              be trusted if the previous boolean is true.
00243 *>
00244 *>     err = 3  Reciprocal condition number: Estimated normwise
00245 *>              reciprocal condition number.  Compared with the threshold
00246 *>              sqrt(n) * dlamch('Epsilon') to determine if the error
00247 *>              estimate is "guaranteed". These reciprocal condition
00248 *>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
00249 *>              appropriately scaled matrix Z.
00250 *>              Let Z = S*A, where S scales each row by a power of the
00251 *>              radix so all absolute row sums of Z are approximately 1.
00252 *>
00253 *>     See Lapack Working Note 165 for further details and extra
00254 *>     cautions.
00255 *> \endverbatim
00256 *>
00257 *> \param[out] ERR_BNDS_COMP
00258 *> \verbatim
00259 *>          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
00260 *>     For each right-hand side, this array contains information about
00261 *>     various error bounds and condition numbers corresponding to the
00262 *>     componentwise relative error, which is defined as follows:
00263 *>
00264 *>     Componentwise relative error in the ith solution vector:
00265 *>                    abs(XTRUE(j,i) - X(j,i))
00266 *>             max_j ----------------------
00267 *>                         abs(X(j,i))
00268 *>
00269 *>     The array is indexed by the right-hand side i (on which the
00270 *>     componentwise relative error depends), and the type of error
00271 *>     information as described below. There currently are up to three
00272 *>     pieces of information returned for each right-hand side. If
00273 *>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
00274 *>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
00275 *>     the first (:,N_ERR_BNDS) entries are returned.
00276 *>
00277 *>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
00278 *>     right-hand side.
00279 *>
00280 *>     The second index in ERR_BNDS_COMP(:,err) contains the following
00281 *>     three fields:
00282 *>     err = 1 "Trust/don't trust" boolean. Trust the answer if the
00283 *>              reciprocal condition number is less than the threshold
00284 *>              sqrt(n) * dlamch('Epsilon').
00285 *>
00286 *>     err = 2 "Guaranteed" error bound: The estimated forward error,
00287 *>              almost certainly within a factor of 10 of the true error
00288 *>              so long as the next entry is greater than the threshold
00289 *>              sqrt(n) * dlamch('Epsilon'). This error bound should only
00290 *>              be trusted if the previous boolean is true.
00291 *>
00292 *>     err = 3  Reciprocal condition number: Estimated componentwise
00293 *>              reciprocal condition number.  Compared with the threshold
00294 *>              sqrt(n) * dlamch('Epsilon') to determine if the error
00295 *>              estimate is "guaranteed". These reciprocal condition
00296 *>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
00297 *>              appropriately scaled matrix Z.
00298 *>              Let Z = S*(A*diag(x)), where x is the solution for the
00299 *>              current right-hand side and S scales each row of
00300 *>              A*diag(x) by a power of the radix so all absolute row
00301 *>              sums of Z are approximately 1.
00302 *>
00303 *>     See Lapack Working Note 165 for further details and extra
00304 *>     cautions.
00305 *> \endverbatim
00306 *>
00307 *> \param[in] NPARAMS
00308 *> \verbatim
00309 *>          NPARAMS is INTEGER
00310 *>     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
00311 *>     PARAMS array is never referenced and default values are used.
00312 *> \endverbatim
00313 *>
00314 *> \param[in,out] PARAMS
00315 *> \verbatim
00316 *>          PARAMS is / output) DOUBLE PRECISION array, dimension NPARAMS
00317 *>     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
00318 *>     that entry will be filled with default value used for that
00319 *>     parameter.  Only positions up to NPARAMS are accessed; defaults
00320 *>     are used for higher-numbered parameters.
00321 *>
00322 *>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
00323 *>            refinement or not.
00324 *>         Default: 1.0D+0
00325 *>            = 0.0 : No refinement is performed, and no error bounds are
00326 *>                    computed.
00327 *>            = 1.0 : Use the double-precision refinement algorithm,
00328 *>                    possibly with doubled-single computations if the
00329 *>                    compilation environment does not support DOUBLE
00330 *>                    PRECISION.
00331 *>              (other values are reserved for future use)
00332 *>
00333 *>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
00334 *>            computations allowed for refinement.
00335 *>         Default: 10
00336 *>         Aggressive: Set to 100 to permit convergence using approximate
00337 *>                     factorizations or factorizations other than LU. If
00338 *>                     the factorization uses a technique other than
00339 *>                     Gaussian elimination, the guarantees in
00340 *>                     err_bnds_norm and err_bnds_comp may no longer be
00341 *>                     trustworthy.
00342 *>
00343 *>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
00344 *>            will attempt to find a solution with small componentwise
00345 *>            relative error in the double-precision algorithm.  Positive
00346 *>            is true, 0.0 is false.
00347 *>         Default: 1.0 (attempt componentwise convergence)
00348 *> \endverbatim
00349 *>
00350 *> \param[out] WORK
00351 *> \verbatim
00352 *>          WORK is COMPLEX*16 array, dimension (2*N)
00353 *> \endverbatim
00354 *>
00355 *> \param[out] RWORK
00356 *> \verbatim
00357 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
00358 *> \endverbatim
00359 *>
00360 *> \param[out] INFO
00361 *> \verbatim
00362 *>          INFO is INTEGER
00363 *>       = 0:  Successful exit. The solution to every right-hand side is
00364 *>         guaranteed.
00365 *>       < 0:  If INFO = -i, the i-th argument had an illegal value
00366 *>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
00367 *>         has been completed, but the factor U is exactly singular, so
00368 *>         the solution and error bounds could not be computed. RCOND = 0
00369 *>         is returned.
00370 *>       = N+J: The solution corresponding to the Jth right-hand side is
00371 *>         not guaranteed. The solutions corresponding to other right-
00372 *>         hand sides K with K > J may not be guaranteed as well, but
00373 *>         only the first such right-hand side is reported. If a small
00374 *>         componentwise error is not requested (PARAMS(3) = 0.0) then
00375 *>         the Jth right-hand side is the first with a normwise error
00376 *>         bound that is not guaranteed (the smallest J such
00377 *>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
00378 *>         the Jth right-hand side is the first with either a normwise or
00379 *>         componentwise error bound that is not guaranteed (the smallest
00380 *>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
00381 *>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
00382 *>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
00383 *>         about all of the right-hand sides check ERR_BNDS_NORM or
00384 *>         ERR_BNDS_COMP.
00385 *> \endverbatim
00386 *
00387 *  Authors:
00388 *  ========
00389 *
00390 *> \author Univ. of Tennessee 
00391 *> \author Univ. of California Berkeley 
00392 *> \author Univ. of Colorado Denver 
00393 *> \author NAG Ltd. 
00394 *
00395 *> \date April 2012
00396 *
00397 *> \ingroup complex16SYcomputational
00398 *
00399 *  =====================================================================
00400       SUBROUTINE ZSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
00401      $                    S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
00402      $                    ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
00403      $                    WORK, RWORK, INFO )
00404 *
00405 *  -- LAPACK computational routine (version 3.4.1) --
00406 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00407 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00408 *     April 2012
00409 *
00410 *     .. Scalar Arguments ..
00411       CHARACTER          UPLO, EQUED
00412       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
00413      $                   N_ERR_BNDS
00414       DOUBLE PRECISION   RCOND
00415 *     ..
00416 *     .. Array Arguments ..
00417       INTEGER            IPIV( * )
00418       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00419      $                   X( LDX, * ), WORK( * )
00420       DOUBLE PRECISION   S( * ), PARAMS( * ), BERR( * ), RWORK( * ),
00421      $                   ERR_BNDS_NORM( NRHS, * ),
00422      $                   ERR_BNDS_COMP( NRHS, * )
00423 *     ..
00424 *
00425 *  ==================================================================
00426 *
00427 *     .. Parameters ..
00428       DOUBLE PRECISION   ZERO, ONE
00429       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00430       DOUBLE PRECISION   ITREF_DEFAULT, ITHRESH_DEFAULT
00431       DOUBLE PRECISION   COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
00432       DOUBLE PRECISION   DZTHRESH_DEFAULT
00433       PARAMETER          ( ITREF_DEFAULT = 1.0D+0 )
00434       PARAMETER          ( ITHRESH_DEFAULT = 10.0D+0 )
00435       PARAMETER          ( COMPONENTWISE_DEFAULT = 1.0D+0 )
00436       PARAMETER          ( RTHRESH_DEFAULT = 0.5D+0 )
00437       PARAMETER          ( DZTHRESH_DEFAULT = 0.25D+0 )
00438       INTEGER            LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
00439      $                   LA_LINRX_CWISE_I
00440       PARAMETER          ( LA_LINRX_ITREF_I = 1,
00441      $                   LA_LINRX_ITHRESH_I = 2 )
00442       PARAMETER          ( LA_LINRX_CWISE_I = 3 )
00443       INTEGER            LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
00444      $                   LA_LINRX_RCOND_I
00445       PARAMETER          ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
00446       PARAMETER          ( LA_LINRX_RCOND_I = 3 )
00447 *     ..
00448 *     .. Local Scalars ..
00449       CHARACTER(1)       NORM
00450       LOGICAL            RCEQU
00451       INTEGER            J, PREC_TYPE, REF_TYPE
00452       INTEGER            N_NORMS
00453       DOUBLE PRECISION   ANORM, RCOND_TMP
00454       DOUBLE PRECISION   ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
00455       LOGICAL            IGNORE_CWISE
00456       INTEGER            ITHRESH
00457       DOUBLE PRECISION   RTHRESH, UNSTABLE_THRESH
00458 *     ..
00459 *     .. External Subroutines ..
00460       EXTERNAL           XERBLA, ZSYCON, ZLA_SYRFSX_EXTENDED
00461 *     ..
00462 *     .. Intrinsic Functions ..
00463       INTRINSIC          MAX, SQRT, TRANSFER
00464 *     ..
00465 *     .. External Functions ..
00466       EXTERNAL           LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC
00467       EXTERNAL           DLAMCH, ZLANSY, ZLA_SYRCOND_X, ZLA_SYRCOND_C
00468       DOUBLE PRECISION   DLAMCH, ZLANSY, ZLA_SYRCOND_X, ZLA_SYRCOND_C
00469       LOGICAL            LSAME
00470       INTEGER            BLAS_FPINFO_X
00471       INTEGER            ILATRANS, ILAPREC
00472 *     ..
00473 *     .. Executable Statements ..
00474 *
00475 *     Check the input parameters.
00476 *
00477       INFO = 0
00478       REF_TYPE = INT( ITREF_DEFAULT )
00479       IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
00480          IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN
00481             PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
00482          ELSE
00483             REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
00484          END IF
00485       END IF
00486 *
00487 *     Set default parameters.
00488 *
00489       ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' )
00490       ITHRESH = INT( ITHRESH_DEFAULT )
00491       RTHRESH = RTHRESH_DEFAULT
00492       UNSTABLE_THRESH = DZTHRESH_DEFAULT
00493       IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0
00494 *
00495       IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
00496          IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN
00497             PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
00498          ELSE
00499             ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
00500          END IF
00501       END IF
00502       IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
00503          IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN
00504             IF ( IGNORE_CWISE ) THEN
00505                PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0
00506             ELSE
00507                PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0
00508             END IF
00509          ELSE
00510             IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0
00511          END IF
00512       END IF
00513       IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
00514          N_NORMS = 0
00515       ELSE IF ( IGNORE_CWISE ) THEN
00516          N_NORMS = 1
00517       ELSE
00518          N_NORMS = 2
00519       END IF
00520 *
00521       RCEQU = LSAME( EQUED, 'Y' )
00522 *
00523 *     Test input parameters.
00524 *
00525       IF ( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00526         INFO = -1
00527       ELSE IF( .NOT.RCEQU .AND. .NOT.LSAME( EQUED, 'N' ) ) THEN
00528         INFO = -2
00529       ELSE IF( N.LT.0 ) THEN
00530         INFO = -3
00531       ELSE IF( NRHS.LT.0 ) THEN
00532         INFO = -4
00533       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00534         INFO = -6
00535       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00536         INFO = -8
00537       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00538         INFO = -12
00539       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00540         INFO = -14
00541       END IF
00542       IF( INFO.NE.0 ) THEN
00543         CALL XERBLA( 'ZSYRFSX', -INFO )
00544         RETURN
00545       END IF
00546 *
00547 *     Quick return if possible.
00548 *
00549       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00550          RCOND = 1.0D+0
00551          DO J = 1, NRHS
00552             BERR( J ) = 0.0D+0
00553             IF ( N_ERR_BNDS .GE. 1 ) THEN
00554                ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
00555                ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
00556             END IF
00557             IF ( N_ERR_BNDS .GE. 2 ) THEN
00558                ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0
00559                ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0
00560             END IF
00561             IF ( N_ERR_BNDS .GE. 3 ) THEN
00562                ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0
00563                ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0
00564             END IF
00565          END DO
00566          RETURN
00567       END IF
00568 *
00569 *     Default to failure.
00570 *
00571       RCOND = 0.0D+0
00572       DO J = 1, NRHS
00573          BERR( J ) = 1.0D+0
00574          IF ( N_ERR_BNDS .GE. 1 ) THEN
00575             ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
00576             ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
00577          END IF
00578          IF ( N_ERR_BNDS .GE. 2 ) THEN
00579             ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
00580             ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
00581          END IF
00582          IF ( N_ERR_BNDS .GE. 3 ) THEN
00583             ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0
00584             ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0
00585          END IF
00586       END DO
00587 *
00588 *     Compute the norm of A and the reciprocal of the condition
00589 *     number of A.
00590 *
00591       NORM = 'I'
00592       ANORM = ZLANSY( NORM, UPLO, N, A, LDA, RWORK )
00593       CALL ZSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK,
00594      $     INFO )
00595 *
00596 *     Perform refinement on each right-hand side
00597 *
00598       IF ( REF_TYPE .NE. 0 ) THEN
00599 
00600          PREC_TYPE = ILAPREC( 'E' )
00601 
00602          CALL ZLA_SYRFSX_EXTENDED( PREC_TYPE, UPLO,  N,
00603      $        NRHS, A, LDA, AF, LDAF, IPIV, RCEQU, S, B,
00604      $        LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP,
00605      $        WORK, RWORK, WORK(N+1),
00606      $        TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), RCOND,
00607      $        ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
00608      $        INFO )
00609       END IF
00610 
00611       ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' )
00612       IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1) THEN
00613 *
00614 *     Compute scaled normwise condition number cond(A*C).
00615 *
00616          IF ( RCEQU ) THEN
00617             RCOND_TMP = ZLA_SYRCOND_C( UPLO, N, A, LDA, AF, LDAF, IPIV,
00618      $           S, .TRUE., INFO, WORK, RWORK )
00619          ELSE
00620             RCOND_TMP = ZLA_SYRCOND_C( UPLO, N, A, LDA, AF, LDAF, IPIV,
00621      $           S, .FALSE., INFO, WORK, RWORK )
00622          END IF
00623          DO J = 1, NRHS
00624 *
00625 *     Cap the error at 1.0.
00626 *
00627             IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
00628      $           .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
00629      $           ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
00630 *
00631 *     Threshold the error (see LAWN).
00632 *
00633             IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
00634                ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
00635                ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0
00636                IF ( INFO .LE. N ) INFO = N + J
00637             ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
00638      $              THEN
00639                ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
00640                ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
00641             END IF
00642 *
00643 *     Save the condition number.
00644 *
00645             IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
00646                ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
00647             END IF
00648          END DO
00649       END IF
00650 
00651       IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2 ) THEN
00652 *
00653 *     Compute componentwise condition number cond(A*diag(Y(:,J))) for
00654 *     each right-hand side using the current solution as an estimate of
00655 *     the true solution.  If the componentwise error estimate is too
00656 *     large, then the solution is a lousy estimate of truth and the
00657 *     estimated RCOND may be too optimistic.  To avoid misleading users,
00658 *     the inverse condition number is set to 0.0 when the estimated
00659 *     cwise error is at least CWISE_WRONG.
00660 *
00661          CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) )
00662          DO J = 1, NRHS
00663             IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
00664      $     THEN
00665                RCOND_TMP = ZLA_SYRCOND_X( UPLO, N, A, LDA, AF, LDAF,
00666      $         IPIV, X(1,J), INFO, WORK, RWORK )
00667             ELSE
00668                RCOND_TMP = 0.0D+0
00669             END IF
00670 *
00671 *     Cap the error at 1.0.
00672 *
00673             IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
00674      $           .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
00675      $           ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
00676 
00677 *
00678 *     Threshold the error (see LAWN).
00679 *
00680             IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
00681                ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
00682                ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0
00683                IF (.NOT. IGNORE_CWISE
00684      $              .AND. INFO.LT.N + J ) INFO = N + J
00685             ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
00686      $              .LT. ERR_LBND ) THEN
00687                ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
00688                ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
00689             END IF
00690 *
00691 *     Save the condition number.
00692 *
00693             IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
00694                ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
00695             END IF
00696 
00697          END DO
00698       END IF
00699 *
00700       RETURN
00701 *
00702 *     End of ZSYRFSX
00703 *
00704       END
 All Files Functions