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