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