Annotation of rpl/lapack/lapack/dgbrfsx.f, revision 1.2
1.1 bertrand 1: SUBROUTINE DGBRFSX( 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, IWORK,
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( * ), IWORK( * )
25: DOUBLE PRECISION 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, * )
30: * ..
31: *
32: * Purpose
33: * =======
34: *
35: * DGBRFSX 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) DOUBLE PRECISION array, dimension (4*N)
300: *
301: * IWORK (workspace) INTEGER array, dimension (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
353: INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
354: INTEGER N_NORMS
355: DOUBLE PRECISION ANORM, RCOND_TMP
356: DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
357: LOGICAL IGNORE_CWISE
358: INTEGER ITHRESH
359: DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH
360: * ..
361: * .. External Subroutines ..
362: EXTERNAL XERBLA, DGBCON
363: EXTERNAL DLA_GBRFSX_EXTENDED
364: * ..
365: * .. Intrinsic Functions ..
366: INTRINSIC MAX, SQRT
367: * ..
368: * .. External Functions ..
369: EXTERNAL LSAME, BLAS_FPINFO_X, ILATRANS, ILAPREC
370: EXTERNAL DLAMCH, DLANGB, DLA_GBRCOND
371: DOUBLE PRECISION DLAMCH, DLANGB, DLA_GBRCOND
372: LOGICAL LSAME
373: INTEGER BLAS_FPINFO_X
374: INTEGER ILATRANS, ILAPREC
375: * ..
376: * .. Executable Statements ..
377: *
378: * Check the input parameters.
379: *
380: INFO = 0
381: TRANS_TYPE = ILATRANS( TRANS )
382: REF_TYPE = INT( ITREF_DEFAULT )
383: IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
384: IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0D+0 ) THEN
385: PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
386: ELSE
387: REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
388: END IF
389: END IF
390: *
391: * Set default parameters.
392: *
393: ILLRCOND_THRESH = DBLE( N ) * DLAMCH( 'Epsilon' )
394: ITHRESH = INT( ITHRESH_DEFAULT )
395: RTHRESH = RTHRESH_DEFAULT
396: UNSTABLE_THRESH = DZTHRESH_DEFAULT
397: IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0D+0
398: *
399: IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
400: IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0D+0 ) THEN
401: PARAMS( LA_LINRX_ITHRESH_I ) = ITHRESH
402: ELSE
403: ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
404: END IF
405: END IF
406: IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
407: IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0D+0 ) THEN
408: IF ( IGNORE_CWISE ) THEN
409: PARAMS( LA_LINRX_CWISE_I ) = 0.0D+0
410: ELSE
411: PARAMS( LA_LINRX_CWISE_I ) = 1.0D+0
412: END IF
413: ELSE
414: IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0D+0
415: END IF
416: END IF
417: IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
418: N_NORMS = 0
419: ELSE IF ( IGNORE_CWISE ) THEN
420: N_NORMS = 1
421: ELSE
422: N_NORMS = 2
423: END IF
424: *
425: NOTRAN = LSAME( TRANS, 'N' )
426: ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
427: COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
428: *
429: * Test input parameters.
430: *
431: IF( TRANS_TYPE.EQ.-1 ) THEN
432: INFO = -1
433: ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND.
434: $ .NOT.LSAME( EQUED, 'N' ) ) THEN
435: INFO = -2
436: ELSE IF( N.LT.0 ) THEN
437: INFO = -3
438: ELSE IF( KL.LT.0 ) THEN
439: INFO = -4
440: ELSE IF( KU.LT.0 ) THEN
441: INFO = -5
442: ELSE IF( NRHS.LT.0 ) THEN
443: INFO = -6
444: ELSE IF( LDAB.LT.KL+KU+1 ) THEN
445: INFO = -8
446: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
447: INFO = -10
448: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
449: INFO = -13
450: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
451: INFO = -15
452: END IF
453: IF( INFO.NE.0 ) THEN
454: CALL XERBLA( 'DGBRFSX', -INFO )
455: RETURN
456: END IF
457: *
458: * Quick return if possible.
459: *
460: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
461: RCOND = 1.0D+0
462: DO J = 1, NRHS
463: BERR( J ) = 0.0D+0
464: IF ( N_ERR_BNDS .GE. 1 ) THEN
465: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
466: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
467: END IF
468: IF ( N_ERR_BNDS .GE. 2 ) THEN
469: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0D+0
470: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0D+0
471: END IF
472: IF ( N_ERR_BNDS .GE. 3 ) THEN
473: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0D+0
474: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0D+0
475: END IF
476: END DO
477: RETURN
478: END IF
479: *
480: * Default to failure.
481: *
482: RCOND = 0.0D+0
483: DO J = 1, NRHS
484: BERR( J ) = 1.0D+0
485: IF ( N_ERR_BNDS .GE. 1 ) THEN
486: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
487: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
488: END IF
489: IF ( N_ERR_BNDS .GE. 2 ) THEN
490: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
491: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
492: END IF
493: IF ( N_ERR_BNDS .GE. 3 ) THEN
494: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0D+0
495: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0D+0
496: END IF
497: END DO
498: *
499: * Compute the norm of A and the reciprocal of the condition
500: * number of A.
501: *
502: IF( NOTRAN ) THEN
503: NORM = 'I'
504: ELSE
505: NORM = '1'
506: END IF
507: ANORM = DLANGB( NORM, N, KL, KU, AB, LDAB, WORK )
508: CALL DGBCON( NORM, N, KL, KU, AFB, LDAFB, IPIV, ANORM, RCOND,
509: $ WORK, IWORK, INFO )
510: *
511: * Perform refinement on each right-hand side
512: *
513: IF (REF_TYPE .NE. 0) THEN
514:
515: PREC_TYPE = ILAPREC( 'E' )
516:
517: IF ( NOTRAN ) THEN
518: CALL DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
519: $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, COLEQU, C, B,
520: $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
521: $ ERR_BNDS_COMP, WORK( N+1 ), WORK( 1 ), WORK( 2*N+1 ),
522: $ WORK( 1 ), RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH,
523: $ IGNORE_CWISE, INFO )
524: ELSE
525: CALL DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU,
526: $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, ROWEQU, R, B,
527: $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
528: $ ERR_BNDS_COMP, WORK( N+1 ), WORK( 1 ), WORK( 2*N+1 ),
529: $ WORK( 1 ), RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH,
530: $ IGNORE_CWISE, INFO )
531: END IF
532: END IF
533:
534: ERR_LBND = MAX( 10.0D+0, SQRT( DBLE( N ) ) ) * DLAMCH( 'Epsilon' )
535: IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN
536: *
537: * Compute scaled normwise condition number cond(A*C).
538: *
539: IF ( COLEQU .AND. NOTRAN ) THEN
540: RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB,
541: $ LDAFB, IPIV, -1, C, INFO, WORK, IWORK )
542: ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN
543: RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB,
544: $ LDAFB, IPIV, -1, R, INFO, WORK, IWORK )
545: ELSE
546: RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB,
547: $ LDAFB, IPIV, 0, R, INFO, WORK, IWORK )
548: END IF
549: DO J = 1, NRHS
550: *
551: * Cap the error at 1.0.
552: *
553: IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
554: $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
555: $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
556: *
557: * Threshold the error (see LAWN).
558: *
559: IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
560: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0D+0
561: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0D+0
562: IF ( INFO .LE. N ) INFO = N + J
563: ELSE IF ( ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND )
564: $ THEN
565: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
566: ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0D+0
567: END IF
568: *
569: * Save the condition number.
570: *
571: IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
572: ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
573: END IF
574:
575: END DO
576: END IF
577:
578: IF (N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2) THEN
579: *
580: * Compute componentwise condition number cond(A*diag(Y(:,J))) for
581: * each right-hand side using the current solution as an estimate of
582: * the true solution. If the componentwise error estimate is too
583: * large, then the solution is a lousy estimate of truth and the
584: * estimated RCOND may be too optimistic. To avoid misleading users,
585: * the inverse condition number is set to 0.0 when the estimated
586: * cwise error is at least CWISE_WRONG.
587: *
588: CWISE_WRONG = SQRT( DLAMCH( 'Epsilon' ) )
589: DO J = 1, NRHS
590: IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
591: $ THEN
592: RCOND_TMP = DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, AFB,
593: $ LDAFB, IPIV, 1, X( 1, J ), INFO, WORK, IWORK )
594: ELSE
595: RCOND_TMP = 0.0D+0
596: END IF
597: *
598: * Cap the error at 1.0.
599: *
600: IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
601: $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0D+0 )
602: $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
603: *
604: * Threshold the error (see LAWN).
605: *
606: IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
607: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0D+0
608: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0D+0
609: IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0D+0
610: $ .AND. INFO.LT.N + J ) INFO = N + J
611: ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
612: $ .LT. ERR_LBND ) THEN
613: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
614: ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0D+0
615: END IF
616: *
617: * Save the condition number.
618: *
619: IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
620: ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
621: END IF
622:
623: END DO
624: END IF
625: *
626: RETURN
627: *
628: * End of DGBRFSX
629: *
630: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbrfsx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack