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