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