LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 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 September 2012
390 *
391 *> \ingroup doubleSYcomputational
392 *
393 * =====================================================================
394  SUBROUTINE dla_syrfsx_extended( PREC_TYPE, UPLO, N, NRHS, A, LDA,
395  $ af, ldaf, ipiv, colequ, c, b, ldb,
396  $ y, ldy, berr_out, n_norms,
397  $ err_bnds_norm, err_bnds_comp, res,
398  $ ayb, dy, y_tail, rcond, ithresh,
399  $ rthresh, dz_ub, ignore_cwise,
400  $ info )
401 *
402 * -- LAPACK computational routine (version 3.4.2) --
403 * -- LAPACK is a software package provided by Univ. of Tennessee, --
404 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
405 * September 2012
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