LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sla_syrfsx_extended.f
Go to the documentation of this file.
1 *> \brief \b SLA_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 SLA_SYRFSX_EXTENDED + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_syrfsx_extended.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_syrfsx_extended.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_syrfsx_extended.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLA_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 * REAL RTHRESH, DZ_UB
35 * ..
36 * .. Array Arguments ..
37 * INTEGER IPIV( * )
38 * REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
39 * $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
40 * REAL 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 *> SLA_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 SSYRFSX 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 REAL 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 REAL 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 SSYTRF.
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 SSYTRF.
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 REAL 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 REAL 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 REAL array, dimension (LDY,NRHS)
166 *> On entry, the solution matrix X, as computed by SSYTRS.
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 REAL 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 SLA_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 REAL 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 REAL 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 .LT. 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 REAL array, dimension (N)
296 *> Workspace to hold the intermediate residual.
297 *> \endverbatim
298 *>
299 *> \param[in] AYB
300 *> \verbatim
301 *> AYB is REAL 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 REAL array, dimension (N)
308 *> Workspace to hold the intermediate solution.
309 *> \endverbatim
310 *>
311 *> \param[in] Y_TAIL
312 *> \verbatim
313 *> Y_TAIL is REAL 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 REAL
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 REAL
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 REAL
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 definte 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 SLA_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 *> \date September 2012
387 *
388 *> \ingroup realSYcomputational
389 *
390 * =====================================================================
391  SUBROUTINE sla_syrfsx_extended( PREC_TYPE, UPLO, N, NRHS, A, LDA,
392  $ af, ldaf, ipiv, colequ, c, b, ldb,
393  $ y, ldy, berr_out, n_norms,
394  $ err_bnds_norm, err_bnds_comp, res,
395  $ ayb, dy, y_tail, rcond, ithresh,
396  $ rthresh, dz_ub, ignore_cwise,
397  $ info )
398 *
399 * -- LAPACK computational routine (version 3.4.2) --
400 * -- LAPACK is a software package provided by Univ. of Tennessee, --
401 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
402 * September 2012
403 *
404 * .. Scalar Arguments ..
405  INTEGER info, lda, ldaf, ldb, ldy, n, nrhs, prec_type,
406  $ n_norms, ithresh
407  CHARACTER uplo
408  LOGICAL colequ, ignore_cwise
409  REAL rthresh, dz_ub
410 * ..
411 * .. Array Arguments ..
412  INTEGER ipiv( * )
413  REAL a( lda, * ), af( ldaf, * ), b( ldb, * ),
414  $ y( ldy, * ), res( * ), dy( * ), y_tail( * )
415  REAL c( * ), ayb( * ), rcond, berr_out( * ),
416  $ err_bnds_norm( nrhs, * ),
417  $ err_bnds_comp( nrhs, * )
418 * ..
419 *
420 * =====================================================================
421 *
422 * .. Local Scalars ..
423  INTEGER uplo2, cnt, i, j, x_state, z_state
424  REAL yk, dyk, ymin, normy, normx, normdx, dxrat,
425  $ dzrat, prevnormdx, prev_dz_z, dxratmax,
426  $ dzratmax, dx_x, dz_z, final_dx_x, final_dz_z,
427  $ eps, hugeval, incr_thresh
428  LOGICAL incr_prec, upper
429 * ..
430 * .. Parameters ..
431  INTEGER unstable_state, working_state, conv_state,
432  $ noprog_state, y_prec_state, base_residual,
433  $ extra_residual, extra_y
434  parameter( unstable_state = 0, working_state = 1,
435  $ conv_state = 2, noprog_state = 3 )
436  parameter( base_residual = 0, extra_residual = 1,
437  $ extra_y = 2 )
438  INTEGER final_nrm_err_i, final_cmp_err_i, berr_i
439  INTEGER rcond_i, nrm_rcond_i, nrm_err_i, cmp_rcond_i
440  INTEGER cmp_err_i, piv_growth_i
441  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
442  $ berr_i = 3 )
443  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
444  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
445  $ piv_growth_i = 9 )
446  INTEGER la_linrx_itref_i, la_linrx_ithresh_i,
447  $ la_linrx_cwise_i
448  parameter( la_linrx_itref_i = 1,
449  $ la_linrx_ithresh_i = 2 )
450  parameter( la_linrx_cwise_i = 3 )
451  INTEGER la_linrx_trust_i, la_linrx_err_i,
452  $ la_linrx_rcond_i
453  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
454  parameter( la_linrx_rcond_i = 3 )
455 * ..
456 * .. External Functions ..
457  LOGICAL lsame
458  EXTERNAL ilauplo
459  INTEGER ilauplo
460 * ..
461 * .. External Subroutines ..
462  EXTERNAL saxpy, scopy, ssytrs, ssymv, blas_ssymv_x,
463  $ blas_ssymv2_x, sla_syamv, sla_wwaddw,
464  $ sla_lin_berr
465  REAL slamch
466 * ..
467 * .. Intrinsic Functions ..
468  INTRINSIC abs, max, min
469 * ..
470 * .. Executable Statements ..
471 *
472  info = 0
473  upper = lsame( uplo, 'U' )
474  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
475  info = -2
476  ELSE IF( n.LT.0 ) THEN
477  info = -3
478  ELSE IF( nrhs.LT.0 ) THEN
479  info = -4
480  ELSE IF( lda.LT.max( 1, n ) ) THEN
481  info = -6
482  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
483  info = -8
484  ELSE IF( ldb.LT.max( 1, n ) ) THEN
485  info = -13
486  ELSE IF( ldy.LT.max( 1, n ) ) THEN
487  info = -15
488  END IF
489  IF( info.NE.0 ) THEN
490  CALL xerbla( 'SLA_SYRFSX_EXTENDED', -info )
491  return
492  END IF
493  eps = slamch( 'Epsilon' )
494  hugeval = slamch( 'Overflow' )
495 * Force HUGEVAL to Inf
496  hugeval = hugeval * hugeval
497 * Using HUGEVAL may lead to spurious underflows.
498  incr_thresh = REAL( n )*eps
499 
500  IF ( lsame( uplo, 'L' ) ) THEN
501  uplo2 = ilauplo( 'L' )
502  ELSE
503  uplo2 = ilauplo( 'U' )
504  ENDIF
505 
506  DO j = 1, nrhs
507  y_prec_state = extra_residual
508  IF ( y_prec_state .EQ. extra_y ) THEN
509  DO i = 1, n
510  y_tail( i ) = 0.0
511  END DO
512  END IF
513 
514  dxrat = 0.0
515  dxratmax = 0.0
516  dzrat = 0.0
517  dzratmax = 0.0
518  final_dx_x = hugeval
519  final_dz_z = hugeval
520  prevnormdx = hugeval
521  prev_dz_z = hugeval
522  dz_z = hugeval
523  dx_x = hugeval
524 
525  x_state = working_state
526  z_state = unstable_state
527  incr_prec = .false.
528 
529  DO cnt = 1, ithresh
530 *
531 * Compute residual RES = B_s - op(A_s) * Y,
532 * op(A) = A, A**T, or A**H depending on TRANS (and type).
533 *
534  CALL scopy( n, b( 1, j ), 1, res, 1 )
535  IF (y_prec_state .EQ. base_residual) THEN
536  CALL ssymv( uplo, n, -1.0, a, lda, y(1,j), 1,
537  $ 1.0, res, 1 )
538  ELSE IF (y_prec_state .EQ. extra_residual) THEN
539  CALL blas_ssymv_x( uplo2, n, -1.0, a, lda,
540  $ y( 1, j ), 1, 1.0, res, 1, prec_type )
541  ELSE
542  CALL blas_ssymv2_x(uplo2, n, -1.0, a, lda,
543  $ y(1, j), y_tail, 1, 1.0, res, 1, prec_type)
544  END IF
545 
546 ! XXX: RES is no longer needed.
547  CALL scopy( n, res, 1, dy, 1 )
548  CALL ssytrs( uplo, n, 1, af, ldaf, ipiv, dy, n, info )
549 *
550 * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
551 *
552  normx = 0.0
553  normy = 0.0
554  normdx = 0.0
555  dz_z = 0.0
556  ymin = hugeval
557 
558  DO i = 1, n
559  yk = abs( y( i, j ) )
560  dyk = abs( dy( i ) )
561 
562  IF ( yk .NE. 0.0 ) THEN
563  dz_z = max( dz_z, dyk / yk )
564  ELSE IF ( dyk .NE. 0.0 ) THEN
565  dz_z = hugeval
566  END IF
567 
568  ymin = min( ymin, yk )
569 
570  normy = max( normy, yk )
571 
572  IF ( colequ ) THEN
573  normx = max( normx, yk * c( i ) )
574  normdx = max( normdx, dyk * c( i ) )
575  ELSE
576  normx = normy
577  normdx = max(normdx, dyk)
578  END IF
579  END DO
580 
581  IF ( normx .NE. 0.0 ) THEN
582  dx_x = normdx / normx
583  ELSE IF ( normdx .EQ. 0.0 ) THEN
584  dx_x = 0.0
585  ELSE
586  dx_x = hugeval
587  END IF
588 
589  dxrat = normdx / prevnormdx
590  dzrat = dz_z / prev_dz_z
591 *
592 * Check termination criteria.
593 *
594  IF ( ymin*rcond .LT. incr_thresh*normy
595  $ .AND. y_prec_state .LT. extra_y )
596  $ incr_prec = .true.
597 
598  IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
599  $ x_state = working_state
600  IF ( x_state .EQ. working_state ) THEN
601  IF ( dx_x .LE. eps ) THEN
602  x_state = conv_state
603  ELSE IF ( dxrat .GT. rthresh ) THEN
604  IF ( y_prec_state .NE. extra_y ) THEN
605  incr_prec = .true.
606  ELSE
607  x_state = noprog_state
608  END IF
609  ELSE
610  IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
611  END IF
612  IF ( x_state .GT. working_state ) final_dx_x = dx_x
613  END IF
614 
615  IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
616  $ z_state = working_state
617  IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
618  $ z_state = working_state
619  IF ( z_state .EQ. working_state ) THEN
620  IF ( dz_z .LE. eps ) THEN
621  z_state = conv_state
622  ELSE IF ( dz_z .GT. dz_ub ) THEN
623  z_state = unstable_state
624  dzratmax = 0.0
625  final_dz_z = hugeval
626  ELSE IF ( dzrat .GT. rthresh ) THEN
627  IF ( y_prec_state .NE. extra_y ) THEN
628  incr_prec = .true.
629  ELSE
630  z_state = noprog_state
631  END IF
632  ELSE
633  IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
634  END IF
635  IF ( z_state .GT. working_state ) final_dz_z = dz_z
636  END IF
637 
638  IF ( x_state.NE.working_state.AND.
639  $ ( ignore_cwise.OR.z_state.NE.working_state ) )
640  $ goto 666
641 
642  IF ( incr_prec ) THEN
643  incr_prec = .false.
644  y_prec_state = y_prec_state + 1
645  DO i = 1, n
646  y_tail( i ) = 0.0
647  END DO
648  END IF
649 
650  prevnormdx = normdx
651  prev_dz_z = dz_z
652 *
653 * Update soluton.
654 *
655  IF (y_prec_state .LT. extra_y) THEN
656  CALL saxpy( n, 1.0, dy, 1, y(1,j), 1 )
657  ELSE
658  CALL sla_wwaddw( n, y(1,j), y_tail, dy )
659  END IF
660 
661  END DO
662 * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
663  666 continue
664 *
665 * Set final_* when cnt hits ithresh.
666 *
667  IF ( x_state .EQ. working_state ) final_dx_x = dx_x
668  IF ( z_state .EQ. working_state ) final_dz_z = dz_z
669 *
670 * Compute error bounds.
671 *
672  IF ( n_norms .GE. 1 ) THEN
673  err_bnds_norm( j, la_linrx_err_i ) =
674  $ final_dx_x / (1 - dxratmax)
675  END IF
676  IF ( n_norms .GE. 2 ) THEN
677  err_bnds_comp( j, la_linrx_err_i ) =
678  $ final_dz_z / (1 - dzratmax)
679  END IF
680 *
681 * Compute componentwise relative backward error from formula
682 * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
683 * where abs(Z) is the componentwise absolute value of the matrix
684 * or vector Z.
685 *
686 * Compute residual RES = B_s - op(A_s) * Y,
687 * op(A) = A, A**T, or A**H depending on TRANS (and type).
688  CALL scopy( n, b( 1, j ), 1, res, 1 )
689  CALL ssymv( uplo, n, -1.0, a, lda, y(1,j), 1, 1.0, res, 1 )
690 
691  DO i = 1, n
692  ayb( i ) = abs( b( i, j ) )
693  END DO
694 *
695 * Compute abs(op(A_s))*abs(Y) + abs(B_s).
696 *
697  CALL sla_syamv( uplo2, n, 1.0,
698  $ a, lda, y(1, j), 1, 1.0, ayb, 1 )
699 
700  CALL sla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
701 *
702 * End of loop for each RHS.
703 *
704  END DO
705 *
706  return
707  END