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