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