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