LAPACK  3.10.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)
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 < 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 <= 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 < 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 *> \ingroup complexGEcomputational
408 *
409 * =====================================================================
410  SUBROUTINE cgerfsx( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
411  $ R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
412  $ ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
413  $ WORK, RWORK, INFO )
414 *
415 * -- LAPACK computational routine --
416 * -- LAPACK is a software package provided by Univ. of Tennessee, --
417 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
418 *
419 * .. Scalar Arguments ..
420  CHARACTER TRANS, EQUED
421  INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
422  $ N_ERR_BNDS
423  REAL RCOND
424 * ..
425 * .. Array Arguments ..
426  INTEGER IPIV( * )
427  COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
428  $ X( LDX , * ), WORK( * )
429  REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
430  $ err_bnds_norm( nrhs, * ),
431  $ err_bnds_comp( nrhs, * ), rwork( * )
432 * ..
433 *
434 * ==================================================================
435 *
436 * .. Parameters ..
437  REAL ZERO, ONE
438  PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
439  REAL ITREF_DEFAULT, ITHRESH_DEFAULT,
440  $ componentwise_default
441  REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT
442  parameter( itref_default = 1.0 )
443  parameter( ithresh_default = 10.0 )
444  parameter( componentwise_default = 1.0 )
445  parameter( rthresh_default = 0.5 )
446  parameter( dzthresh_default = 0.25 )
447  INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
448  $ la_linrx_cwise_i
449  parameter( la_linrx_itref_i = 1,
450  $ la_linrx_ithresh_i = 2 )
451  parameter( la_linrx_cwise_i = 3 )
452  INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
453  $ LA_LINRX_RCOND_I
454  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
455  parameter( la_linrx_rcond_i = 3 )
456 * ..
457 * .. Local Scalars ..
458  CHARACTER(1) NORM
459  LOGICAL ROWEQU, COLEQU, NOTRAN
460  INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
461  INTEGER N_NORMS
462  REAL ANORM, RCOND_TMP
463  REAL ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
464  LOGICAL IGNORE_CWISE
465  INTEGER ITHRESH
466  REAL RTHRESH, UNSTABLE_THRESH
467 * ..
468 * .. External Subroutines ..
470 * ..
471 * .. Intrinsic Functions ..
472  INTRINSIC max, sqrt, transfer
473 * ..
474 * .. External Functions ..
475  EXTERNAL lsame, ilatrans, ilaprec
477  REAL SLAMCH, CLANGE, CLA_GERCOND_X, CLA_GERCOND_C
478  LOGICAL LSAME
479  INTEGER ILATRANS, ILAPREC
480 * ..
481 * .. Executable Statements ..
482 *
483 * Check the input parameters.
484 *
485  info = 0
486  trans_type = ilatrans( trans )
487  ref_type = int( itref_default )
488  IF ( nparams .GE. la_linrx_itref_i ) THEN
489  IF ( params( la_linrx_itref_i ) .LT. 0.0 ) THEN
490  params( la_linrx_itref_i ) = itref_default
491  ELSE
492  ref_type = params( la_linrx_itref_i )
493  END IF
494  END IF
495 *
496 * Set default parameters.
497 *
498  illrcond_thresh = real( n ) * slamch( 'Epsilon' )
499  ithresh = int( ithresh_default )
500  rthresh = rthresh_default
501  unstable_thresh = dzthresh_default
502  ignore_cwise = componentwise_default .EQ. 0.0
503 *
504  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
505  IF ( params( la_linrx_ithresh_i ).LT.0.0 ) THEN
506  params(la_linrx_ithresh_i) = ithresh
507  ELSE
508  ithresh = int( params( la_linrx_ithresh_i ) )
509  END IF
510  END IF
511  IF ( nparams.GE.la_linrx_cwise_i ) THEN
512  IF ( params( la_linrx_cwise_i ).LT.0.0 ) THEN
513  IF ( ignore_cwise ) THEN
514  params( la_linrx_cwise_i ) = 0.0
515  ELSE
516  params( la_linrx_cwise_i ) = 1.0
517  END IF
518  ELSE
519  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0
520  END IF
521  END IF
522  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
523  n_norms = 0
524  ELSE IF ( ignore_cwise ) THEN
525  n_norms = 1
526  ELSE
527  n_norms = 2
528  END IF
529 *
530  notran = lsame( trans, 'N' )
531  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
532  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
533 *
534 * Test input parameters.
535 *
536  IF( trans_type.EQ.-1 ) THEN
537  info = -1
538  ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
539  $ .NOT.lsame( equed, 'N' ) ) THEN
540  info = -2
541  ELSE IF( n.LT.0 ) THEN
542  info = -3
543  ELSE IF( nrhs.LT.0 ) THEN
544  info = -4
545  ELSE IF( lda.LT.max( 1, n ) ) THEN
546  info = -6
547  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
548  info = -8
549  ELSE IF( ldb.LT.max( 1, n ) ) THEN
550  info = -13
551  ELSE IF( ldx.LT.max( 1, n ) ) THEN
552  info = -15
553  END IF
554  IF( info.NE.0 ) THEN
555  CALL xerbla( 'CGERFSX', -info )
556  RETURN
557  END IF
558 *
559 * Quick return if possible.
560 *
561  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
562  rcond = 1.0
563  DO j = 1, nrhs
564  berr( j ) = 0.0
565  IF ( n_err_bnds .GE. 1 ) THEN
566  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
567  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
568  END IF
569  IF ( n_err_bnds .GE. 2 ) THEN
570  err_bnds_norm( j, la_linrx_err_i ) = 0.0
571  err_bnds_comp( j, la_linrx_err_i ) = 0.0
572  END IF
573  IF ( n_err_bnds .GE. 3 ) THEN
574  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0
575  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
576  END IF
577  END DO
578  RETURN
579  END IF
580 *
581 * Default to failure.
582 *
583  rcond = 0.0
584  DO j = 1, nrhs
585  berr( j ) = 1.0
586  IF ( n_err_bnds .GE. 1 ) THEN
587  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
588  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
589  END IF
590  IF ( n_err_bnds .GE. 2 ) THEN
591  err_bnds_norm( j, la_linrx_err_i ) = 1.0
592  err_bnds_comp( j, la_linrx_err_i ) = 1.0
593  END IF
594  IF ( n_err_bnds .GE. 3 ) THEN
595  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
596  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0
597  END IF
598  END DO
599 *
600 * Compute the norm of A and the reciprocal of the condition
601 * number of A.
602 *
603  IF( notran ) THEN
604  norm = 'I'
605  ELSE
606  norm = '1'
607  END IF
608  anorm = clange( norm, n, n, a, lda, rwork )
609  CALL cgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, info )
610 *
611 * Perform refinement on each right-hand side
612 *
613  IF ( ref_type .NE. 0 ) THEN
614 
615  prec_type = ilaprec( 'D' )
616 
617  IF ( notran ) THEN
618  CALL cla_gerfsx_extended( prec_type, trans_type, n,
619  $ nrhs, a, lda, af, ldaf, ipiv, colequ, c, b,
620  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
621  $ err_bnds_comp, work, rwork, work(n+1),
622  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
623  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
624  $ info )
625  ELSE
626  CALL cla_gerfsx_extended( prec_type, trans_type, n,
627  $ nrhs, a, lda, af, ldaf, ipiv, rowequ, r, b,
628  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
629  $ err_bnds_comp, work, rwork, work(n+1),
630  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
631  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
632  $ info )
633  END IF
634  END IF
635 
636  err_lbnd = max( 10.0, sqrt( real( n ) ) ) * slamch( 'Epsilon' )
637  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
638 *
639 * Compute scaled normwise condition number cond(A*C).
640 *
641  IF ( colequ .AND. notran ) THEN
642  rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
643  $ c, .true., info, work, rwork )
644  ELSE IF ( rowequ .AND. .NOT. notran ) THEN
645  rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
646  $ r, .true., info, work, rwork )
647  ELSE
648  rcond_tmp = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
649  $ c, .false., info, work, rwork )
650  END IF
651  DO j = 1, nrhs
652 *
653 * Cap the error at 1.0.
654 *
655  IF ( n_err_bnds .GE. la_linrx_err_i
656  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0 )
657  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0
658 *
659 * Threshold the error (see LAWN).
660 *
661  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
662  err_bnds_norm( j, la_linrx_err_i ) = 1.0
663  err_bnds_norm( j, la_linrx_trust_i ) = 0.0
664  IF ( info .LE. n ) info = n + j
665  ELSE IF (err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd)
666  $ THEN
667  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
668  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
669  END IF
670 *
671 * Save the condition number.
672 *
673  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
674  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
675  END IF
676  END DO
677  END IF
678 
679  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
680 *
681 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
682 * each right-hand side using the current solution as an estimate of
683 * the true solution. If the componentwise error estimate is too
684 * large, then the solution is a lousy estimate of truth and the
685 * estimated RCOND may be too optimistic. To avoid misleading users,
686 * the inverse condition number is set to 0.0 when the estimated
687 * cwise error is at least CWISE_WRONG.
688 *
689  cwise_wrong = sqrt( slamch( 'Epsilon' ) )
690  DO j = 1, nrhs
691  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
692  $ THEN
693  rcond_tmp = cla_gercond_x( trans, n, a, lda, af, ldaf,
694  $ ipiv, x(1,j), info, work, rwork )
695  ELSE
696  rcond_tmp = 0.0
697  END IF
698 *
699 * Cap the error at 1.0.
700 *
701  IF ( n_err_bnds .GE. la_linrx_err_i
702  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0 )
703  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0
704 *
705 * Threshold the error (see LAWN).
706 *
707  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
708  err_bnds_comp( j, la_linrx_err_i ) = 1.0
709  err_bnds_comp( j, la_linrx_trust_i ) = 0.0
710  IF ( params( la_linrx_cwise_i ) .EQ. 1.0
711  $ .AND. info.LT.n + j ) info = n + j
712  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
713  $ .LT. err_lbnd ) THEN
714  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
715  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
716  END IF
717 *
718 * Save the condition number.
719 *
720  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
721  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
722  END IF
723 
724  END DO
725  END IF
726 *
727  RETURN
728 *
729 * End of CGERFSX
730 *
731  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:58
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:58
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:115
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:414
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.
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
subroutine cgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
CGECON
Definition: cgecon.f:124
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 slamch(CMACH)
SLAMCH
Definition: slamch.f:68