1: SUBROUTINE ZLA_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
21: LOGICAL COLEQU, IGNORE_CWISE
22: INTEGER ITHRESH
23: DOUBLE PRECISION RTHRESH, DZ_UB
24: * ..
25: * .. Array Arguments
26: INTEGER IPIV( * )
27: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
28: $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
29: DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
30: $ ERRS_N( NRHS, * ), ERRS_C( NRHS, * )
31: * ..
32: *
33: * Purpose
34: * =======
35: *
36: * ZLA_GERFSX_EXTENDED improves the computed solution to a system of
37: * linear equations by performing extra-precise iterative refinement
38: * and provides error bounds and backward error estimates for the solution.
39: * This subroutine is called by ZGERFSX to perform iterative refinement.
40: * In addition to normwise error bound, the code provides maximum
41: * componentwise error bound if possible. See comments for ERR_BNDS_NORM
42: * and ERR_BNDS_COMP for details of the error bounds. Note that this
43: * subroutine is only resonsible for setting the second fields of
44: * ERR_BNDS_NORM and ERR_BNDS_COMP.
45: *
46: * Arguments
47: * =========
48: *
49: * PREC_TYPE (input) INTEGER
50: * Specifies the intermediate precision to be used in refinement.
51: * The value is defined by ILAPREC(P) where P is a CHARACTER and
52: * P = 'S': Single
53: * = 'D': Double
54: * = 'I': Indigenous
55: * = 'X', 'E': Extra
56: *
57: * TRANS_TYPE (input) INTEGER
58: * Specifies the transposition operation on A.
59: * The value is defined by ILATRANS(T) where T is a CHARACTER and
60: * T = 'N': No transpose
61: * = 'T': Transpose
62: * = 'C': Conjugate transpose
63: *
64: * N (input) INTEGER
65: * The number of linear equations, i.e., the order of the
66: * matrix A. N >= 0.
67: *
68: * NRHS (input) INTEGER
69: * The number of right-hand-sides, i.e., the number of columns of the
70: * matrix B.
71: *
72: * A (input) COMPLEX*16 array, dimension (LDA,N)
73: * On entry, the N-by-N matrix A.
74: *
75: * LDA (input) INTEGER
76: * The leading dimension of the array A. LDA >= max(1,N).
77: *
78: * AF (input) COMPLEX*16 array, dimension (LDAF,N)
79: * The factors L and U from the factorization
80: * A = P*L*U as computed by ZGETRF.
81: *
82: * LDAF (input) INTEGER
83: * The leading dimension of the array AF. LDAF >= max(1,N).
84: *
85: * IPIV (input) INTEGER array, dimension (N)
86: * The pivot indices from the factorization A = P*L*U
87: * as computed by ZGETRF; row i of the matrix was interchanged
88: * with row IPIV(i).
89: *
90: * COLEQU (input) LOGICAL
91: * If .TRUE. then column equilibration was done to A before calling
92: * this routine. This is needed to compute the solution and error
93: * bounds correctly.
94: *
95: * C (input) DOUBLE PRECISION array, dimension (N)
96: * The column scale factors for A. If COLEQU = .FALSE., C
97: * is not accessed. If C is input, each element of C should be a power
98: * of the radix to ensure a reliable solution and error estimates.
99: * Scaling by powers of the radix does not cause rounding errors unless
100: * the result underflows or overflows. Rounding errors during scaling
101: * lead to refining with a matrix that is not equivalent to the
102: * input matrix, producing error estimates that may not be
103: * reliable.
104: *
105: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
106: * The right-hand-side matrix B.
107: *
108: * LDB (input) INTEGER
109: * The leading dimension of the array B. LDB >= max(1,N).
110: *
111: * Y (input/output) COMPLEX*16 array, dimension (LDY,NRHS)
112: * On entry, the solution matrix X, as computed by ZGETRS.
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 ZLA_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) COMPLEX*16 array, dimension (N)
226: * Workspace to hold the intermediate residual.
227: *
228: * AYB (input) DOUBLE PRECISION array, dimension (N)
229: * Workspace.
230: *
231: * DY (input) COMPLEX*16 array, dimension (N)
232: * Workspace to hold the intermediate solution.
233: *
234: * Y_TAIL (input) COMPLEX*16 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 ZGETRS 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: COMPLEX*16 ZDUM
291: * ..
292: * .. Parameters ..
293: INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
294: $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
295: $ EXTRA_Y
296: PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1,
297: $ CONV_STATE = 2,
298: $ NOPROG_STATE = 3 )
299: PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1,
300: $ EXTRA_Y = 2 )
301: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
302: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
303: INTEGER CMP_ERR_I, PIV_GROWTH_I
304: PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
305: $ BERR_I = 3 )
306: PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
307: PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
308: $ PIV_GROWTH_I = 9 )
309: INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
310: $ LA_LINRX_CWISE_I
311: PARAMETER ( LA_LINRX_ITREF_I = 1,
312: $ LA_LINRX_ITHRESH_I = 2 )
313: PARAMETER ( LA_LINRX_CWISE_I = 3 )
314: INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
315: $ LA_LINRX_RCOND_I
316: PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
317: PARAMETER ( LA_LINRX_RCOND_I = 3 )
318: * ..
319: * .. External Subroutines ..
320: EXTERNAL ZAXPY, ZCOPY, ZGETRS, ZGEMV, BLAS_ZGEMV_X,
321: $ BLAS_ZGEMV2_X, ZLA_GEAMV, ZLA_WWADDW, DLAMCH,
322: $ CHLA_TRANSTYPE, ZLA_LIN_BERR
323: DOUBLE PRECISION DLAMCH
324: CHARACTER CHLA_TRANSTYPE
325: * ..
326: * .. Intrinsic Functions ..
327: INTRINSIC ABS, MAX, MIN
328: * ..
329: * .. Statement Functions ..
330: DOUBLE PRECISION CABS1
331: * ..
332: * .. Statement Function Definitions ..
333: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
334: * ..
335: * .. Executable Statements ..
336: *
337: IF ( INFO.NE.0 ) RETURN
338: TRANS = CHLA_TRANSTYPE(TRANS_TYPE)
339: EPS = DLAMCH( 'Epsilon' )
340: HUGEVAL = DLAMCH( 'Overflow' )
341: * Force HUGEVAL to Inf
342: HUGEVAL = HUGEVAL * HUGEVAL
343: * Using HUGEVAL may lead to spurious underflows.
344: INCR_THRESH = DBLE( N ) * EPS
345: *
346: DO J = 1, NRHS
347: Y_PREC_STATE = EXTRA_RESIDUAL
348: IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN
349: DO I = 1, N
350: Y_TAIL( I ) = 0.0D+0
351: END DO
352: END IF
353:
354: DXRAT = 0.0D+0
355: DXRATMAX = 0.0D+0
356: DZRAT = 0.0D+0
357: DZRATMAX = 0.0D+0
358: FINAL_DX_X = HUGEVAL
359: FINAL_DZ_Z = HUGEVAL
360: PREVNORMDX = HUGEVAL
361: PREV_DZ_Z = HUGEVAL
362: DZ_Z = HUGEVAL
363: DX_X = HUGEVAL
364:
365: X_STATE = WORKING_STATE
366: Z_STATE = UNSTABLE_STATE
367: INCR_PREC = .FALSE.
368:
369: DO CNT = 1, ITHRESH
370: *
371: * Compute residual RES = B_s - op(A_s) * Y,
372: * op(A) = A, A**T, or A**H depending on TRANS (and type).
373: *
374: CALL ZCOPY( N, B( 1, J ), 1, RES, 1 )
375: IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN
376: CALL ZGEMV( TRANS, N, N, (-1.0D+0,0.0D+0), A, LDA,
377: $ Y( 1, J ), 1, (1.0D+0,0.0D+0), RES, 1)
378: ELSE IF (Y_PREC_STATE .EQ. EXTRA_RESIDUAL) THEN
379: CALL BLAS_ZGEMV_X( TRANS_TYPE, N, N, (-1.0D+0,0.0D+0), A,
380: $ LDA, Y( 1, J ), 1, (1.0D+0,0.0D+0),
381: $ RES, 1, PREC_TYPE )
382: ELSE
383: CALL BLAS_ZGEMV2_X( TRANS_TYPE, N, N, (-1.0D+0,0.0D+0),
384: $ A, LDA, Y(1, J), Y_TAIL, 1, (1.0D+0,0.0D+0), RES, 1,
385: $ PREC_TYPE)
386: END IF
387:
388: ! XXX: RES is no longer needed.
389: CALL ZCOPY( N, RES, 1, DY, 1 )
390: CALL ZGETRS( TRANS, N, 1, AF, LDAF, IPIV, DY, N, INFO )
391: *
392: * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
393: *
394: NORMX = 0.0D+0
395: NORMY = 0.0D+0
396: NORMDX = 0.0D+0
397: DZ_Z = 0.0D+0
398: YMIN = HUGEVAL
399: *
400: DO I = 1, N
401: YK = CABS1( Y( I, J ) )
402: DYK = CABS1( DY( I ) )
403:
404: IF ( YK .NE. 0.0D+0 ) THEN
405: DZ_Z = MAX( DZ_Z, DYK / YK )
406: ELSE IF ( DYK .NE. 0.0D+0 ) THEN
407: DZ_Z = HUGEVAL
408: END IF
409:
410: YMIN = MIN( YMIN, YK )
411:
412: NORMY = MAX( NORMY, YK )
413:
414: IF ( COLEQU ) THEN
415: NORMX = MAX( NORMX, YK * C( I ) )
416: NORMDX = MAX( NORMDX, DYK * C( I ) )
417: ELSE
418: NORMX = NORMY
419: NORMDX = MAX(NORMDX, DYK)
420: END IF
421: END DO
422:
423: IF ( NORMX .NE. 0.0D+0 ) THEN
424: DX_X = NORMDX / NORMX
425: ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN
426: DX_X = 0.0D+0
427: ELSE
428: DX_X = HUGEVAL
429: END IF
430:
431: DXRAT = NORMDX / PREVNORMDX
432: DZRAT = DZ_Z / PREV_DZ_Z
433: *
434: * Check termination criteria
435: *
436: IF (.NOT.IGNORE_CWISE
437: $ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY
438: $ .AND. Y_PREC_STATE .LT. EXTRA_Y )
439: $ INCR_PREC = .TRUE.
440:
441: IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH )
442: $ X_STATE = WORKING_STATE
443: IF ( X_STATE .EQ. WORKING_STATE ) THEN
444: IF (DX_X .LE. EPS) THEN
445: X_STATE = CONV_STATE
446: ELSE IF ( DXRAT .GT. RTHRESH ) THEN
447: IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
448: INCR_PREC = .TRUE.
449: ELSE
450: X_STATE = NOPROG_STATE
451: END IF
452: ELSE
453: IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT
454: END IF
455: IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X
456: END IF
457:
458: IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB )
459: $ Z_STATE = WORKING_STATE
460: IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH )
461: $ Z_STATE = WORKING_STATE
462: IF ( Z_STATE .EQ. WORKING_STATE ) THEN
463: IF ( DZ_Z .LE. EPS ) THEN
464: Z_STATE = CONV_STATE
465: ELSE IF ( DZ_Z .GT. DZ_UB ) THEN
466: Z_STATE = UNSTABLE_STATE
467: DZRATMAX = 0.0D+0
468: FINAL_DZ_Z = HUGEVAL
469: ELSE IF ( DZRAT .GT. RTHRESH ) THEN
470: IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
471: INCR_PREC = .TRUE.
472: ELSE
473: Z_STATE = NOPROG_STATE
474: END IF
475: ELSE
476: IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT
477: END IF
478: IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
479: END IF
480: *
481: * Exit if both normwise and componentwise stopped working,
482: * but if componentwise is unstable, let it go at least two
483: * iterations.
484: *
485: IF ( X_STATE.NE.WORKING_STATE ) THEN
486: IF ( IGNORE_CWISE ) GOTO 666
487: IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE )
488: $ GOTO 666
489: IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666
490: END IF
491:
492: IF ( INCR_PREC ) THEN
493: INCR_PREC = .FALSE.
494: Y_PREC_STATE = Y_PREC_STATE + 1
495: DO I = 1, N
496: Y_TAIL( I ) = 0.0D+0
497: END DO
498: END IF
499:
500: PREVNORMDX = NORMDX
501: PREV_DZ_Z = DZ_Z
502: *
503: * Update soluton.
504: *
505: IF ( Y_PREC_STATE .LT. EXTRA_Y ) THEN
506: CALL ZAXPY( N, (1.0D+0,0.0D+0), DY, 1, Y(1,J), 1 )
507: ELSE
508: CALL ZLA_WWADDW( N, Y( 1, J ), Y_TAIL, DY )
509: END IF
510:
511: END DO
512: * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
513: 666 CONTINUE
514: *
515: * Set final_* when cnt hits ithresh
516: *
517: IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X
518: IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
519: *
520: * Compute error bounds
521: *
522: IF (N_NORMS .GE. 1) THEN
523: ERRS_N( J, LA_LINRX_ERR_I ) = FINAL_DX_X / (1 - DXRATMAX)
524:
525: END IF
526: IF ( N_NORMS .GE. 2 ) THEN
527: ERRS_C( J, LA_LINRX_ERR_I ) = FINAL_DZ_Z / (1 - DZRATMAX)
528: END IF
529: *
530: * Compute componentwise relative backward error from formula
531: * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
532: * where abs(Z) is the componentwise absolute value of the matrix
533: * or vector Z.
534: *
535: * Compute residual RES = B_s - op(A_s) * Y,
536: * op(A) = A, A**T, or A**H depending on TRANS (and type).
537: *
538: CALL ZCOPY( N, B( 1, J ), 1, RES, 1 )
539: CALL ZGEMV( TRANS, N, N, (-1.0D+0,0.0D+0), A, LDA, Y(1,J), 1,
540: $ (1.0D+0,0.0D+0), RES, 1 )
541:
542: DO I = 1, N
543: AYB( I ) = CABS1( B( I, J ) )
544: END DO
545: *
546: * Compute abs(op(A_s))*abs(Y) + abs(B_s).
547: *
548: CALL ZLA_GEAMV ( TRANS_TYPE, N, N, 1.0D+0,
549: $ A, LDA, Y(1, J), 1, 1.0D+0, AYB, 1 )
550:
551: CALL ZLA_LIN_BERR ( N, N, 1, RES, AYB, BERR_OUT( J ) )
552: *
553: * End of loop for each RHS.
554: *
555: END DO
556: *
557: RETURN
558: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_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