SRC/zporfsx.f

   1 *> \brief \b ZPORFSX
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZPORFSX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zporfsx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zporfsx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zporfsx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B,
  22 *                           LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
  23 *                           ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
  24 *                           WORK, RWORK, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       CHARACTER          UPLO, EQUED
  28 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
  29 *      $                   N_ERR_BNDS
  30 *       DOUBLE PRECISION   RCOND
  31 *       ..
  32 *       .. Array Arguments ..
  33 *       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  34 *      $                   X( LDX, * ), WORK( * )
  35 *       DOUBLE PRECISION   RWORK( * ), S( * ), PARAMS(*), BERR( * ),
  36 *      $                   ERR_BNDS_NORM( NRHS, * ),
  37 *      $                   ERR_BNDS_COMP( NRHS, * )
  38 *       ..
  39 *
  40 *
  41 *> \par Purpose:
  42 *  =============
  43 *>
  44 *> \verbatim
  45 *>
  46 *>    ZPORFSX improves the computed solution to a system of linear
  47 *>    equations when the coefficient matrix is symmetric positive
  48 *>    definite, and provides error bounds and backward error estimates
  49 *>    for the solution.  In addition to normwise error bound, the code
  50 *>    provides maximum componentwise error bound if possible.  See
  51 *>    comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
  52 *>    error bounds.
  53 *>
  54 *>    The original system of linear equations may have been equilibrated
  55 *>    before calling this routine, as described by arguments EQUED and S
  56 *>    below. In this case, the solution and error bounds returned are
  57 *>    for the original unequilibrated system.
  58 *> \endverbatim
  59 *
  60 *  Arguments:
  61 *  ==========
  62 *
  63 *> \verbatim
  64 *>     Some optional parameters are bundled in the PARAMS array.  These
  65 *>     settings determine how refinement is performed, but often the
  66 *>     defaults are acceptable.  If the defaults are acceptable, users
  67 *>     can pass NPARAMS = 0 which prevents the source code from accessing
  68 *>     the PARAMS argument.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] UPLO
  72 *> \verbatim
  73 *>          UPLO is CHARACTER*1
  74 *>       = 'U':  Upper triangle of A is stored;
  75 *>       = 'L':  Lower triangle of A is stored.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] EQUED
  79 *> \verbatim
  80 *>          EQUED is CHARACTER*1
  81 *>     Specifies the form of equilibration that was done to A
  82 *>     before calling this routine. This is needed to compute
  83 *>     the solution and error bounds correctly.
  84 *>       = 'N':  No equilibration
  85 *>       = 'Y':  Both row and column equilibration, i.e., A has been
  86 *>               replaced by diag(S) * A * diag(S).
  87 *>               The right hand side B has been changed accordingly.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] N
  91 *> \verbatim
  92 *>          N is INTEGER
  93 *>     The order of the matrix A.  N >= 0.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] NRHS
  97 *> \verbatim
  98 *>          NRHS is INTEGER
  99 *>     The number of right hand sides, i.e., the number of columns
 100 *>     of the matrices B and X.  NRHS >= 0.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] A
 104 *> \verbatim
 105 *>          A is COMPLEX*16 array, dimension (LDA,N)
 106 *>     The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
 107 *>     upper triangular part of A contains the upper triangular part
 108 *>     of the matrix A, and the strictly lower triangular part of A
 109 *>     is not referenced.  If UPLO = 'L', the leading N-by-N lower
 110 *>     triangular part of A contains the lower triangular part of
 111 *>     the matrix A, and the strictly upper triangular part of A is
 112 *>     not referenced.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] LDA
 116 *> \verbatim
 117 *>          LDA is INTEGER
 118 *>     The leading dimension of the array A.  LDA >= max(1,N).
 119 *> \endverbatim
 120 *>
 121 *> \param[in] AF
 122 *> \verbatim
 123 *>          AF is COMPLEX*16 array, dimension (LDAF,N)
 124 *>     The triangular factor U or L from the Cholesky factorization
 125 *>     A = U**T*U or A = L*L**T, as computed by DPOTRF.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] LDAF
 129 *> \verbatim
 130 *>          LDAF is INTEGER
 131 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
 132 *> \endverbatim
 133 *>
 134 *> \param[in,out] S
 135 *> \verbatim
 136 *>          S is DOUBLE PRECISION array, dimension (N)
 137 *>     The row scale factors for A.  If EQUED = 'Y', A is multiplied on
 138 *>     the left and right by diag(S).  S is an input argument if FACT =
 139 *>     'F'; otherwise, S is an output argument.  If FACT = 'F' and EQUED
 140 *>     = 'Y', each element of S must be positive.  If S is output, each
 141 *>     element of S is a power of the radix. If S is input, each element
 142 *>     of S should be a power of the radix to ensure a reliable solution
 143 *>     and error estimates. Scaling by powers of the radix does not cause
 144 *>     rounding errors unless the result underflows or overflows.
 145 *>     Rounding errors during scaling lead to refining with a matrix that
 146 *>     is not equivalent to the input matrix, producing error estimates
 147 *>     that may not be reliable.
 148 *> \endverbatim
 149 *>
 150 *> \param[in] B
 151 *> \verbatim
 152 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 153 *>     The right hand side matrix B.
 154 *> \endverbatim
 155 *>
 156 *> \param[in] LDB
 157 *> \verbatim
 158 *>          LDB is INTEGER
 159 *>     The leading dimension of the array B.  LDB >= max(1,N).
 160 *> \endverbatim
 161 *>
 162 *> \param[in,out] X
 163 *> \verbatim
 164 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
 165 *>     On entry, the solution matrix X, as computed by DGETRS.
 166 *>     On exit, the improved solution matrix X.
 167 *> \endverbatim
 168 *>
 169 *> \param[in] LDX
 170 *> \verbatim
 171 *>          LDX is INTEGER
 172 *>     The leading dimension of the array X.  LDX >= max(1,N).
 173 *> \endverbatim
 174 *>
 175 *> \param[out] RCOND
 176 *> \verbatim
 177 *>          RCOND is DOUBLE PRECISION
 178 *>     Reciprocal scaled condition number.  This is an estimate of the
 179 *>     reciprocal Skeel condition number of the matrix A after
 180 *>     equilibration (if done).  If this is less than the machine
 181 *>     precision (in particular, if it is zero), the matrix is singular
 182 *>     to working precision.  Note that the error may still be small even
 183 *>     if this number is very small and the matrix appears ill-
 184 *>     conditioned.
 185 *> \endverbatim
 186 *>
 187 *> \param[out] BERR
 188 *> \verbatim
 189 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 190 *>     Componentwise relative backward error.  This is the
 191 *>     componentwise relative backward error of each solution vector X(j)
 192 *>     (i.e., the smallest relative change in any element of A or B that
 193 *>     makes X(j) an exact solution).
 194 *> \endverbatim
 195 *>
 196 *> \param[in] N_ERR_BNDS
 197 *> \verbatim
 198 *>          N_ERR_BNDS is INTEGER
 199 *>     Number of error bounds to return for each right hand side
 200 *>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
 201 *>     ERR_BNDS_COMP below.
 202 *> \endverbatim
 203 *>
 204 *> \param[out] ERR_BNDS_NORM
 205 *> \verbatim
 206 *>          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
 207 *>     For each right-hand side, this array contains information about
 208 *>     various error bounds and condition numbers corresponding to the
 209 *>     normwise relative error, which is defined as follows:
 210 *>
 211 *>     Normwise relative error in the ith solution vector:
 212 *>             max_j (abs(XTRUE(j,i) - X(j,i)))
 213 *>            ------------------------------
 214 *>                  max_j abs(X(j,i))
 215 *>
 216 *>     The array is indexed by the type of error information as described
 217 *>     below. There currently are up to three pieces of information
 218 *>     returned.
 219 *>
 220 *>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
 221 *>     right-hand side.
 222 *>
 223 *>     The second index in ERR_BNDS_NORM(:,err) contains the following
 224 *>     three fields:
 225 *>     err = 1 "Trust/don't trust" boolean. Trust the answer if the
 226 *>              reciprocal condition number is less than the threshold
 227 *>              sqrt(n) * dlamch('Epsilon').
 228 *>
 229 *>     err = 2 "Guaranteed" error bound: The estimated forward error,
 230 *>              almost certainly within a factor of 10 of the true error
 231 *>              so long as the next entry is greater than the threshold
 232 *>              sqrt(n) * dlamch('Epsilon'). This error bound should only
 233 *>              be trusted if the previous boolean is true.
 234 *>
 235 *>     err = 3  Reciprocal condition number: Estimated normwise
 236 *>              reciprocal condition number.  Compared with the threshold
 237 *>              sqrt(n) * dlamch('Epsilon') to determine if the error
 238 *>              estimate is "guaranteed". These reciprocal condition
 239 *>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
 240 *>              appropriately scaled matrix Z.
 241 *>              Let Z = S*A, where S scales each row by a power of the
 242 *>              radix so all absolute row sums of Z are approximately 1.
 243 *>
 244 *>     See Lapack Working Note 165 for further details and extra
 245 *>     cautions.
 246 *> \endverbatim
 247 *>
 248 *> \param[out] ERR_BNDS_COMP
 249 *> \verbatim
 250 *>          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
 251 *>     For each right-hand side, this array contains information about
 252 *>     various error bounds and condition numbers corresponding to the
 253 *>     componentwise relative error, which is defined as follows:
 254 *>
 255 *>     Componentwise relative error in the ith solution vector:
 256 *>                    abs(XTRUE(j,i) - X(j,i))
 257 *>             max_j ----------------------
 258 *>                         abs(X(j,i))
 259 *>
 260 *>     The array is indexed by the right-hand side i (on which the
 261 *>     componentwise relative error depends), and the type of error
 262 *>     information as described below. There currently are up to three
 263 *>     pieces of information returned for each right-hand side. If
 264 *>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
 265 *>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
 266 *>     the first (:,N_ERR_BNDS) entries are returned.
 267 *>
 268 *>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
 269 *>     right-hand side.
 270 *>
 271 *>     The second index in ERR_BNDS_COMP(:,err) contains the following
 272 *>     three fields:
 273 *>     err = 1 "Trust/don't trust" boolean. Trust the answer if the
 274 *>              reciprocal condition number is less than the threshold
 275 *>              sqrt(n) * dlamch('Epsilon').
 276 *>
 277 *>     err = 2 "Guaranteed" error bound: The estimated forward error,
 278 *>              almost certainly within a factor of 10 of the true error
 279 *>              so long as the next entry is greater than the threshold
 280 *>              sqrt(n) * dlamch('Epsilon'). This error bound should only
 281 *>              be trusted if the previous boolean is true.
 282 *>
 283 *>     err = 3  Reciprocal condition number: Estimated componentwise
 284 *>              reciprocal condition number.  Compared with the threshold
 285 *>              sqrt(n) * dlamch('Epsilon') to determine if the error
 286 *>              estimate is "guaranteed". These reciprocal condition
 287 *>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
 288 *>              appropriately scaled matrix Z.
 289 *>              Let Z = S*(A*diag(x)), where x is the solution for the
 290 *>              current right-hand side and S scales each row of
 291 *>              A*diag(x) by a power of the radix so all absolute row
 292 *>              sums of Z are approximately 1.
 293 *>
 294 *>     See Lapack Working Note 165 for further details and extra
 295 *>     cautions.
 296 *> \endverbatim
 297 *>
 298 *> \param[in] NPARAMS
 299 *> \verbatim
 300 *>          NPARAMS is INTEGER
 301 *>     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
 302 *>     PARAMS array is never referenced and default values are used.
 303 *> \endverbatim
 304 *>
 305 *> \param[in,out] PARAMS
 306 *> \verbatim
 307 *>          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
 308 *>     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
 309 *>     that entry will be filled with default value used for that
 310 *>     parameter.  Only positions up to NPARAMS are accessed; defaults
 311 *>     are used for higher-numbered parameters.
 312 *>
 313 *>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
 314 *>            refinement or not.
 315 *>         Default: 1.0D+0
 316 *>            = 0.0 : No refinement is performed, and no error bounds are
 317 *>                    computed.
 318 *>            = 1.0 : Use the double-precision refinement algorithm,
 319 *>                    possibly with doubled-single computations if the
 320 *>                    compilation environment does not support DOUBLE
 321 *>                    PRECISION.
 322 *>              (other values are reserved for future use)
 323 *>
 324 *>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
 325 *>            computations allowed for refinement.
 326 *>         Default: 10
 327 *>         Aggressive: Set to 100 to permit convergence using approximate
 328 *>                     factorizations or factorizations other than LU. If
 329 *>                     the factorization uses a technique other than
 330 *>                     Gaussian elimination, the guarantees in
 331 *>                     err_bnds_norm and err_bnds_comp may no longer be
 332 *>                     trustworthy.
 333 *>
 334 *>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
 335 *>            will attempt to find a solution with small componentwise
 336 *>            relative error in the double-precision algorithm.  Positive
 337 *>            is true, 0.0 is false.
 338 *>         Default: 1.0 (attempt componentwise convergence)
 339 *> \endverbatim
 340 *>
 341 *> \param[out] WORK
 342 *> \verbatim
 343 *>          WORK is COMPLEX*16 array, dimension (2*N)
 344 *> \endverbatim
 345 *>
 346 *> \param[out] RWORK
 347 *> \verbatim
 348 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
 349 *> \endverbatim
 350 *>
 351 *> \param[out] INFO
 352 *> \verbatim
 353 *>          INFO is INTEGER
 354 *>       = 0:  Successful exit. The solution to every right-hand side is
 355 *>         guaranteed.
 356 *>       < 0:  If INFO = -i, the i-th argument had an illegal value
 357 *>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
 358 *>         has been completed, but the factor U is exactly singular, so
 359 *>         the solution and error bounds could not be computed. RCOND = 0
 360 *>         is returned.
 361 *>       = N+J: The solution corresponding to the Jth right-hand side is
 362 *>         not guaranteed. The solutions corresponding to other right-
 363 *>         hand sides K with K > J may not be guaranteed as well, but
 364 *>         only the first such right-hand side is reported. If a small
 365 *>         componentwise error is not requested (PARAMS(3) = 0.0) then
 366 *>         the Jth right-hand side is the first with a normwise error
 367 *>         bound that is not guaranteed (the smallest J such
 368 *>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
 369 *>         the Jth right-hand side is the first with either a normwise or
 370 *>         componentwise error bound that is not guaranteed (the smallest
 371 *>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
 372 *>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
 373 *>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
 374 *>         about all of the right-hand sides check ERR_BNDS_NORM or
 375 *>         ERR_BNDS_COMP.
 376 *> \endverbatim
 377 *
 378 *  Authors:
 379 *  ========
 380 *
 381 *> \author Univ. of Tennessee
 382 *> \author Univ. of California Berkeley
 383 *> \author Univ. of Colorado Denver
 384 *> \author NAG Ltd.
 385 *
 386 *> \date April 2012
 387 *
 388 *> \ingroup complex16POcomputational
 389 *
 390 *  =====================================================================
 391       SUBROUTINE ZPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B,
 392      $                    LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
 393      $                    ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
 394      $                    WORK, RWORK, INFO )
 395 *
 396 *  -- LAPACK computational routine (version 3.4.1) --
 397 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 398 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 399 *     April 2012
 400 *
 401 *     .. Scalar Arguments ..
 402       CHARACTER          UPLO, EQUED
 403       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
 404      $                   N_ERR_BNDS
 405       DOUBLE PRECISION   RCOND
 406 *     ..
 407 *     .. Array Arguments ..
 408       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
 409      $                   X( LDX, * ), WORK( * )
 410       DOUBLE PRECISION   RWORK( * ), S( * ), PARAMS(*), BERR( * ),
 411      $                   ERR_BNDS_NORM( NRHS, * ),
 412      $                   ERR_BNDS_COMP( NRHS, * )
 413 *     ..
 414 *
 415 *  ==================================================================
 416 *
 417 *     .. Parameters ..
 418       DOUBLE PRECISION   ZERO, ONE
 419       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 420       DOUBLE PRECISION   ITREF_DEFAULT, ITHRESH_DEFAULT
 421       DOUBLE PRECISION   COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
 422       DOUBLE PRECISION   DZTHRESH_DEFAULT
 423       PARAMETER          ( ITREF_DEFAULT = 1.0D+0 )
 424       PARAMETER          ( ITHRESH_DEFAULT = 10.0D+0 )
 425       PARAMETER          ( COMPONENTWISE_DEFAULT = 1.0D+0 )
 426       PARAMETER          ( RTHRESH_DEFAULT = 0.5D+0 )
 427       PARAMETER          ( DZTHRESH_DEFAULT = 0.25D+0 )
 428       INTEGER            LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
 429      $                   LA_LINRX_CWISE_I
 430       PARAMETER          ( LA_LINRX_ITREF_I = 1,
 431      $                   LA_LINRX_ITHRESH_I = 2 )
 432       PARAMETER          ( LA_LINRX_CWISE_I = 3 )
 433       INTEGER            LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
 434      $                   LA_LINRX_RCOND_I
 435       PARAMETER          ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
 436       PARAMETER          ( LA_LINRX_RCOND_I = 3 )
 437 *     ..
 438 *     .. Local Scalars ..
 439       CHARACTER(1)       NORM
 440       LOGICAL            RCEQU
 441       INTEGER            J, PREC_TYPE, REF_TYPE
 442       INTEGER            N_NORMS
 443       DOUBLE PRECISION   ANORM, RCOND_TMP
 444       DOUBLE PRECISION   ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
 445       LOGICAL            IGNORE_CWISE
 446       INTEGER            ITHRESH
 447       DOUBLE PRECISION   RTHRESH, UNSTABLE_THRESH
 448 *     ..
 449 *     .. External Subroutines ..
 450       EXTERNAL           XERBLA, ZPOCON, ZLA_PORFSX_EXTENDED
 451 *     ..
 452 *     .. Intrinsic Functions ..
 453       INTRINSIC          MAX, SQRT, TRANSFER
 454 *     ..
 455 *     .. External Functions ..
 456       EXTERNAL           LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC
 457       EXTERNAL           DLAMCH, ZLANHE, ZLA_PORCOND_X, ZLA_PORCOND_C
 458       DOUBLE PRECISION   DLAMCH, ZLANHE, ZLA_PORCOND_X, ZLA_PORCOND_C
 459       LOGICAL            LSAME
 460       INTEGER            BLAS_FPINFO_X
 461       INTEGER            ILATRANS, ILAPREC
 462 *     ..
 463 *     .. Executable Statements ..
 464 *
 465 *     Check the input parameters.
 466 *
 467       INFO = 0
 468       REF_TYPE = INT( ITREF_DEFAULT )
 469       IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
 470          IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN
 471             PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
 472          ELSE
 473             REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
 474          END IF
 475       END IF
 476 *
 477 *     Set default parameters.
 478 *
 479       ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' )
 480       ITHRESH = INT( ITHRESH_DEFAULT )
 481       RTHRESH = RTHRESH_DEFAULT
 482       UNSTABLE_THRESH = DZTHRESH_DEFAULT
 483       IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0
 484 *
 485       IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
 486          IF ( PARAMS(LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN
 487             PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
 488          ELSE
 489             ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
 490          END IF
 491       END IF
 492       IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
 493          IF ( PARAMS(LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN
 494             IF ( IGNORE_CWISE ) THEN
 495                PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0
 496             ELSE
 497                PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0
 498             END IF
 499          ELSE
 500             IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0
 501          END IF
 502       END IF
 503       IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
 504          N_NORMS = 0
 505       ELSE IF ( IGNORE_CWISE ) THEN
 506          N_NORMS = 1
 507       ELSE
 508          N_NORMS = 2
 509       END IF
 510 *
 511       RCEQU = LSAME( EQUED, 'Y' )
 512 *
 513 *     Test input parameters.
 514 *
 515       IF (.NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 516         INFO = -1
 517       ELSE IF( .NOT.RCEQU .AND. .NOT.LSAME( EQUED, 'N' ) ) THEN
 518         INFO = -2
 519       ELSE IF( N.LT.0 ) THEN
 520         INFO = -3
 521       ELSE IF( NRHS.LT.0 ) THEN
 522         INFO = -4
 523       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 524         INFO = -6
 525       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 526         INFO = -8
 527       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 528         INFO = -11
 529       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 530         INFO = -13
 531       END IF
 532       IF( INFO.NE.0 ) THEN
 533         CALL XERBLA( 'ZPORFSX', -INFO )
 534         RETURN
 535       END IF
 536 *
 537 *     Quick return if possible.
 538 *
 539       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
 540          RCOND = 1.0D+0
 541          DO J = 1, NRHS
 542             BERR( J ) = 0.0D+0
 543             IF ( N_ERR_BNDS .GE. 1 ) THEN
 544                ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
 545                ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
 546             END IF
 547             IF ( N_ERR_BNDS .GE. 2 ) THEN
 548                ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0
 549                ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0
 550             END IF
 551             IF ( N_ERR_BNDS .GE. 3 ) THEN
 552                ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0
 553                ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0
 554             END IF
 555          END DO
 556          RETURN
 557       END IF
 558 *
 559 *     Default to failure.
 560 *
 561       RCOND = 0.0D+0
 562       DO J = 1, NRHS
 563          BERR( J ) = 1.0D+0
 564          IF ( N_ERR_BNDS .GE. 1 ) THEN
 565             ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
 566             ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
 567          END IF
 568          IF ( N_ERR_BNDS .GE. 2 ) THEN
 569             ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
 570             ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
 571          END IF
 572          IF ( N_ERR_BNDS .GE. 3 ) THEN
 573             ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0
 574             ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0
 575          END IF
 576       END DO
 577 *
 578 *     Compute the norm of A and the reciprocal of the condition
 579 *     number of A.
 580 *
 581       NORM = 'I'
 582       ANORM = ZLANHE( NORM, UPLO, N, A, LDA, RWORK )
 583       CALL ZPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK,
 584      $     INFO )
 585 *
 586 *     Perform refinement on each right-hand side
 587 *
 588       IF ( REF_TYPE .NE. 0 ) THEN
 589
 590          PREC_TYPE = ILAPREC( 'E' )
 591
 592          CALL ZLA_PORFSX_EXTENDED( PREC_TYPE, UPLO, N,
 593      $        NRHS, A, LDA, AF, LDAF, RCEQU, S, B,
 594      $        LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM, ERR_BNDS_COMP,
 595      $        WORK, RWORK, WORK(N+1),
 596      $        TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N), RCOND,
 597      $        ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
 598      $        INFO )
 599       END IF
 600
 601       ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' )
 602       IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN
 603 *
 604 *     Compute scaled normwise condition number cond(A*C).
 605 *
 606          IF ( RCEQU ) THEN
 607             RCOND_TMP = ZLA_PORCOND_C( UPLO, N, A, LDA, AF, LDAF,
 608      $           S, .TRUE., INFO, WORK, RWORK )
 609          ELSE
 610             RCOND_TMP = ZLA_PORCOND_C( UPLO, N, A, LDA, AF, LDAF,
 611      $           S, .FALSE., INFO, WORK, RWORK )
 612          END IF
 613          DO J = 1, NRHS
 614 *
 615 *     Cap the error at 1.0.
 616 *
 617             IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
 618      $           .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
 619      $           ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
 620 *
 621 *     Threshold the error (see LAWN).
 622 *
 623             IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
 624                ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
 625                ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0
 626                IF ( INFO .LE. N ) INFO = N + J
 627             ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
 628      $     THEN
 629                ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
 630                ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
 631             END IF
 632 *
 633 *     Save the condition number.
 634 *
 635             IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
 636                ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
 637             END IF
 638
 639          END DO
 640       END IF
 641
 642       IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN
 643 *
 644 *     Compute componentwise condition number cond(A*diag(Y(:,J))) for
 645 *     each right-hand side using the current solution as an estimate of
 646 *     the true solution.  If the componentwise error estimate is too
 647 *     large, then the solution is a lousy estimate of truth and the
 648 *     estimated RCOND may be too optimistic.  To avoid misleading users,
 649 *     the inverse condition number is set to 0.0 when the estimated
 650 *     cwise error is at least CWISE_WRONG.
 651 *
 652          CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) )
 653          DO J = 1, NRHS
 654             IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
 655      $     THEN
 656                RCOND_TMP = ZLA_PORCOND_X( UPLO, N, A, LDA, AF, LDAF,
 657      $         X(1,J), INFO, WORK, RWORK )
 658             ELSE
 659                RCOND_TMP = 0.0D+0
 660             END IF
 661 *
 662 *     Cap the error at 1.0.
 663 *
 664             IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
 665      $           .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
 666      $           ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
 667 *
 668 *     Threshold the error (see LAWN).
 669 *
 670             IF (RCOND_TMP .LT. ILLRCOND_THRESH) THEN
 671                ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
 672                ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0
 673                IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0
 674      $              .AND. INFO.LT.N + J ) INFO = N + J
 675             ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
 676      $              .LT. ERR_LBND ) THEN
 677                ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
 678                ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
 679             END IF
 680 *
 681 *     Save the condition number.
 682 *
 683             IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
 684                ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
 685             END IF
 686
 687          END DO
 688       END IF
 689 *
 690       RETURN
 691 *
 692 *     End of ZPORFSX
 693 *
 694       END