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