1: *> \brief \b DLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric indefinite matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.
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 DLA_SYRFSX_EXTENDED 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 December 2016
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.7.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: * December 2016
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, UPPER
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: INFO = 0
476: UPPER = LSAME( UPLO, 'U' )
477: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
478: INFO = -2
479: ELSE IF( N.LT.0 ) THEN
480: INFO = -3
481: ELSE IF( NRHS.LT.0 ) THEN
482: INFO = -4
483: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
484: INFO = -6
485: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
486: INFO = -8
487: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
488: INFO = -13
489: ELSE IF( LDY.LT.MAX( 1, N ) ) THEN
490: INFO = -15
491: END IF
492: IF( INFO.NE.0 ) THEN
493: CALL XERBLA( 'DLA_SYRFSX_EXTENDED', -INFO )
494: RETURN
495: END IF
496: EPS = DLAMCH( 'Epsilon' )
497: HUGEVAL = DLAMCH( 'Overflow' )
498: * Force HUGEVAL to Inf
499: HUGEVAL = HUGEVAL * HUGEVAL
500: * Using HUGEVAL may lead to spurious underflows.
501: INCR_THRESH = DBLE( N )*EPS
502:
503: IF ( LSAME ( UPLO, 'L' ) ) THEN
504: UPLO2 = ILAUPLO( 'L' )
505: ELSE
506: UPLO2 = ILAUPLO( 'U' )
507: ENDIF
508:
509: DO J = 1, NRHS
510: Y_PREC_STATE = EXTRA_RESIDUAL
511: IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN
512: DO I = 1, N
513: Y_TAIL( I ) = 0.0D+0
514: END DO
515: END IF
516:
517: DXRAT = 0.0D+0
518: DXRATMAX = 0.0D+0
519: DZRAT = 0.0D+0
520: DZRATMAX = 0.0D+0
521: FINAL_DX_X = HUGEVAL
522: FINAL_DZ_Z = HUGEVAL
523: PREVNORMDX = HUGEVAL
524: PREV_DZ_Z = HUGEVAL
525: DZ_Z = HUGEVAL
526: DX_X = HUGEVAL
527:
528: X_STATE = WORKING_STATE
529: Z_STATE = UNSTABLE_STATE
530: INCR_PREC = .FALSE.
531:
532: DO CNT = 1, ITHRESH
533: *
534: * Compute residual RES = B_s - op(A_s) * Y,
535: * op(A) = A, A**T, or A**H depending on TRANS (and type).
536: *
537: CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
538: IF (Y_PREC_STATE .EQ. BASE_RESIDUAL) THEN
539: CALL DSYMV( UPLO, N, -1.0D+0, A, LDA, Y(1,J), 1,
540: $ 1.0D+0, RES, 1 )
541: ELSE IF (Y_PREC_STATE .EQ. EXTRA_RESIDUAL) THEN
542: CALL BLAS_DSYMV_X( UPLO2, N, -1.0D+0, A, LDA,
543: $ Y( 1, J ), 1, 1.0D+0, RES, 1, PREC_TYPE )
544: ELSE
545: CALL BLAS_DSYMV2_X(UPLO2, N, -1.0D+0, A, LDA,
546: $ Y(1, J), Y_TAIL, 1, 1.0D+0, RES, 1, PREC_TYPE)
547: END IF
548:
549: ! XXX: RES is no longer needed.
550: CALL DCOPY( N, RES, 1, DY, 1 )
551: CALL DSYTRS( UPLO, N, 1, AF, LDAF, IPIV, DY, N, INFO )
552: *
553: * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
554: *
555: NORMX = 0.0D+0
556: NORMY = 0.0D+0
557: NORMDX = 0.0D+0
558: DZ_Z = 0.0D+0
559: YMIN = HUGEVAL
560:
561: DO I = 1, N
562: YK = ABS( Y( I, J ) )
563: DYK = ABS( DY( I ) )
564:
565: IF ( YK .NE. 0.0D+0 ) THEN
566: DZ_Z = MAX( DZ_Z, DYK / YK )
567: ELSE IF ( DYK .NE. 0.0D+0 ) THEN
568: DZ_Z = HUGEVAL
569: END IF
570:
571: YMIN = MIN( YMIN, YK )
572:
573: NORMY = MAX( NORMY, YK )
574:
575: IF ( COLEQU ) THEN
576: NORMX = MAX( NORMX, YK * C( I ) )
577: NORMDX = MAX( NORMDX, DYK * C( I ) )
578: ELSE
579: NORMX = NORMY
580: NORMDX = MAX(NORMDX, DYK)
581: END IF
582: END DO
583:
584: IF ( NORMX .NE. 0.0D+0 ) THEN
585: DX_X = NORMDX / NORMX
586: ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN
587: DX_X = 0.0D+0
588: ELSE
589: DX_X = HUGEVAL
590: END IF
591:
592: DXRAT = NORMDX / PREVNORMDX
593: DZRAT = DZ_Z / PREV_DZ_Z
594: *
595: * Check termination criteria.
596: *
597: IF ( YMIN*RCOND .LT. INCR_THRESH*NORMY
598: $ .AND. Y_PREC_STATE .LT. EXTRA_Y )
599: $ INCR_PREC = .TRUE.
600:
601: IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH )
602: $ X_STATE = WORKING_STATE
603: IF ( X_STATE .EQ. WORKING_STATE ) THEN
604: IF ( DX_X .LE. EPS ) THEN
605: X_STATE = CONV_STATE
606: ELSE IF ( DXRAT .GT. RTHRESH ) THEN
607: IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
608: INCR_PREC = .TRUE.
609: ELSE
610: X_STATE = NOPROG_STATE
611: END IF
612: ELSE
613: IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT
614: END IF
615: IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X
616: END IF
617:
618: IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB )
619: $ Z_STATE = WORKING_STATE
620: IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH )
621: $ Z_STATE = WORKING_STATE
622: IF ( Z_STATE .EQ. WORKING_STATE ) THEN
623: IF ( DZ_Z .LE. EPS ) THEN
624: Z_STATE = CONV_STATE
625: ELSE IF ( DZ_Z .GT. DZ_UB ) THEN
626: Z_STATE = UNSTABLE_STATE
627: DZRATMAX = 0.0D+0
628: FINAL_DZ_Z = HUGEVAL
629: ELSE IF ( DZRAT .GT. RTHRESH ) THEN
630: IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
631: INCR_PREC = .TRUE.
632: ELSE
633: Z_STATE = NOPROG_STATE
634: END IF
635: ELSE
636: IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT
637: END IF
638: IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
639: END IF
640:
641: IF ( X_STATE.NE.WORKING_STATE.AND.
642: $ ( IGNORE_CWISE.OR.Z_STATE.NE.WORKING_STATE ) )
643: $ GOTO 666
644:
645: IF ( INCR_PREC ) THEN
646: INCR_PREC = .FALSE.
647: Y_PREC_STATE = Y_PREC_STATE + 1
648: DO I = 1, N
649: Y_TAIL( I ) = 0.0D+0
650: END DO
651: END IF
652:
653: PREVNORMDX = NORMDX
654: PREV_DZ_Z = DZ_Z
655: *
656: * Update soluton.
657: *
658: IF (Y_PREC_STATE .LT. EXTRA_Y) THEN
659: CALL DAXPY( N, 1.0D+0, DY, 1, Y(1,J), 1 )
660: ELSE
661: CALL DLA_WWADDW( N, Y(1,J), Y_TAIL, DY )
662: END IF
663:
664: END DO
665: * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
666: 666 CONTINUE
667: *
668: * Set final_* when cnt hits ithresh.
669: *
670: IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X
671: IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
672: *
673: * Compute error bounds.
674: *
675: IF ( N_NORMS .GE. 1 ) THEN
676: ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) =
677: $ FINAL_DX_X / (1 - DXRATMAX)
678: END IF
679: IF ( N_NORMS .GE. 2 ) THEN
680: ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) =
681: $ FINAL_DZ_Z / (1 - DZRATMAX)
682: END IF
683: *
684: * Compute componentwise relative backward error from formula
685: * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
686: * where abs(Z) is the componentwise absolute value of the matrix
687: * or vector Z.
688: *
689: * Compute residual RES = B_s - op(A_s) * Y,
690: * op(A) = A, A**T, or A**H depending on TRANS (and type).
691: CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
692: CALL DSYMV( UPLO, N, -1.0D+0, A, LDA, Y(1,J), 1, 1.0D+0, RES,
693: $ 1 )
694:
695: DO I = 1, N
696: AYB( I ) = ABS( B( I, J ) )
697: END DO
698: *
699: * Compute abs(op(A_s))*abs(Y) + abs(B_s).
700: *
701: CALL DLA_SYAMV( UPLO2, N, 1.0D+0,
702: $ A, LDA, Y(1, J), 1, 1.0D+0, AYB, 1 )
703:
704: CALL DLA_LIN_BERR( N, N, 1, RES, AYB, BERR_OUT( J ) )
705: *
706: * End of loop for each RHS.
707: *
708: END DO
709: *
710: RETURN
711: 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