File:
[local] /
rpl /
lapack /
lapack /
zgbrfsx.f
Revision
1.4:
download - view:
text,
annotated -
select for diffs -
revision graph
Tue Dec 21 13:53:42 2010 UTC (13 years, 4 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_3,
rpl-4_1_2,
rpl-4_1_1,
rpl-4_1_0,
rpl-4_0_24,
rpl-4_0_22,
rpl-4_0_21,
rpl-4_0_20,
rpl-4_0,
HEAD
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGBRFSX( TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB,
2: $ LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND,
3: $ BERR, N_ERR_BNDS, ERR_BNDS_NORM,
4: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK,
5: $ INFO )
6: *
7: * -- LAPACK routine (version 3.2.2) --
8: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
9: * -- Jason Riedy of Univ. of California Berkeley. --
10: * -- June 2010 --
11: *
12: * -- LAPACK is a software package provided by Univ. of Tennessee, --
13: * -- Univ. of California Berkeley and NAG Ltd. --
14: *
15: IMPLICIT NONE
16: * ..
17: * .. Scalar Arguments ..
18: CHARACTER TRANS, EQUED
19: INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS,
20: $ NPARAMS, N_ERR_BNDS
21: DOUBLE PRECISION RCOND
22: * ..
23: * .. Array Arguments ..
24: INTEGER IPIV( * )
25: COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
26: $ X( LDX , * ),WORK( * )
27: DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
28: $ ERR_BNDS_NORM( NRHS, * ),
29: $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
30: * ..
31: *
32: * Purpose
33: * =======
34: *
35: * ZGBRFSX improves the computed solution to a system of linear
36: * equations and provides error bounds and backward error estimates
37: * for the solution. In addition to normwise error bound, the code
38: * provides maximum componentwise error bound if possible. See
39: * comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
40: * error bounds.
41: *
42: * The original system of linear equations may have been equilibrated
43: * before calling this routine, as described by arguments EQUED, R
44: * and C below. In this case, the solution and error bounds returned
45: * are for the original unequilibrated system.
46: *
47: * Arguments
48: * =========
49: *
50: * Some optional parameters are bundled in the PARAMS array. These
51: * settings determine how refinement is performed, but often the
52: * defaults are acceptable. If the defaults are acceptable, users
53: * can pass NPARAMS = 0 which prevents the source code from accessing
54: * the PARAMS argument.
55: *
56: * TRANS (input) CHARACTER*1
57: * Specifies the form of the system of equations:
58: * = 'N': A * X = B (No transpose)
59: * = 'T': A**T * X = B (Transpose)
60: * = 'C': A**H * X = B (Conjugate transpose = Transpose)
61: *
62: * EQUED (input) CHARACTER*1
63: * Specifies the form of equilibration that was done to A
64: * before calling this routine. This is needed to compute
65: * the solution and error bounds correctly.
66: * = 'N': No equilibration
67: * = 'R': Row equilibration, i.e., A has been premultiplied by
68: * diag(R).
69: * = 'C': Column equilibration, i.e., A has been postmultiplied
70: * by diag(C).
71: * = 'B': Both row and column equilibration, i.e., A has been
72: * replaced by diag(R) * A * diag(C).
73: * The right hand side B has been changed accordingly.
74: *
75: * N (input) INTEGER
76: * The order of the matrix A. N >= 0.
77: *
78: * KL (input) INTEGER
79: * The number of subdiagonals within the band of A. KL >= 0.
80: *
81: * KU (input) INTEGER
82: * The number of superdiagonals within the band of A. KU >= 0.
83: *
84: * NRHS (input) INTEGER
85: * The number of right hand sides, i.e., the number of columns
86: * of the matrices B and X. NRHS >= 0.
87: *
88: * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
89: * The original band matrix A, stored in rows 1 to KL+KU+1.
90: * The j-th column of A is stored in the j-th column of the
91: * array AB as follows:
92: * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
93: *
94: * LDAB (input) INTEGER
95: * The leading dimension of the array AB. LDAB >= KL+KU+1.
96: *
97: * AFB (input) DOUBLE PRECISION array, dimension (LDAFB,N)
98: * Details of the LU factorization of the band matrix A, as
99: * computed by DGBTRF. U is stored as an upper triangular band
100: * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
101: * the multipliers used during the factorization are stored in
102: * rows KL+KU+2 to 2*KL+KU+1.
103: *
104: * LDAFB (input) INTEGER
105: * The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1.
106: *
107: * IPIV (input) INTEGER array, dimension (N)
108: * The pivot indices from DGETRF; for 1<=i<=N, row i of the
109: * matrix was interchanged with row IPIV(i).
110: *
111: * R (input or output) DOUBLE PRECISION array, dimension (N)
112: * The row scale factors for A. If EQUED = 'R' or 'B', A is
113: * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
114: * is not accessed. R is an input argument if FACT = 'F';
115: * otherwise, R is an output argument. If FACT = 'F' and
116: * EQUED = 'R' or 'B', each element of R must be positive.
117: * If R is output, each element of R is a power of the radix.
118: * If R is input, each element of R should be a power of the radix
119: * to ensure a reliable solution and error estimates. Scaling by
120: * powers of the radix does not cause rounding errors unless the
121: * result underflows or overflows. Rounding errors during scaling
122: * lead to refining with a matrix that is not equivalent to the
123: * input matrix, producing error estimates that may not be
124: * reliable.
125: *
126: * C (input or output) DOUBLE PRECISION array, dimension (N)
127: * The column scale factors for A. If EQUED = 'C' or 'B', A is
128: * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
129: * is not accessed. C is an input argument if FACT = 'F';
130: * otherwise, C is an output argument. If FACT = 'F' and
131: * EQUED = 'C' or 'B', each element of C must be positive.
132: * If C is output, each element of C is a power of the radix.
133: * If C is input, each element of C should be a power of the radix
134: * to ensure a reliable solution and error estimates. Scaling by
135: * powers of the radix does not cause rounding errors unless the
136: * result underflows or overflows. Rounding errors during scaling
137: * lead to refining with a matrix that is not equivalent to the
138: * input matrix, producing error estimates that may not be
139: * reliable.
140: *
141: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
142: * The right hand side matrix B.
143: *
144: * LDB (input) INTEGER
145: * The leading dimension of the array B. LDB >= max(1,N).
146: *
147: * X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS)
148: * On entry, the solution matrix X, as computed by DGETRS.
149: * On exit, the improved solution matrix X.
150: *
151: * LDX (input) INTEGER
152: * The leading dimension of the array X. LDX >= max(1,N).
153: *
154: * RCOND (output) DOUBLE PRECISION
155: * Reciprocal scaled condition number. This is an estimate of the
156: * reciprocal Skeel condition number of the matrix A after
157: * equilibration (if done). If this is less than the machine
158: * precision (in particular, if it is zero), the matrix is singular
159: * to working precision. Note that the error may still be small even
160: * if this number is very small and the matrix appears ill-
161: * conditioned.
162: *
163: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
164: * Componentwise relative backward error. This is the
165: * componentwise relative backward error of each solution vector X(j)
166: * (i.e., the smallest relative change in any element of A or B that
167: * makes X(j) an exact solution).
168: *
169: * N_ERR_BNDS (input) INTEGER
170: * Number of error bounds to return for each right hand side
171: * and each type (normwise or componentwise). See ERR_BNDS_NORM and
172: * ERR_BNDS_COMP below.
173: *
174: * ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
175: * For each right-hand side, this array contains information about
176: * various error bounds and condition numbers corresponding to the
177: * normwise relative error, which is defined as follows:
178: *
179: * Normwise relative error in the ith solution vector:
180: * max_j (abs(XTRUE(j,i) - X(j,i)))
181: * ------------------------------
182: * max_j abs(X(j,i))
183: *
184: * The array is indexed by the type of error information as described
185: * below. There currently are up to three pieces of information
186: * returned.
187: *
188: * The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
189: * right-hand side.
190: *
191: * The second index in ERR_BNDS_NORM(:,err) contains the following
192: * three fields:
193: * err = 1 "Trust/don't trust" boolean. Trust the answer if the
194: * reciprocal condition number is less than the threshold
195: * sqrt(n) * dlamch('Epsilon').
196: *
197: * err = 2 "Guaranteed" error bound: The estimated forward error,
198: * almost certainly within a factor of 10 of the true error
199: * so long as the next entry is greater than the threshold
200: * sqrt(n) * dlamch('Epsilon'). This error bound should only
201: * be trusted if the previous boolean is true.
202: *
203: * err = 3 Reciprocal condition number: Estimated normwise
204: * reciprocal condition number. Compared with the threshold
205: * sqrt(n) * dlamch('Epsilon') to determine if the error
206: * estimate is "guaranteed". These reciprocal condition
207: * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
208: * appropriately scaled matrix Z.
209: * Let Z = S*A, where S scales each row by a power of the
210: * radix so all absolute row sums of Z are approximately 1.
211: *
212: * See Lapack Working Note 165 for further details and extra
213: * cautions.
214: *
215: * ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
216: * For each right-hand side, this array contains information about
217: * various error bounds and condition numbers corresponding to the
218: * componentwise relative error, which is defined as follows:
219: *
220: * Componentwise relative error in the ith solution vector:
221: * abs(XTRUE(j,i) - X(j,i))
222: * max_j ----------------------
223: * abs(X(j,i))
224: *
225: * The array is indexed by the right-hand side i (on which the
226: * componentwise relative error depends), and the type of error
227: * information as described below. There currently are up to three
228: * pieces of information returned for each right-hand side. If
229: * componentwise accuracy is not requested (PARAMS(3) = 0.0), then
230: * ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
231: * the first (:,N_ERR_BNDS) entries are returned.
232: *
233: * The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
234: * right-hand side.
235: *
236: * The second index in ERR_BNDS_COMP(:,err) contains the following
237: * three fields:
238: * err = 1 "Trust/don't trust" boolean. Trust the answer if the
239: * reciprocal condition number is less than the threshold
240: * sqrt(n) * dlamch('Epsilon').
241: *
242: * err = 2 "Guaranteed" error bound: The estimated forward error,
243: * almost certainly within a factor of 10 of the true error
244: * so long as the next entry is greater than the threshold
245: * sqrt(n) * dlamch('Epsilon'). This error bound should only
246: * be trusted if the previous boolean is true.
247: *
248: * err = 3 Reciprocal condition number: Estimated componentwise
249: * reciprocal condition number. Compared with the threshold
250: * sqrt(n) * dlamch('Epsilon') to determine if the error
251: * estimate is "guaranteed". These reciprocal condition
252: * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
253: * appropriately scaled matrix Z.
254: * Let Z = S*(A*diag(x)), where x is the solution for the
255: * current right-hand side and S scales each row of
256: * A*diag(x) by a power of the radix so all absolute row
257: * sums of Z are approximately 1.
258: *
259: * See Lapack Working Note 165 for further details and extra
260: * cautions.
261: *
262: * NPARAMS (input) INTEGER
263: * Specifies the number of parameters set in PARAMS. If .LE. 0, the
264: * PARAMS array is never referenced and default values are used.
265: *
266: * PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS
267: * Specifies algorithm parameters. If an entry is .LT. 0.0, then
268: * that entry will be filled with default value used for that
269: * parameter. Only positions up to NPARAMS are accessed; defaults
270: * are used for higher-numbered parameters.
271: *
272: * PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
273: * refinement or not.
274: * Default: 1.0D+0
275: * = 0.0 : No refinement is performed, and no error bounds are
276: * computed.
277: * = 1.0 : Use the double-precision refinement algorithm,
278: * possibly with doubled-single computations if the
279: * compilation environment does not support DOUBLE
280: * PRECISION.
281: * (other values are reserved for future use)
282: *
283: * PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
284: * computations allowed for refinement.
285: * Default: 10
286: * Aggressive: Set to 100 to permit convergence using approximate
287: * factorizations or factorizations other than LU. If
288: * the factorization uses a technique other than
289: * Gaussian elimination, the guarantees in
290: * err_bnds_norm and err_bnds_comp may no longer be
291: * trustworthy.
292: *
293: * PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
294: * will attempt to find a solution with small componentwise
295: * relative error in the double-precision algorithm. Positive
296: * is true, 0.0 is false.
297: * Default: 1.0 (attempt componentwise convergence)
298: *
299: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
300: *
301: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
302: *
303: * INFO (output) INTEGER
304: * = 0: Successful exit. The solution to every right-hand side is
305: * guaranteed.
306: * < 0: If INFO = -i, the i-th argument had an illegal value
307: * > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
308: * has been completed, but the factor U is exactly singular, so
309: * the solution and error bounds could not be computed. RCOND = 0
310: * is returned.
311: * = N+J: The solution corresponding to the Jth right-hand side is
312: * not guaranteed. The solutions corresponding to other right-
313: * hand sides K with K > J may not be guaranteed as well, but
314: * only the first such right-hand side is reported. If a small
315: * componentwise error is not requested (PARAMS(3) = 0.0) then
316: * the Jth right-hand side is the first with a normwise error
317: * bound that is not guaranteed (the smallest J such
318: * that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
319: * the Jth right-hand side is the first with either a normwise or
320: * componentwise error bound that is not guaranteed (the smallest
321: * J such that either ERR_BNDS_NORM(J,1) = 0.0 or
322: * ERR_BNDS_COMP(J,1) = 0.0). See the definition of
323: * ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
324: * about all of the right-hand sides check ERR_BNDS_NORM or
325: * ERR_BNDS_COMP.
326: *
327: * ==================================================================
328: *
329: * .. Parameters ..
330: DOUBLE PRECISION ZERO, ONE
331: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
332: DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
333: DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
334: DOUBLE PRECISION DZTHRESH_DEFAULT
335: PARAMETER ( ITREF_DEFAULT = 1.0D+0 )
336: PARAMETER ( ITHRESH_DEFAULT = 10.0D+0 )
337: PARAMETER ( COMPONENTWISE_DEFAULT = 1.0D+0 )
338: PARAMETER ( RTHRESH_DEFAULT = 0.5D+0 )
339: PARAMETER ( DZTHRESH_DEFAULT = 0.25D+0 )
340: INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
341: $ LA_LINRX_CWISE_I
342: PARAMETER ( LA_LINRX_ITREF_I = 1,
343: $ LA_LINRX_ITHRESH_I = 2 )
344: PARAMETER ( LA_LINRX_CWISE_I = 3 )
345: INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
346: $ LA_LINRX_RCOND_I
347: PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
348: PARAMETER ( LA_LINRX_RCOND_I = 3 )
349: * ..
350: * .. Local Scalars ..
351: CHARACTER(1) NORM
352: LOGICAL ROWEQU, COLEQU, NOTRAN, IGNORE_CWISE
353: INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE, N_NORMS,
354: $ ITHRESH
355: DOUBLE PRECISION ANORM, RCOND_TMP, ILLRCOND_THRESH, ERR_LBND,
356: $ CWISE_WRONG, RTHRESH, UNSTABLE_THRESH
357: * ..
358: * .. External Subroutines ..
359: EXTERNAL XERBLA, ZGBCON, ZLA_GBRFSX_EXTENDED
360: * ..
361: * .. Intrinsic Functions ..
362: INTRINSIC MAX, SQRT, TRANSFER
363: * ..
364: * .. External Functions ..
365: EXTERNAL LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC
366: EXTERNAL DLAMCH, ZLANGB, ZLA_GBRCOND_X, ZLA_GBRCOND_C
367: DOUBLE PRECISION DLAMCH, ZLANGB, ZLA_GBRCOND_X, ZLA_GBRCOND_C
368: LOGICAL LSAME
369: INTEGER BLAS_FPINFO_X
370: INTEGER ILATRANS, ILAPREC
371: * ..
372: * .. Executable Statements ..
373: *
374: * Check the input parameters.
375: *
376: INFO = 0
377: TRANS_TYPE = ILATRANS( TRANS )
378: REF_TYPE = INT( ITREF_DEFAULT )
379: IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
380: IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN
381: PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
382: ELSE
383: REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
384: END IF
385: END IF
386: *
387: * Set default parameters.
388: *
389: ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' )
390: ITHRESH = INT( ITHRESH_DEFAULT )
391: RTHRESH = RTHRESH_DEFAULT
392: UNSTABLE_THRESH = DZTHRESH_DEFAULT
393: IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0
394: *
395: IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
396: IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN
397: PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
398: ELSE
399: ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
400: END IF
401: END IF
402: IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
403: IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN
404: IF ( IGNORE_CWISE ) THEN
405: PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0
406: ELSE
407: PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0
408: END IF
409: ELSE
410: IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0
411: END IF
412: END IF
413: IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
414: N_NORMS = 0
415: ELSE IF ( IGNORE_CWISE ) THEN
416: N_NORMS = 1
417: ELSE
418: N_NORMS = 2
419: END IF
420: *
421: NOTRAN = LSAME( TRANS, 'N' )
422: ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
423: COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
424: *
425: * Test input parameters.
426: *
427: IF( TRANS_TYPE.EQ.-1 ) THEN
428: INFO = -1
429: ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND.
430: $ .NOT.LSAME( EQUED, 'N' ) ) THEN
431: INFO = -2
432: ELSE IF( N.LT.0 ) THEN
433: INFO = -3
434: ELSE IF( KL.LT.0 ) THEN
435: INFO = -4
436: ELSE IF( KU.LT.0 ) THEN
437: INFO = -5
438: ELSE IF( NRHS.LT.0 ) THEN
439: INFO = -6
440: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
441: INFO = -8
442: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
443: INFO = -10
444: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
445: INFO = -13
446: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
447: INFO = -15
448: END IF
449: IF( INFO.NE.0 ) THEN
450: CALL XERBLA( 'ZGBRFSX', -INFO )
451: RETURN
452: END IF
453: *
454: * Quick return if possible.
455: *
456: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
457: RCOND = 1.0D+0
458: DO J = 1, NRHS
459: BERR( J ) = 0.0D+0
460: IF ( N_ERR_BNDS .GE. 1 ) THEN
461: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
462: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
463: END IF
464: IF ( N_ERR_BNDS .GE. 2 ) THEN
465: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0
466: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0
467: END IF
468: IF ( N_ERR_BNDS .GE. 3 ) THEN
469: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0
470: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0
471: END IF
472: END DO
473: RETURN
474: END IF
475: *
476: * Default to failure.
477: *
478: RCOND = 0.0D+0
479: DO J = 1, NRHS
480: BERR( J ) = 1.0D+0
481: IF ( N_ERR_BNDS .GE. 1 ) THEN
482: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
483: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
484: END IF
485: IF ( N_ERR_BNDS .GE. 2 ) THEN
486: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
487: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
488: END IF
489: IF ( N_ERR_BNDS .GE. 3 ) THEN
490: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0
491: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0
492: END IF
493: END DO
494: *
495: * Compute the norm of A and the reciprocal of the condition
496: * number of A.
497: *
498: IF( NOTRAN ) THEN
499: NORM = 'I'
500: ELSE
501: NORM = '1'
502: END IF
503: ANORM = ZLANGB( NORM, N, KL, KU, AB, LDAB, RWORK )
504: CALL ZGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND,
505: $ WORK, RWORK, INFO )
506: *
507: * Perform refinement on each right-hand side
508: *
509: IF ( REF_TYPE .NE. 0 ) THEN
510:
511: PREC_TYPE = ILAPREC( 'E' )
512:
513: IF ( NOTRAN ) THEN
514: CALL ZLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
515: $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B,
516: $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
517: $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
518: $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
519: $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
520: $ INFO )
521: ELSE
522: CALL ZLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
523: $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, ROWEQU, R, B,
524: $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
525: $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
526: $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
527: $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
528: $ INFO )
529: END IF
530: END IF
531:
532: ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' )
533: IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1) THEN
534: *
535: * Compute scaled normwise condition number cond(A*C).
536: *
537: IF ( COLEQU .AND. NOTRAN ) THEN
538: RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
539: $ LDAFB, IPIV, C, .TRUE., INFO, WORK, RWORK )
540: ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN
541: RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
542: $ LDAFB, IPIV, R, .TRUE., INFO, WORK, RWORK )
543: ELSE
544: RCOND_TMP = ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
545: $ LDAFB, IPIV, C, .FALSE., INFO, WORK, RWORK )
546: END IF
547: DO J = 1, NRHS
548: *
549: * Cap the error at 1.0.
550: *
551: IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
552: $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0)
553: $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
554: *
555: * Threshold the error (see LAWN).
556: *
557: IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
558: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
559: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0
560: IF ( INFO .LE. N ) INFO = N + J
561: ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
562: $ THEN
563: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
564: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
565: END IF
566: *
567: * Save the condition number.
568: *
569: IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
570: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
571: END IF
572:
573: END DO
574: END IF
575:
576: IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN
577: *
578: * Compute componentwise condition number cond(A*diag(Y(:,J))) for
579: * each right-hand side using the current solution as an estimate of
580: * the true solution. If the componentwise error estimate is too
581: * large, then the solution is a lousy estimate of truth and the
582: * estimated RCOND may be too optimistic. To avoid misleading users,
583: * the inverse condition number is set to 0.0 when the estimated
584: * cwise error is at least CWISE_WRONG.
585: *
586: CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) )
587: DO J = 1, NRHS
588: IF (ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
589: $ THEN
590: RCOND_TMP = ZLA_GBRCOND_X( TRANS, N, KL, KU, AB, LDAB,
591: $ AFB, LDAFB, IPIV, X( 1, J ), INFO, WORK, RWORK )
592: ELSE
593: RCOND_TMP = 0.0D+0
594: END IF
595: *
596: * Cap the error at 1.0.
597: *
598: IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
599: $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
600: $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
601: *
602: * Threshold the error (see LAWN).
603: *
604: IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
605: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
606: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0
607: IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0
608: $ .AND. INFO.LT.N + J ) INFO = N + J
609: ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
610: $ .LT. ERR_LBND ) THEN
611: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
612: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
613: END IF
614: *
615: * Save the condition number.
616: *
617: IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
618: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
619: END IF
620:
621: END DO
622: END IF
623: *
624: RETURN
625: *
626: * End of ZGBRFSX
627: *
628: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgbrfsx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack