LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dla_syrfsx_extended.f
Go to the documentation of this file.
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 responsible 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 P
71*> = 'S': Single
72*> = 'D': Double
73*> = 'I': Indigenous
74*> = 'X' or '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 (LDY,NRHS)
166*> On entry, the solution matrix X, as computed by DSYTRS.
167*> On exit, the improved solution matrix Y.
168*> \endverbatim
169*>
170*> \param[in] LDY
171*> \verbatim
172*> LDY is INTEGER
173*> The leading dimension of the array Y. LDY >= max(1,N).
174*> \endverbatim
175*>
176*> \param[out] BERR_OUT
177*> \verbatim
178*> BERR_OUT is DOUBLE PRECISION array, dimension (NRHS)
179*> On exit, BERR_OUT(j) contains the componentwise relative backward
180*> error for right-hand-side j from the formula
181*> max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
182*> where abs(Z) is the componentwise absolute value of the matrix
183*> or vector Z. This is computed by DLA_LIN_BERR.
184*> \endverbatim
185*>
186*> \param[in] N_NORMS
187*> \verbatim
188*> N_NORMS is INTEGER
189*> Determines which error bounds to return (see ERR_BNDS_NORM
190*> and ERR_BNDS_COMP).
191*> If N_NORMS >= 1 return normwise error bounds.
192*> If N_NORMS >= 2 return componentwise error bounds.
193*> \endverbatim
194*>
195*> \param[in,out] ERR_BNDS_NORM
196*> \verbatim
197*> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
198*> For each right-hand side, this array contains information about
199*> various error bounds and condition numbers corresponding to the
200*> normwise relative error, which is defined as follows:
201*>
202*> Normwise relative error in the ith solution vector:
203*> max_j (abs(XTRUE(j,i) - X(j,i)))
204*> ------------------------------
205*> max_j abs(X(j,i))
206*>
207*> The array is indexed by the type of error information as described
208*> below. There currently are up to three pieces of information
209*> returned.
210*>
211*> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
212*> right-hand side.
213*>
214*> The second index in ERR_BNDS_NORM(:,err) contains the following
215*> three fields:
216*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
217*> reciprocal condition number is less than the threshold
218*> sqrt(n) * slamch('Epsilon').
219*>
220*> err = 2 "Guaranteed" error bound: The estimated forward error,
221*> almost certainly within a factor of 10 of the true error
222*> so long as the next entry is greater than the threshold
223*> sqrt(n) * slamch('Epsilon'). This error bound should only
224*> be trusted if the previous boolean is true.
225*>
226*> err = 3 Reciprocal condition number: Estimated normwise
227*> reciprocal condition number. Compared with the threshold
228*> sqrt(n) * slamch('Epsilon') to determine if the error
229*> estimate is "guaranteed". These reciprocal condition
230*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
231*> appropriately scaled matrix Z.
232*> Let Z = S*A, where S scales each row by a power of the
233*> radix so all absolute row sums of Z are approximately 1.
234*>
235*> This subroutine is only responsible for setting the second field
236*> above.
237*> See Lapack Working Note 165 for further details and extra
238*> cautions.
239*> \endverbatim
240*>
241*> \param[in,out] ERR_BNDS_COMP
242*> \verbatim
243*> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
244*> For each right-hand side, this array contains information about
245*> various error bounds and condition numbers corresponding to the
246*> componentwise relative error, which is defined as follows:
247*>
248*> Componentwise relative error in the ith solution vector:
249*> abs(XTRUE(j,i) - X(j,i))
250*> max_j ----------------------
251*> abs(X(j,i))
252*>
253*> The array is indexed by the right-hand side i (on which the
254*> componentwise relative error depends), and the type of error
255*> information as described below. There currently are up to three
256*> pieces of information returned for each right-hand side. If
257*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
258*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
259*> the first (:,N_ERR_BNDS) entries are returned.
260*>
261*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
262*> right-hand side.
263*>
264*> The second index in ERR_BNDS_COMP(:,err) contains the following
265*> three fields:
266*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
267*> reciprocal condition number is less than the threshold
268*> sqrt(n) * slamch('Epsilon').
269*>
270*> err = 2 "Guaranteed" error bound: The estimated forward error,
271*> almost certainly within a factor of 10 of the true error
272*> so long as the next entry is greater than the threshold
273*> sqrt(n) * slamch('Epsilon'). This error bound should only
274*> be trusted if the previous boolean is true.
275*>
276*> err = 3 Reciprocal condition number: Estimated componentwise
277*> reciprocal condition number. Compared with the threshold
278*> sqrt(n) * slamch('Epsilon') to determine if the error
279*> estimate is "guaranteed". These reciprocal condition
280*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
281*> appropriately scaled matrix Z.
282*> Let Z = S*(A*diag(x)), where x is the solution for the
283*> current right-hand side and S scales each row of
284*> A*diag(x) by a power of the radix so all absolute row
285*> sums of Z are approximately 1.
286*>
287*> This subroutine is only responsible for setting the second field
288*> above.
289*> See Lapack Working Note 165 for further details and extra
290*> cautions.
291*> \endverbatim
292*>
293*> \param[in] RES
294*> \verbatim
295*> RES is DOUBLE PRECISION array, dimension (N)
296*> Workspace to hold the intermediate residual.
297*> \endverbatim
298*>
299*> \param[in] AYB
300*> \verbatim
301*> AYB is DOUBLE PRECISION array, dimension (N)
302*> Workspace. This can be the same workspace passed for Y_TAIL.
303*> \endverbatim
304*>
305*> \param[in] DY
306*> \verbatim
307*> DY is DOUBLE PRECISION array, dimension (N)
308*> Workspace to hold the intermediate solution.
309*> \endverbatim
310*>
311*> \param[in] Y_TAIL
312*> \verbatim
313*> Y_TAIL is DOUBLE PRECISION array, dimension (N)
314*> Workspace to hold the trailing bits of the intermediate solution.
315*> \endverbatim
316*>
317*> \param[in] RCOND
318*> \verbatim
319*> RCOND is DOUBLE PRECISION
320*> Reciprocal scaled condition number. This is an estimate of the
321*> reciprocal Skeel condition number of the matrix A after
322*> equilibration (if done). If this is less than the machine
323*> precision (in particular, if it is zero), the matrix is singular
324*> to working precision. Note that the error may still be small even
325*> if this number is very small and the matrix appears ill-
326*> conditioned.
327*> \endverbatim
328*>
329*> \param[in] ITHRESH
330*> \verbatim
331*> ITHRESH is INTEGER
332*> The maximum number of residual computations allowed for
333*> refinement. The default is 10. For 'aggressive' set to 100 to
334*> permit convergence using approximate factorizations or
335*> factorizations other than LU. If the factorization uses a
336*> technique other than Gaussian elimination, the guarantees in
337*> ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
338*> \endverbatim
339*>
340*> \param[in] RTHRESH
341*> \verbatim
342*> RTHRESH is DOUBLE PRECISION
343*> Determines when to stop refinement if the error estimate stops
344*> decreasing. Refinement will stop when the next solution no longer
345*> satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
346*> the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
347*> default value is 0.5. For 'aggressive' set to 0.9 to permit
348*> convergence on extremely ill-conditioned matrices. See LAWN 165
349*> for more details.
350*> \endverbatim
351*>
352*> \param[in] DZ_UB
353*> \verbatim
354*> DZ_UB is DOUBLE PRECISION
355*> Determines when to start considering componentwise convergence.
356*> Componentwise convergence is only considered after each component
357*> of the solution Y is stable, which we define as the relative
358*> change in each component being less than DZ_UB. The default value
359*> is 0.25, requiring the first bit to be stable. See LAWN 165 for
360*> more details.
361*> \endverbatim
362*>
363*> \param[in] IGNORE_CWISE
364*> \verbatim
365*> IGNORE_CWISE is LOGICAL
366*> If .TRUE. then ignore componentwise convergence. Default value
367*> is .FALSE..
368*> \endverbatim
369*>
370*> \param[out] INFO
371*> \verbatim
372*> INFO is INTEGER
373*> = 0: Successful exit.
374*> < 0: if INFO = -i, the ith argument to DLA_SYRFSX_EXTENDED had an illegal
375*> value
376*> \endverbatim
377*
378* Authors:
379* ========
380*
381*> \author Univ. of Tennessee
382*> \author Univ. of California Berkeley
383*> \author Univ. of Colorado Denver
384*> \author NAG Ltd.
385*
386*> \ingroup la_herfsx_extended
387*
388* =====================================================================
389 SUBROUTINE dla_syrfsx_extended( PREC_TYPE, UPLO, N, NRHS, A, LDA,
390 $ AF, LDAF, IPIV, COLEQU, C, B, LDB,
391 $ Y, LDY, BERR_OUT, N_NORMS,
392 $ ERR_BNDS_NORM, ERR_BNDS_COMP, RES,
393 $ AYB, DY, Y_TAIL, RCOND, ITHRESH,
394 $ RTHRESH, DZ_UB, IGNORE_CWISE,
395 $ INFO )
396*
397* -- LAPACK computational routine --
398* -- LAPACK is a software package provided by Univ. of Tennessee, --
399* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
400*
401* .. Scalar Arguments ..
402 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
403 $ N_NORMS, ITHRESH
404 CHARACTER UPLO
405 LOGICAL COLEQU, IGNORE_CWISE
406 DOUBLE PRECISION RTHRESH, DZ_UB
407* ..
408* .. Array Arguments ..
409 INTEGER IPIV( * )
410 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
411 $ y( ldy, * ), res( * ), dy( * ), y_tail( * )
412 DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
413 $ ERR_BNDS_NORM( NRHS, * ),
414 $ ERR_BNDS_COMP( NRHS, * )
415* ..
416*
417* =====================================================================
418*
419* .. Local Scalars ..
420 INTEGER UPLO2, CNT, I, J, X_STATE, Z_STATE
421 DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
422 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
423 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
424 $ EPS, HUGEVAL, INCR_THRESH
425 LOGICAL INCR_PREC, UPPER
426* ..
427* .. Parameters ..
428 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
429 $ NOPROG_STATE, Y_PREC_STATE, BASE_RESIDUAL,
430 $ EXTRA_RESIDUAL, EXTRA_Y
431 parameter( unstable_state = 0, working_state = 1,
432 $ conv_state = 2, noprog_state = 3 )
433 parameter( base_residual = 0, extra_residual = 1,
434 $ extra_y = 2 )
435 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
436 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
437 INTEGER CMP_ERR_I, PIV_GROWTH_I
438 PARAMETER ( FINAL_NRM_ERR_I = 1, final_cmp_err_i = 2,
439 $ berr_i = 3 )
440 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
441 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
442 $ piv_growth_i = 9 )
443 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
444 $ la_linrx_cwise_i
445 parameter( la_linrx_itref_i = 1,
446 $ la_linrx_ithresh_i = 2 )
447 parameter( la_linrx_cwise_i = 3 )
448 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
449 $ la_linrx_rcond_i
450 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
451 parameter( la_linrx_rcond_i = 3 )
452* ..
453* .. External Functions ..
454 LOGICAL LSAME
455 EXTERNAL ILAUPLO
456 INTEGER ILAUPLO
457* ..
458* .. External Subroutines ..
459 EXTERNAL daxpy, dcopy, dsytrs, dsymv, blas_dsymv_x,
460 $ blas_dsymv2_x, dla_syamv, dla_wwaddw,
462 DOUBLE PRECISION DLAMCH
463* ..
464* .. Intrinsic Functions ..
465 INTRINSIC abs, max, min
466* ..
467* .. Executable Statements ..
468*
469 info = 0
470 upper = lsame( uplo, 'U' )
471 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
472 info = -2
473 ELSE IF( n.LT.0 ) THEN
474 info = -3
475 ELSE IF( nrhs.LT.0 ) THEN
476 info = -4
477 ELSE IF( lda.LT.max( 1, n ) ) THEN
478 info = -6
479 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
480 info = -8
481 ELSE IF( ldb.LT.max( 1, n ) ) THEN
482 info = -13
483 ELSE IF( ldy.LT.max( 1, n ) ) THEN
484 info = -15
485 END IF
486 IF( info.NE.0 ) THEN
487 CALL xerbla( 'DLA_SYRFSX_EXTENDED', -info )
488 RETURN
489 END IF
490 eps = dlamch( 'Epsilon' )
491 hugeval = dlamch( 'Overflow' )
492* Force HUGEVAL to Inf
493 hugeval = hugeval * hugeval
494* Using HUGEVAL may lead to spurious underflows.
495 incr_thresh = dble( n )*eps
496
497 IF ( lsame( uplo, 'L' ) ) THEN
498 uplo2 = ilauplo( 'L' )
499 ELSE
500 uplo2 = ilauplo( 'U' )
501 ENDIF
502
503 DO j = 1, nrhs
504 y_prec_state = extra_residual
505 IF ( y_prec_state .EQ. extra_y ) THEN
506 DO i = 1, n
507 y_tail( i ) = 0.0d+0
508 END DO
509 END IF
510
511 dxrat = 0.0d+0
512 dxratmax = 0.0d+0
513 dzrat = 0.0d+0
514 dzratmax = 0.0d+0
515 final_dx_x = hugeval
516 final_dz_z = hugeval
517 prevnormdx = hugeval
518 prev_dz_z = hugeval
519 dz_z = hugeval
520 dx_x = hugeval
521
522 x_state = working_state
523 z_state = unstable_state
524 incr_prec = .false.
525
526 DO cnt = 1, ithresh
527*
528* Compute residual RES = B_s - op(A_s) * Y,
529* op(A) = A, A**T, or A**H depending on TRANS (and type).
530*
531 CALL dcopy( n, b( 1, j ), 1, res, 1 )
532 IF (y_prec_state .EQ. base_residual) THEN
533 CALL dsymv( uplo, n, -1.0d+0, a, lda, y(1,j), 1,
534 $ 1.0d+0, res, 1 )
535 ELSE IF (y_prec_state .EQ. extra_residual) THEN
536 CALL blas_dsymv_x( uplo2, n, -1.0d+0, a, lda,
537 $ y( 1, j ), 1, 1.0d+0, res, 1, prec_type )
538 ELSE
539 CALL blas_dsymv2_x(uplo2, n, -1.0d+0, a, lda,
540 $ y(1, j), y_tail, 1, 1.0d+0, res, 1, prec_type)
541 END IF
542
543! XXX: RES is no longer needed.
544 CALL dcopy( n, res, 1, dy, 1 )
545 CALL dsytrs( uplo, n, 1, af, ldaf, ipiv, dy, n, info )
546*
547* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
548*
549 normx = 0.0d+0
550 normy = 0.0d+0
551 normdx = 0.0d+0
552 dz_z = 0.0d+0
553 ymin = hugeval
554
555 DO i = 1, n
556 yk = abs( y( i, j ) )
557 dyk = abs( dy( i ) )
558
559 IF ( yk .NE. 0.0d+0 ) THEN
560 dz_z = max( dz_z, dyk / yk )
561 ELSE IF ( dyk .NE. 0.0d+0 ) THEN
562 dz_z = hugeval
563 END IF
564
565 ymin = min( ymin, yk )
566
567 normy = max( normy, yk )
568
569 IF ( colequ ) THEN
570 normx = max( normx, yk * c( i ) )
571 normdx = max( normdx, dyk * c( i ) )
572 ELSE
573 normx = normy
574 normdx = max(normdx, dyk)
575 END IF
576 END DO
577
578 IF ( normx .NE. 0.0d+0 ) THEN
579 dx_x = normdx / normx
580 ELSE IF ( normdx .EQ. 0.0d+0 ) THEN
581 dx_x = 0.0d+0
582 ELSE
583 dx_x = hugeval
584 END IF
585
586 dxrat = normdx / prevnormdx
587 dzrat = dz_z / prev_dz_z
588*
589* Check termination criteria.
590*
591 IF ( ymin*rcond .LT. incr_thresh*normy
592 $ .AND. y_prec_state .LT. extra_y )
593 $ incr_prec = .true.
594
595 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
596 $ x_state = working_state
597 IF ( x_state .EQ. working_state ) THEN
598 IF ( dx_x .LE. eps ) THEN
599 x_state = conv_state
600 ELSE IF ( dxrat .GT. rthresh ) THEN
601 IF ( y_prec_state .NE. extra_y ) THEN
602 incr_prec = .true.
603 ELSE
604 x_state = noprog_state
605 END IF
606 ELSE
607 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
608 END IF
609 IF ( x_state .GT. working_state ) final_dx_x = dx_x
610 END IF
611
612 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
613 $ z_state = working_state
614 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
615 $ z_state = working_state
616 IF ( z_state .EQ. working_state ) THEN
617 IF ( dz_z .LE. eps ) THEN
618 z_state = conv_state
619 ELSE IF ( dz_z .GT. dz_ub ) THEN
620 z_state = unstable_state
621 dzratmax = 0.0d+0
622 final_dz_z = hugeval
623 ELSE IF ( dzrat .GT. rthresh ) THEN
624 IF ( y_prec_state .NE. extra_y ) THEN
625 incr_prec = .true.
626 ELSE
627 z_state = noprog_state
628 END IF
629 ELSE
630 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
631 END IF
632 IF ( z_state .GT. working_state ) final_dz_z = dz_z
633 END IF
634
635 IF ( x_state.NE.working_state.AND.
636 $ ( ignore_cwise.OR.z_state.NE.working_state ) )
637 $ GOTO 666
638
639 IF ( incr_prec ) THEN
640 incr_prec = .false.
641 y_prec_state = y_prec_state + 1
642 DO i = 1, n
643 y_tail( i ) = 0.0d+0
644 END DO
645 END IF
646
647 prevnormdx = normdx
648 prev_dz_z = dz_z
649*
650* Update solution.
651*
652 IF (y_prec_state .LT. extra_y) THEN
653 CALL daxpy( n, 1.0d+0, dy, 1, y(1,j), 1 )
654 ELSE
655 CALL dla_wwaddw( n, y(1,j), y_tail, dy )
656 END IF
657
658 END DO
659* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
660 666 CONTINUE
661*
662* Set final_* when cnt hits ithresh.
663*
664 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
665 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
666*
667* Compute error bounds.
668*
669 IF ( n_norms .GE. 1 ) THEN
670 err_bnds_norm( j, la_linrx_err_i ) =
671 $ final_dx_x / (1 - dxratmax)
672 END IF
673 IF ( n_norms .GE. 2 ) THEN
674 err_bnds_comp( j, la_linrx_err_i ) =
675 $ final_dz_z / (1 - dzratmax)
676 END IF
677*
678* Compute componentwise relative backward error from formula
679* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
680* where abs(Z) is the componentwise absolute value of the matrix
681* or vector Z.
682*
683* Compute residual RES = B_s - op(A_s) * Y,
684* op(A) = A, A**T, or A**H depending on TRANS (and type).
685 CALL dcopy( n, b( 1, j ), 1, res, 1 )
686 CALL dsymv( uplo, n, -1.0d+0, a, lda, y(1,j), 1, 1.0d+0, res,
687 $ 1 )
688
689 DO i = 1, n
690 ayb( i ) = abs( b( i, j ) )
691 END DO
692*
693* Compute abs(op(A_s))*abs(Y) + abs(B_s).
694*
695 CALL dla_syamv( uplo2, n, 1.0d+0,
696 $ a, lda, y(1, j), 1, 1.0d+0, ayb, 1 )
697
698 CALL dla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
699*
700* End of loop for each RHS.
701*
702 END DO
703*
704 RETURN
705*
706* End of DLA_SYRFSX_EXTENDED
707*
708 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
subroutine dsytrs(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
DSYTRS
Definition dsytrs.f:120
subroutine dla_syamv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bou...
Definition dla_syamv.f:177
subroutine dla_syrfsx_extended(prec_type, uplo, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)
DLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric inde...
subroutine dla_lin_berr(n, nz, nrhs, res, ayb, berr)
DLA_LIN_BERR computes a component-wise relative backward error.
subroutine dla_wwaddw(n, x, y, w)
DLA_WWADDW adds a vector into a doubled-single vector.
Definition dla_wwaddw.f:81