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