LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
zgerfsx.f
Go to the documentation of this file.
1 *> \brief \b ZGERFSX
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/zgerfsx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgerfsx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgerfsx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGERFSX( 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 * DOUBLE PRECISION RCOND
31 * ..
32 * .. Array Arguments ..
33 * INTEGER IPIV( * )
34 * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
35 * \$ X( LDX , * ), WORK( * )
36 * DOUBLE PRECISION 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 *> ZGERFSX 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*16 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*16 array, dimension (LDAF,N)
124 *> The factors L and U from the factorization A = P*L*U
125 *> as computed by ZGETRF.
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 ZGETRF; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 array, dimension (LDX,NRHS)
186 *> On entry, the solution matrix X, as computed by ZGETRS.
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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) * dlamch('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) * dlamch('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) * dlamch('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 DOUBLE PRECISION 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.0D+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*16 array, dimension (2*N)
365 *> \endverbatim
366 *>
367 *> \param[out] RWORK
368 *> \verbatim
369 *> RWORK is DOUBLE PRECISION 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 November 2011
408 *
409 *> \ingroup complex16GEcomputational
410 *
411 * =====================================================================
412  SUBROUTINE zgerfsx( 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.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  DOUBLE PRECISION RCOND
427 * ..
428 * .. Array Arguments ..
429  INTEGER IPIV( * )
430  COMPLEX*16 A( lda, * ), AF( ldaf, * ), B( ldb, * ),
431  \$ x( ldx , * ), work( * )
432  DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
433  \$ err_bnds_norm( nrhs, * ),
434  \$ err_bnds_comp( nrhs, * ), rwork( * )
435 * ..
436 *
437 * ==================================================================
438 *
439 * .. Parameters ..
440  DOUBLE PRECISION ZERO, ONE
441  parameter ( zero = 0.0d+0, one = 1.0d+0 )
442  DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
443  DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
444  DOUBLE PRECISION DZTHRESH_DEFAULT
445  parameter ( itref_default = 1.0d+0 )
446  parameter ( ithresh_default = 10.0d+0 )
447  parameter ( componentwise_default = 1.0d+0 )
448  parameter ( rthresh_default = 0.5d+0 )
449  parameter ( dzthresh_default = 0.25d+0 )
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  DOUBLE PRECISION ANORM, RCOND_TMP
466  DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
467  LOGICAL IGNORE_CWISE
468  INTEGER ITHRESH
469  DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH
470 * ..
471 * .. External Subroutines ..
473 * ..
474 * .. Intrinsic Functions ..
475  INTRINSIC max, sqrt, transfer
476 * ..
477 * .. External Functions ..
478  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
480  DOUBLE PRECISION DLAMCH, ZLANGE, ZLA_GERCOND_X, ZLA_GERCOND_C
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.0d+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 = dble( n ) * dlamch( 'Epsilon' )
503  ithresh = int( ithresh_default )
504  rthresh = rthresh_default
505  unstable_thresh = dzthresh_default
506  ignore_cwise = componentwise_default .EQ. 0.0d+0
507 *
508  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
509  IF ( params( la_linrx_ithresh_i ).LT.0.0d+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.0d+0 ) THEN
517  IF ( ignore_cwise ) THEN
518  params( la_linrx_cwise_i ) = 0.0d+0
519  ELSE
520  params( la_linrx_cwise_i ) = 1.0d+0
521  END IF
522  ELSE
523  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+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( 'ZGERFSX', -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.0d+0
567  DO j = 1, nrhs
568  berr( j ) = 0.0d+0
569  IF ( n_err_bnds .GE. 1 ) THEN
570  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
571  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
572  END IF
573  IF ( n_err_bnds .GE. 2 ) THEN
574  err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
575  err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
576  END IF
577  IF ( n_err_bnds .GE. 3 ) THEN
578  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
579  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
580  END IF
581  END DO
582  RETURN
583  END IF
584 *
585 * Default to failure.
586 *
587  rcond = 0.0d+0
588  DO j = 1, nrhs
589  berr( j ) = 1.0d+0
590  IF ( n_err_bnds .GE. 1 ) THEN
591  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
592  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
593  END IF
594  IF ( n_err_bnds .GE. 2 ) THEN
595  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
596  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
597  END IF
598  IF ( n_err_bnds .GE. 3 ) THEN
599  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
600  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+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 = zlange( norm, n, n, a, lda, rwork )
613  CALL zgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, info )
614 *
615 * Perform refinement on each right-hand side
616 *
617  IF ( ref_type .NE. 0 ) THEN
618
619  prec_type = ilaprec( 'E' )
620
621  IF ( notran ) THEN
622  CALL zla_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, rwork, work(n+1),
626  \$ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
627  \$ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
628  \$ info )
629  ELSE
630  CALL zla_gerfsx_extended( prec_type, trans_type, n,
631  \$ nrhs, a, lda, af, ldaf, ipiv, rowequ, r, b,
632  \$ ldb, x, ldx, berr, n_norms, err_bnds_norm,
633  \$ err_bnds_comp, work, rwork, work(n+1),
634  \$ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
635  \$ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
636  \$ info )
637  END IF
638  END IF
639
640  err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
641  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
642 *
643 * Compute scaled normwise condition number cond(A*C).
644 *
645  IF ( colequ .AND. notran ) THEN
646  rcond_tmp = zla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
647  \$ c, .true., info, work, rwork )
648  ELSE IF ( rowequ .AND. .NOT. notran ) THEN
649  rcond_tmp = zla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
650  \$ r, .true., info, work, rwork )
651  ELSE
652  rcond_tmp = zla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
653  \$ c, .false., info, work, rwork )
654  END IF
655  DO j = 1, nrhs
656 *
657 * Cap the error at 1.0.
658 *
659  IF ( n_err_bnds .GE. la_linrx_err_i
660  \$ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0 )
661  \$ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
662 *
663 * Threshold the error (see LAWN).
664 *
665  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
666  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
667  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
668  IF ( info .LE. n ) info = n + j
669  ELSE IF (err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd)
670  \$ THEN
671  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
672  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
673  END IF
674 *
675 * Save the condition number.
676 *
677  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
678  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
679  END IF
680  END DO
681  END IF
682
683  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
684 *
685 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
686 * each right-hand side using the current solution as an estimate of
687 * the true solution. If the componentwise error estimate is too
688 * large, then the solution is a lousy estimate of truth and the
689 * estimated RCOND may be too optimistic. To avoid misleading users,
690 * the inverse condition number is set to 0.0 when the estimated
691 * cwise error is at least CWISE_WRONG.
692 *
693  cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
694  DO j = 1, nrhs
695  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
696  \$ THEN
697  rcond_tmp = zla_gercond_x( trans, n, a, lda, af, ldaf,
698  \$ ipiv, x(1,j), info, work, rwork )
699  ELSE
700  rcond_tmp = 0.0d+0
701  END IF
702 *
703 * Cap the error at 1.0.
704 *
705  IF ( n_err_bnds .GE. la_linrx_err_i
706  \$ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
707  \$ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
708 *
709 * Threshold the error (see LAWN).
710 *
711  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
712  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
713  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
714  IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
715  \$ .AND. info.LT.n + j ) info = n + j
716  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
717  \$ .LT. err_lbnd ) THEN
718  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
719  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
720  END IF
721 *
722 * Save the condition number.
723 *
724  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
725  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
726  END IF
727
728  END DO
729  END IF
730 *
731  RETURN
732 *
733 * End of ZGERFSX
734 *
735  END
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
double precision function zla_gercond_c(TRANS, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
ZLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices...
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zgerfsx(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)
ZGERFSX
Definition: zgerfsx.f:416
subroutine zla_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)
ZLA_GERFSX_EXTENDED
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine zgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZGECON
Definition: zgecon.f:126
double precision function zla_gercond_x(TRANS, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
ZLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices...
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55