LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
cgerfsx.f
Go to the documentation of this file.
1 *> \brief \b CGERFSX
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGERFSX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgerfsx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgerfsx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgerfsx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGERFSX( 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, RWORK, 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( * )
34 * COMPLEX 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, * ), RWORK( * )
39 * ..
40 *
41 *
42 *> \par Purpose:
43 * =============
44 *>
45 *> \verbatim
46 *>
47 *> CGERFSX 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 COMPLEX 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 COMPLEX array, dimension (LDAF,N)
124 *> The factors L and U from the factorization A = P*L*U
125 *> as computed by CGETRF.
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 CGETRF; 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 COMPLEX 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 COMPLEX array, dimension (LDX,NRHS)
186 *> On entry, the solution matrix X, as computed by CGETRS.
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 COMPLEX array, dimension (2*N)
365 *> \endverbatim
366 *>
367 *> \param[out] RWORK
368 *> \verbatim
369 *> RWORK is REAL array, dimension (2*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 December 2016
408 *
409 *> \ingroup complexGEcomputational
410 *
411 * =====================================================================
412  SUBROUTINE cgerfsx( 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, RWORK, INFO )
416 *
417 * -- LAPACK computational routine (version 3.7.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 * December 2016
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( * )
430  COMPLEX 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, * ), rwork( * )
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, transfer
476 * ..
477 * .. External Functions ..
478  EXTERNAL lsame, ilatrans, ilaprec
480  REAL SLAMCH, CLANGE, CLA_GERCOND_X, CLA_GERCOND_C
481  LOGICAL LSAME
482  INTEGER ILATRANS, ILAPREC
483 * ..
484 * .. Executable Statements ..
485 *
486 * Check the input parameters.
487 *
488  info = 0
489  trans_type = ilatrans( trans )
490  ref_type = int( itref_default )
491  IF ( nparams .GE. la_linrx_itref_i ) THEN
492  IF ( params( la_linrx_itref_i ) .LT. 0.0 ) THEN
493  params( la_linrx_itref_i ) = itref_default
494  ELSE
495  ref_type = params( la_linrx_itref_i )
496  END IF
497  END IF
498 *
499 * Set default parameters.
500 *
501  illrcond_thresh = REAL( N ) * SLAMCH( 'Epsilon' )
502  ithresh = int( ithresh_default )
503  rthresh = rthresh_default
504  unstable_thresh = dzthresh_default
505  ignore_cwise = componentwise_default .EQ. 0.0
506 *
507  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
508  IF ( params( la_linrx_ithresh_i ).LT.0.0 ) THEN
509  params(la_linrx_ithresh_i) = ithresh
510  ELSE
511  ithresh = int( params( la_linrx_ithresh_i ) )
512  END IF
513  END IF
514  IF ( nparams.GE.la_linrx_cwise_i ) THEN
515  IF ( params( la_linrx_cwise_i ).LT.0.0 ) THEN
516  IF ( ignore_cwise ) THEN
517  params( la_linrx_cwise_i ) = 0.0
518  ELSE
519  params( la_linrx_cwise_i ) = 1.0
520  END IF
521  ELSE
522  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0
523  END IF
524  END IF
525  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
526  n_norms = 0
527  ELSE IF ( ignore_cwise ) THEN
528  n_norms = 1
529  ELSE
530  n_norms = 2
531  END IF
532 *
533  notran = lsame( trans, 'N' )
534  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
535  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
536 *
537 * Test input parameters.
538 *
539  IF( trans_type.EQ.-1 ) THEN
540  info = -1
541  ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
542  $ .NOT.lsame( equed, 'N' ) ) THEN
543  info = -2
544  ELSE IF( n.LT.0 ) THEN
545  info = -3
546  ELSE IF( nrhs.LT.0 ) THEN
547  info = -4
548  ELSE IF( lda.LT.max( 1, n ) ) THEN
549  info = -6
550  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
551  info = -8
552  ELSE IF( ldb.LT.max( 1, n ) ) THEN
553  info = -13
554  ELSE IF( ldx.LT.max( 1, n ) ) THEN
555  info = -15
556  END IF
557  IF( info.NE.0 ) THEN
558  CALL xerbla( 'CGERFSX', -info )
559  RETURN
560  END IF
561 *
562 * Quick return if possible.
563 *
564  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
565  rcond = 1.0
566  DO j = 1, nrhs
567  berr( j ) = 0.0
568  IF ( n_err_bnds .GE. 1 ) THEN
569  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
570  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
571  END IF
572  IF ( n_err_bnds .GE. 2 ) THEN
573  err_bnds_norm( j, la_linrx_err_i ) = 0.0
574  err_bnds_comp( j, la_linrx_err_i ) = 0.0
575  END IF
576  IF ( n_err_bnds .GE. 3 ) THEN
577  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0
578  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
579  END IF
580  END DO
581  RETURN
582  END IF
583 *
584 * Default to failure.
585 *
586  rcond = 0.0
587  DO j = 1, nrhs
588  berr( j ) = 1.0
589  IF ( n_err_bnds .GE. 1 ) THEN
590  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
591  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
592  END IF
593  IF ( n_err_bnds .GE. 2 ) THEN
594  err_bnds_norm( j, la_linrx_err_i ) = 1.0
595  err_bnds_comp( j, la_linrx_err_i ) = 1.0
596  END IF
597  IF ( n_err_bnds .GE. 3 ) THEN
598  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
599  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0
600  END IF
601  END DO
602 *
603 * Compute the norm of A and the reciprocal of the condition
604 * number of A.
605 *
606  IF( notran ) THEN
607  norm = 'I'
608  ELSE
609  norm = '1'
610  END IF
611  anorm = clange( norm, n, n, a, lda, rwork )
612  CALL cgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, info )
613 *
614 * Perform refinement on each right-hand side
615 *
616  IF ( ref_type .NE. 0 ) THEN
617 
618  prec_type = ilaprec( 'D' )
619 
620  IF ( notran ) THEN
621  CALL cla_gerfsx_extended( prec_type, trans_type, n,
622  $ nrhs, a, lda, af, ldaf, ipiv, colequ, c, b,
623  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
624  $ err_bnds_comp, work, rwork, work(n+1),
625  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
626  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
627  $ info )
628  ELSE
629  CALL cla_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, rwork, work(n+1),
633  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
634  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
635  $ info )
636  END IF
637  END IF
638 
639  err_lbnd = max( 10.0, sqrt( REAL( N ) ) ) * slamch( 'Epsilon' )
640  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
641 *
642 * Compute scaled normwise condition number cond(A*C).
643 *
644  IF ( colequ .AND. notran ) THEN
645  rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
646  $ c, .true., info, work, rwork )
647  ELSE IF ( rowequ .AND. .NOT. notran ) THEN
648  rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
649  $ r, .true., info, work, rwork )
650  ELSE
651  rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
652  $ c, .false., info, work, rwork )
653  END IF
654  DO j = 1, nrhs
655 *
656 * Cap the error at 1.0.
657 *
658  IF ( n_err_bnds .GE. la_linrx_err_i
659  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0 )
660  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0
661 *
662 * Threshold the error (see LAWN).
663 *
664  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
665  err_bnds_norm( j, la_linrx_err_i ) = 1.0
666  err_bnds_norm( j, la_linrx_trust_i ) = 0.0
667  IF ( info .LE. n ) info = n + j
668  ELSE IF (err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd)
669  $ THEN
670  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
671  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
672  END IF
673 *
674 * Save the condition number.
675 *
676  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
677  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
678  END IF
679  END DO
680  END IF
681 
682  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
683 *
684 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
685 * each right-hand side using the current solution as an estimate of
686 * the true solution. If the componentwise error estimate is too
687 * large, then the solution is a lousy estimate of truth and the
688 * estimated RCOND may be too optimistic. To avoid misleading users,
689 * the inverse condition number is set to 0.0 when the estimated
690 * cwise error is at least CWISE_WRONG.
691 *
692  cwise_wrong = sqrt( slamch( 'Epsilon' ) )
693  DO j = 1, nrhs
694  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
695  $ THEN
696  rcond_tmp = cla_gercond_x( trans, n, a, lda, af, ldaf,
697  $ ipiv, x(1,j), info, work, rwork )
698  ELSE
699  rcond_tmp = 0.0
700  END IF
701 *
702 * Cap the error at 1.0.
703 *
704  IF ( n_err_bnds .GE. la_linrx_err_i
705  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0 )
706  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0
707 *
708 * Threshold the error (see LAWN).
709 *
710  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
711  err_bnds_comp( j, la_linrx_err_i ) = 1.0
712  err_bnds_comp( j, la_linrx_trust_i ) = 0.0
713  IF ( params( la_linrx_cwise_i ) .EQ. 1.0
714  $ .AND. info.LT.n + j ) info = n + j
715  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
716  $ .LT. err_lbnd ) THEN
717  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
718  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
719  END IF
720 *
721 * Save the condition number.
722 *
723  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
724  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
725  END IF
726 
727  END DO
728  END IF
729 *
730  RETURN
731 *
732 * End of CGERFSX
733 *
734  END
subroutine cla_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)
CLA_GERFSX_EXTENDED
real function cla_gercond_c(TRANS, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
CLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices...
real function cla_gercond_x(TRANS, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
CLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices...
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
subroutine cgerfsx(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, RWORK, INFO)
CGERFSX
Definition: cgerfsx.f:416
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine cgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
CGECON
Definition: cgecon.f:126