Annotation of rpl/lapack/lapack/dla_gerfsx_extended.f, revision 1.2
1.1 bertrand 1: SUBROUTINE DLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, NRHS, A,
2: $ LDA, AF, LDAF, IPIV, COLEQU, C, B,
3: $ LDB, Y, LDY, BERR_OUT, N_NORMS,
4: $ ERRS_N, ERRS_C, RES, AYB, DY,
5: $ Y_TAIL, RCOND, ITHRESH, RTHRESH,
6: $ DZ_UB, IGNORE_CWISE, INFO )
7: *
8: * -- LAPACK routine (version 3.2.1) --
9: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
10: * -- Jason Riedy of Univ. of California Berkeley. --
11: * -- April 2009 --
12: *
13: * -- LAPACK is a software package provided by Univ. of Tennessee, --
14: * -- Univ. of California Berkeley and NAG Ltd. --
15: *
16: IMPLICIT NONE
17: * ..
18: * .. Scalar Arguments ..
19: INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
20: $ TRANS_TYPE, N_NORMS, ITHRESH
21: LOGICAL COLEQU, IGNORE_CWISE
22: DOUBLE PRECISION RTHRESH, DZ_UB
23: * ..
24: * .. Array Arguments ..
25: INTEGER IPIV( * )
26: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
27: $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
28: DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
29: $ ERRS_N( NRHS, * ), ERRS_C( NRHS, * )
30: * ..
31: *
32: * Purpose
33: * =======
34: *
35: * DLA_GERFSX_EXTENDED improves the computed solution to a system of
36: * linear equations by performing extra-precise iterative refinement
37: * and provides error bounds and backward error estimates for the solution.
38: * This subroutine is called by DGERFSX to perform iterative refinement.
39: * In addition to normwise error bound, the code provides maximum
40: * componentwise error bound if possible. See comments for ERR_BNDS_NORM
41: * and ERR_BNDS_COMP for details of the error bounds. Note that this
42: * subroutine is only resonsible for setting the second fields of
43: * ERR_BNDS_NORM and ERR_BNDS_COMP.
44: *
45: * Arguments
46: * =========
47: *
48: * PREC_TYPE (input) INTEGER
49: * Specifies the intermediate precision to be used in refinement.
50: * The value is defined by ILAPREC(P) where P is a CHARACTER and
51: * P = 'S': Single
52: * = 'D': Double
53: * = 'I': Indigenous
54: * = 'X', 'E': Extra
55: *
56: * TRANS_TYPE (input) INTEGER
57: * Specifies the transposition operation on A.
58: * The value is defined by ILATRANS(T) where T is a CHARACTER and
59: * T = 'N': No transpose
60: * = 'T': Transpose
61: * = 'C': Conjugate transpose
62: *
63: * N (input) INTEGER
64: * The number of linear equations, i.e., the order of the
65: * matrix A. N >= 0.
66: *
67: * NRHS (input) INTEGER
68: * The number of right-hand-sides, i.e., the number of columns of the
69: * matrix B.
70: *
71: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
72: * On entry, the N-by-N matrix A.
73: *
74: * LDA (input) INTEGER
75: * The leading dimension of the array A. LDA >= max(1,N).
76: *
77: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
78: * The factors L and U from the factorization
79: * A = P*L*U as computed by DGETRF.
80: *
81: * LDAF (input) INTEGER
82: * The leading dimension of the array AF. LDAF >= max(1,N).
83: *
84: * IPIV (input) INTEGER array, dimension (N)
85: * The pivot indices from the factorization A = P*L*U
86: * as computed by DGETRF; row i of the matrix was interchanged
87: * with row IPIV(i).
88: *
89: * COLEQU (input) LOGICAL
90: * If .TRUE. then column equilibration was done to A before calling
91: * this routine. This is needed to compute the solution and error
92: * bounds correctly.
93: *
94: * C (input) DOUBLE PRECISION array, dimension (N)
95: * The column scale factors for A. If COLEQU = .FALSE., C
96: * is not accessed. If C is input, each element of C should be a power
97: * of the radix to ensure a reliable solution and error estimates.
98: * Scaling by powers of the radix does not cause rounding errors unless
99: * the result underflows or overflows. Rounding errors during scaling
100: * lead to refining with a matrix that is not equivalent to the
101: * input matrix, producing error estimates that may not be
102: * reliable.
103: *
104: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
105: * The right-hand-side matrix B.
106: *
107: * LDB (input) INTEGER
108: * The leading dimension of the array B. LDB >= max(1,N).
109: *
110: * Y (input/output) DOUBLE PRECISION array, dimension
111: * (LDY,NRHS)
112: * On entry, the solution matrix X, as computed by DGETRS.
113: * On exit, the improved solution matrix Y.
114: *
115: * LDY (input) INTEGER
116: * The leading dimension of the array Y. LDY >= max(1,N).
117: *
118: * BERR_OUT (output) DOUBLE PRECISION array, dimension (NRHS)
119: * On exit, BERR_OUT(j) contains the componentwise relative backward
120: * error for right-hand-side j from the formula
121: * max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
122: * where abs(Z) is the componentwise absolute value of the matrix
123: * or vector Z. This is computed by DLA_LIN_BERR.
124: *
125: * N_NORMS (input) INTEGER
126: * Determines which error bounds to return (see ERR_BNDS_NORM
127: * and ERR_BNDS_COMP).
128: * If N_NORMS >= 1 return normwise error bounds.
129: * If N_NORMS >= 2 return componentwise error bounds.
130: *
131: * ERR_BNDS_NORM (input/output) DOUBLE PRECISION array, dimension
132: * (NRHS, N_ERR_BNDS)
133: * For each right-hand side, this array contains information about
134: * various error bounds and condition numbers corresponding to the
135: * normwise relative error, which is defined as follows:
136: *
137: * Normwise relative error in the ith solution vector:
138: * max_j (abs(XTRUE(j,i) - X(j,i)))
139: * ------------------------------
140: * max_j abs(X(j,i))
141: *
142: * The array is indexed by the type of error information as described
143: * below. There currently are up to three pieces of information
144: * returned.
145: *
146: * The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
147: * right-hand side.
148: *
149: * The second index in ERR_BNDS_NORM(:,err) contains the following
150: * three fields:
151: * err = 1 "Trust/don't trust" boolean. Trust the answer if the
152: * reciprocal condition number is less than the threshold
153: * sqrt(n) * slamch('Epsilon').
154: *
155: * err = 2 "Guaranteed" error bound: The estimated forward error,
156: * almost certainly within a factor of 10 of the true error
157: * so long as the next entry is greater than the threshold
158: * sqrt(n) * slamch('Epsilon'). This error bound should only
159: * be trusted if the previous boolean is true.
160: *
161: * err = 3 Reciprocal condition number: Estimated normwise
162: * reciprocal condition number. Compared with the threshold
163: * sqrt(n) * slamch('Epsilon') to determine if the error
164: * estimate is "guaranteed". These reciprocal condition
165: * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
166: * appropriately scaled matrix Z.
167: * Let Z = S*A, where S scales each row by a power of the
168: * radix so all absolute row sums of Z are approximately 1.
169: *
170: * This subroutine is only responsible for setting the second field
171: * above.
172: * See Lapack Working Note 165 for further details and extra
173: * cautions.
174: *
175: * ERR_BNDS_COMP (input/output) DOUBLE PRECISION array, dimension
176: * (NRHS, N_ERR_BNDS)
177: * For each right-hand side, this array contains information about
178: * various error bounds and condition numbers corresponding to the
179: * componentwise relative error, which is defined as follows:
180: *
181: * Componentwise relative error in the ith solution vector:
182: * abs(XTRUE(j,i) - X(j,i))
183: * max_j ----------------------
184: * abs(X(j,i))
185: *
186: * The array is indexed by the right-hand side i (on which the
187: * componentwise relative error depends), and the type of error
188: * information as described below. There currently are up to three
189: * pieces of information returned for each right-hand side. If
190: * componentwise accuracy is not requested (PARAMS(3) = 0.0), then
191: * ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
192: * the first (:,N_ERR_BNDS) entries are returned.
193: *
194: * The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
195: * right-hand side.
196: *
197: * The second index in ERR_BNDS_COMP(:,err) contains the following
198: * three fields:
199: * err = 1 "Trust/don't trust" boolean. Trust the answer if the
200: * reciprocal condition number is less than the threshold
201: * sqrt(n) * slamch('Epsilon').
202: *
203: * err = 2 "Guaranteed" error bound: The estimated forward error,
204: * almost certainly within a factor of 10 of the true error
205: * so long as the next entry is greater than the threshold
206: * sqrt(n) * slamch('Epsilon'). This error bound should only
207: * be trusted if the previous boolean is true.
208: *
209: * err = 3 Reciprocal condition number: Estimated componentwise
210: * reciprocal condition number. Compared with the threshold
211: * sqrt(n) * slamch('Epsilon') to determine if the error
212: * estimate is "guaranteed". These reciprocal condition
213: * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
214: * appropriately scaled matrix Z.
215: * Let Z = S*(A*diag(x)), where x is the solution for the
216: * current right-hand side and S scales each row of
217: * A*diag(x) by a power of the radix so all absolute row
218: * sums of Z are approximately 1.
219: *
220: * This subroutine is only responsible for setting the second field
221: * above.
222: * See Lapack Working Note 165 for further details and extra
223: * cautions.
224: *
225: * RES (input) DOUBLE PRECISION array, dimension (N)
226: * Workspace to hold the intermediate residual.
227: *
228: * AYB (input) DOUBLE PRECISION array, dimension (N)
229: * Workspace. This can be the same workspace passed for Y_TAIL.
230: *
231: * DY (input) DOUBLE PRECISION array, dimension (N)
232: * Workspace to hold the intermediate solution.
233: *
234: * Y_TAIL (input) DOUBLE PRECISION array, dimension (N)
235: * Workspace to hold the trailing bits of the intermediate solution.
236: *
237: * RCOND (input) DOUBLE PRECISION
238: * Reciprocal scaled condition number. This is an estimate of the
239: * reciprocal Skeel condition number of the matrix A after
240: * equilibration (if done). If this is less than the machine
241: * precision (in particular, if it is zero), the matrix is singular
242: * to working precision. Note that the error may still be small even
243: * if this number is very small and the matrix appears ill-
244: * conditioned.
245: *
246: * ITHRESH (input) INTEGER
247: * The maximum number of residual computations allowed for
248: * refinement. The default is 10. For 'aggressive' set to 100 to
249: * permit convergence using approximate factorizations or
250: * factorizations other than LU. If the factorization uses a
251: * technique other than Gaussian elimination, the guarantees in
252: * ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
253: *
254: * RTHRESH (input) DOUBLE PRECISION
255: * Determines when to stop refinement if the error estimate stops
256: * decreasing. Refinement will stop when the next solution no longer
257: * satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
258: * the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
259: * default value is 0.5. For 'aggressive' set to 0.9 to permit
260: * convergence on extremely ill-conditioned matrices. See LAWN 165
261: * for more details.
262: *
263: * DZ_UB (input) DOUBLE PRECISION
264: * Determines when to start considering componentwise convergence.
265: * Componentwise convergence is only considered after each component
266: * of the solution Y is stable, which we definte as the relative
267: * change in each component being less than DZ_UB. The default value
268: * is 0.25, requiring the first bit to be stable. See LAWN 165 for
269: * more details.
270: *
271: * IGNORE_CWISE (input) LOGICAL
272: * If .TRUE. then ignore componentwise convergence. Default value
273: * is .FALSE..
274: *
275: * INFO (output) INTEGER
276: * = 0: Successful exit.
277: * < 0: if INFO = -i, the ith argument to DGETRS had an illegal
278: * value
279: *
280: * =====================================================================
281: *
282: * .. Local Scalars ..
283: CHARACTER TRANS
284: INTEGER CNT, I, J, X_STATE, Z_STATE, Y_PREC_STATE
285: DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
286: $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
287: $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
288: $ EPS, HUGEVAL, INCR_THRESH
289: LOGICAL INCR_PREC
290: * ..
291: * .. Parameters ..
292: INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
293: $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
294: $ EXTRA_Y
295: PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1,
296: $ CONV_STATE = 2, NOPROG_STATE = 3 )
297: PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1,
298: $ EXTRA_Y = 2 )
299: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
300: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
301: INTEGER CMP_ERR_I, PIV_GROWTH_I
302: PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
303: $ BERR_I = 3 )
304: PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
305: PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
306: $ PIV_GROWTH_I = 9 )
307: INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
308: $ LA_LINRX_CWISE_I
309: PARAMETER ( LA_LINRX_ITREF_I = 1,
310: $ LA_LINRX_ITHRESH_I = 2 )
311: PARAMETER ( LA_LINRX_CWISE_I = 3 )
312: INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
313: $ LA_LINRX_RCOND_I
314: PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
315: PARAMETER ( LA_LINRX_RCOND_I = 3 )
316: * ..
317: * .. External Subroutines ..
318: EXTERNAL DAXPY, DCOPY, DGETRS, DGEMV, BLAS_DGEMV_X,
319: $ BLAS_DGEMV2_X, DLA_GEAMV, DLA_WWADDW, DLAMCH,
320: $ CHLA_TRANSTYPE, DLA_LIN_BERR
321: DOUBLE PRECISION DLAMCH
322: CHARACTER CHLA_TRANSTYPE
323: * ..
324: * .. Intrinsic Functions ..
325: INTRINSIC ABS, MAX, MIN
326: * ..
327: * .. Executable Statements ..
328: *
329: IF ( INFO.NE.0 ) RETURN
330: TRANS = CHLA_TRANSTYPE(TRANS_TYPE)
331: EPS = DLAMCH( 'Epsilon' )
332: HUGEVAL = DLAMCH( 'Overflow' )
333: * Force HUGEVAL to Inf
334: HUGEVAL = HUGEVAL * HUGEVAL
335: * Using HUGEVAL may lead to spurious underflows.
336: INCR_THRESH = DBLE( N ) * EPS
337: *
338: DO J = 1, NRHS
339: Y_PREC_STATE = EXTRA_RESIDUAL
340: IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN
341: DO I = 1, N
342: Y_TAIL( I ) = 0.0D+0
343: END DO
344: END IF
345:
346: DXRAT = 0.0D+0
347: DXRATMAX = 0.0D+0
348: DZRAT = 0.0D+0
349: DZRATMAX = 0.0D+0
350: FINAL_DX_X = HUGEVAL
351: FINAL_DZ_Z = HUGEVAL
352: PREVNORMDX = HUGEVAL
353: PREV_DZ_Z = HUGEVAL
354: DZ_Z = HUGEVAL
355: DX_X = HUGEVAL
356:
357: X_STATE = WORKING_STATE
358: Z_STATE = UNSTABLE_STATE
359: INCR_PREC = .FALSE.
360:
361: DO CNT = 1, ITHRESH
362: *
363: * Compute residual RES = B_s - op(A_s) * Y,
364: * op(A) = A, A**T, or A**H depending on TRANS (and type).
365: *
366: CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
367: IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN
368: CALL DGEMV( TRANS, N, N, -1.0D+0, A, LDA, Y( 1, J ), 1,
369: $ 1.0D+0, RES, 1 )
370: ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN
371: CALL BLAS_DGEMV_X( TRANS_TYPE, N, N, -1.0D+0, A, LDA,
372: $ Y( 1, J ), 1, 1.0D+0, RES, 1, PREC_TYPE )
373: ELSE
374: CALL BLAS_DGEMV2_X( TRANS_TYPE, N, N, -1.0D+0, A, LDA,
375: $ Y( 1, J ), Y_TAIL, 1, 1.0D+0, RES, 1, PREC_TYPE )
376: END IF
377:
378: ! XXX: RES is no longer needed.
379: CALL DCOPY( N, RES, 1, DY, 1 )
380: CALL DGETRS( TRANS, N, 1, AF, LDAF, IPIV, DY, N, INFO )
381: *
382: * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
383: *
384: NORMX = 0.0D+0
385: NORMY = 0.0D+0
386: NORMDX = 0.0D+0
387: DZ_Z = 0.0D+0
388: YMIN = HUGEVAL
389: *
390: DO I = 1, N
391: YK = ABS( Y( I, J ) )
392: DYK = ABS( DY( I ) )
393:
394: IF ( YK .NE. 0.0D+0 ) THEN
395: DZ_Z = MAX( DZ_Z, DYK / YK )
396: ELSE IF ( DYK .NE. 0.0D+0 ) THEN
397: DZ_Z = HUGEVAL
398: END IF
399:
400: YMIN = MIN( YMIN, YK )
401:
402: NORMY = MAX( NORMY, YK )
403:
404: IF ( COLEQU ) THEN
405: NORMX = MAX( NORMX, YK * C( I ) )
406: NORMDX = MAX( NORMDX, DYK * C( I ) )
407: ELSE
408: NORMX = NORMY
409: NORMDX = MAX( NORMDX, DYK )
410: END IF
411: END DO
412:
413: IF ( NORMX .NE. 0.0D+0 ) THEN
414: DX_X = NORMDX / NORMX
415: ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN
416: DX_X = 0.0D+0
417: ELSE
418: DX_X = HUGEVAL
419: END IF
420:
421: DXRAT = NORMDX / PREVNORMDX
422: DZRAT = DZ_Z / PREV_DZ_Z
423: *
424: * Check termination criteria
425: *
426: IF (.NOT.IGNORE_CWISE
427: $ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY
428: $ .AND. Y_PREC_STATE .LT. EXTRA_Y)
429: $ INCR_PREC = .TRUE.
430:
431: IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH )
432: $ X_STATE = WORKING_STATE
433: IF ( X_STATE .EQ. WORKING_STATE ) THEN
434: IF ( DX_X .LE. EPS ) THEN
435: X_STATE = CONV_STATE
436: ELSE IF ( DXRAT .GT. RTHRESH ) THEN
437: IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
438: INCR_PREC = .TRUE.
439: ELSE
440: X_STATE = NOPROG_STATE
441: END IF
442: ELSE
443: IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT
444: END IF
445: IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X
446: END IF
447:
448: IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB )
449: $ Z_STATE = WORKING_STATE
450: IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH )
451: $ Z_STATE = WORKING_STATE
452: IF ( Z_STATE .EQ. WORKING_STATE ) THEN
453: IF ( DZ_Z .LE. EPS ) THEN
454: Z_STATE = CONV_STATE
455: ELSE IF ( DZ_Z .GT. DZ_UB ) THEN
456: Z_STATE = UNSTABLE_STATE
457: DZRATMAX = 0.0D+0
458: FINAL_DZ_Z = HUGEVAL
459: ELSE IF ( DZRAT .GT. RTHRESH ) THEN
460: IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
461: INCR_PREC = .TRUE.
462: ELSE
463: Z_STATE = NOPROG_STATE
464: END IF
465: ELSE
466: IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT
467: END IF
468: IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
469: END IF
470: *
471: * Exit if both normwise and componentwise stopped working,
472: * but if componentwise is unstable, let it go at least two
473: * iterations.
474: *
475: IF ( X_STATE.NE.WORKING_STATE ) THEN
476: IF ( IGNORE_CWISE) GOTO 666
477: IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE )
478: $ GOTO 666
479: IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666
480: END IF
481:
482: IF ( INCR_PREC ) THEN
483: INCR_PREC = .FALSE.
484: Y_PREC_STATE = Y_PREC_STATE + 1
485: DO I = 1, N
486: Y_TAIL( I ) = 0.0D+0
487: END DO
488: END IF
489:
490: PREVNORMDX = NORMDX
491: PREV_DZ_Z = DZ_Z
492: *
493: * Update soluton.
494: *
495: IF ( Y_PREC_STATE .LT. EXTRA_Y ) THEN
496: CALL DAXPY( N, 1.0D+0, DY, 1, Y( 1, J ), 1 )
497: ELSE
498: CALL DLA_WWADDW( N, Y( 1, J ), Y_TAIL, DY )
499: END IF
500:
501: END DO
502: * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
503: 666 CONTINUE
504: *
505: * Set final_* when cnt hits ithresh.
506: *
507: IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X
508: IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
509: *
510: * Compute error bounds
511: *
512: IF (N_NORMS .GE. 1) THEN
513: ERRS_N( J, LA_LINRX_ERR_I ) = FINAL_DX_X / (1 - DXRATMAX)
514: END IF
515: IF ( N_NORMS .GE. 2 ) THEN
516: ERRS_C( J, LA_LINRX_ERR_I ) = FINAL_DZ_Z / (1 - DZRATMAX)
517: END IF
518: *
519: * Compute componentwise relative backward error from formula
520: * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
521: * where abs(Z) is the componentwise absolute value of the matrix
522: * or vector Z.
523: *
524: * Compute residual RES = B_s - op(A_s) * Y,
525: * op(A) = A, A**T, or A**H depending on TRANS (and type).
526: *
527: CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
528: CALL DGEMV( TRANS, N, N, -1.0D+0, A, LDA, Y(1,J), 1, 1.0D+0,
529: $ RES, 1 )
530:
531: DO I = 1, N
532: AYB( I ) = ABS( B( I, J ) )
533: END DO
534: *
535: * Compute abs(op(A_s))*abs(Y) + abs(B_s).
536: *
537: CALL DLA_GEAMV ( TRANS_TYPE, N, N, 1.0D+0,
538: $ A, LDA, Y(1, J), 1, 1.0D+0, AYB, 1 )
539:
540: CALL DLA_LIN_BERR ( N, N, 1, RES, AYB, BERR_OUT( J ) )
541: *
542: * End of loop for each RHS.
543: *
544: END DO
545: *
546: RETURN
547: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_gerfsx_extended.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack