LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
sla_gerfsx_extended.f
Go to the documentation of this file.
1*> \brief \b SLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general 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*> Download SLA_GERFSX_EXTENDED + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_gerfsx_extended.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_gerfsx_extended.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_gerfsx_extended.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, NRHS, A,
20* LDA, AF, LDAF, IPIV, COLEQU, C, B,
21* LDB, Y, LDY, BERR_OUT, N_NORMS,
22* ERRS_N, ERRS_C, RES,
23* AYB, DY, Y_TAIL, RCOND, ITHRESH,
24* RTHRESH, DZ_UB, IGNORE_CWISE,
25* INFO )
26*
27* .. Scalar Arguments ..
28* INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
29* $ TRANS_TYPE, N_NORMS, ITHRESH
30* LOGICAL COLEQU, IGNORE_CWISE
31* REAL RTHRESH, DZ_UB
32* ..
33* .. Array Arguments ..
34* INTEGER IPIV( * )
35* REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
36* $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
37* REAL C( * ), AYB( * ), RCOND, BERR_OUT( * ),
38* $ ERRS_N( NRHS, * ),
39* $ ERRS_C( NRHS, * )
40* ..
41*
42*
43*> \par Purpose:
44* =============
45*>
46*> \verbatim
47*>
48*> SLA_GERFSX_EXTENDED improves the computed solution to a system of
49*> linear equations by performing extra-precise iterative refinement
50*> and provides error bounds and backward error estimates for the solution.
51*> This subroutine is called by SGERFSX to perform iterative refinement.
52*> In addition to normwise error bound, the code provides maximum
53*> componentwise error bound if possible. See comments for ERRS_N
54*> and ERRS_C for details of the error bounds. Note that this
55*> subroutine is only responsible for setting the second fields of
56*> ERRS_N and ERRS_C.
57*> \endverbatim
58*
59* Arguments:
60* ==========
61*
62*> \param[in] PREC_TYPE
63*> \verbatim
64*> PREC_TYPE is INTEGER
65*> Specifies the intermediate precision to be used in refinement.
66*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
67*> = 'S': Single
68*> = 'D': Double
69*> = 'I': Indigenous
70*> = 'X' or 'E': Extra
71*> \endverbatim
72*>
73*> \param[in] TRANS_TYPE
74*> \verbatim
75*> TRANS_TYPE is INTEGER
76*> Specifies the transposition operation on A.
77*> The value is defined by ILATRANS(T) where T is a CHARACTER and T
78*> = 'N': No transpose
79*> = 'T': Transpose
80*> = 'C': Conjugate transpose
81*> \endverbatim
82*>
83*> \param[in] N
84*> \verbatim
85*> N is INTEGER
86*> The number of linear equations, i.e., the order of the
87*> matrix A. N >= 0.
88*> \endverbatim
89*>
90*> \param[in] NRHS
91*> \verbatim
92*> NRHS is INTEGER
93*> The number of right-hand-sides, i.e., the number of columns of the
94*> matrix B.
95*> \endverbatim
96*>
97*> \param[in] A
98*> \verbatim
99*> A is REAL array, dimension (LDA,N)
100*> On entry, the N-by-N matrix A.
101*> \endverbatim
102*>
103*> \param[in] LDA
104*> \verbatim
105*> LDA is INTEGER
106*> The leading dimension of the array A. LDA >= max(1,N).
107*> \endverbatim
108*>
109*> \param[in] AF
110*> \verbatim
111*> AF is REAL array, dimension (LDAF,N)
112*> The factors L and U from the factorization
113*> A = P*L*U as computed by SGETRF.
114*> \endverbatim
115*>
116*> \param[in] LDAF
117*> \verbatim
118*> LDAF is INTEGER
119*> The leading dimension of the array AF. LDAF >= max(1,N).
120*> \endverbatim
121*>
122*> \param[in] IPIV
123*> \verbatim
124*> IPIV is INTEGER array, dimension (N)
125*> The pivot indices from the factorization A = P*L*U
126*> as computed by SGETRF; row i of the matrix was interchanged
127*> with row IPIV(i).
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 SGETRS.
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 ERRS_N
190*> and ERRS_C).
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] ERRS_N
196*> \verbatim
197*> ERRS_N 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 ERRS_N(i,:) corresponds to the ith
212*> right-hand side.
213*>
214*> The second index in ERRS_N(:,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] ERRS_C
242*> \verbatim
243*> ERRS_C 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*> ERRS_C 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 ERRS_C(i,:) corresponds to the ith
262*> right-hand side.
263*>
264*> The second index in ERRS_C(:,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*> ERRS_N and ERRS_C 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 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 SGETRS 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_gerfsx_extended
387*
388* =====================================================================
389 SUBROUTINE sla_gerfsx_extended( PREC_TYPE, TRANS_TYPE, N, NRHS,
390 $ A,
391 $ LDA, AF, LDAF, IPIV, COLEQU, C, B,
392 $ LDB, Y, LDY, BERR_OUT, N_NORMS,
393 $ ERRS_N, ERRS_C, RES,
394 $ AYB, DY, Y_TAIL, RCOND, ITHRESH,
395 $ RTHRESH, DZ_UB, IGNORE_CWISE,
396 $ INFO )
397*
398* -- LAPACK computational routine --
399* -- LAPACK is a software package provided by Univ. of Tennessee, --
400* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
401*
402* .. Scalar Arguments ..
403 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
404 $ TRANS_TYPE, N_NORMS, ITHRESH
405 LOGICAL COLEQU, IGNORE_CWISE
406 REAL RTHRESH, DZ_UB
407* ..
408* .. Array Arguments ..
409 INTEGER IPIV( * )
410 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
411 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
412 REAL C( * ), AYB( * ), RCOND, BERR_OUT( * ),
413 $ errs_n( nrhs, * ),
414 $ errs_c( nrhs, * )
415* ..
416*
417* =====================================================================
418*
419* .. Local Scalars ..
420 CHARACTER TRANS
421 INTEGER CNT, I, J, X_STATE, Z_STATE, Y_PREC_STATE
422 REAL YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
423 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
424 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
425 $ EPS, HUGEVAL, INCR_THRESH
426 LOGICAL INCR_PREC
427* ..
428* .. Parameters ..
429 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
430 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
431 $ EXTRA_Y
432 PARAMETER ( UNSTABLE_STATE = 0, working_state = 1,
433 $ conv_state = 2, noprog_state = 3 )
434 parameter( base_residual = 0, extra_residual = 1,
435 $ extra_y = 2 )
436 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
437 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
438 INTEGER CMP_ERR_I, PIV_GROWTH_I
439 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
440 $ berr_i = 3 )
441 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
442 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
443 $ piv_growth_i = 9 )
444 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
445 $ la_linrx_cwise_i
446 parameter( la_linrx_itref_i = 1,
447 $ la_linrx_ithresh_i = 2 )
448 parameter( la_linrx_cwise_i = 3 )
449 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
450 $ LA_LINRX_RCOND_I
451 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
452 parameter( la_linrx_rcond_i = 3 )
453* ..
454* .. External Subroutines ..
455 EXTERNAL saxpy, scopy, sgetrs, sgemv,
456 $ blas_sgemv_x,
457 $ blas_sgemv2_x, sla_geamv, sla_wwaddw, slamch,
459 REAL SLAMCH
460 CHARACTER CHLA_TRANSTYPE
461* ..
462* .. Intrinsic Functions ..
463 INTRINSIC abs, max, min
464* ..
465* .. Executable Statements ..
466*
467 IF ( info.NE.0 ) RETURN
468 trans = chla_transtype(trans_type)
469 eps = slamch( 'Epsilon' )
470 hugeval = slamch( 'Overflow' )
471* Force HUGEVAL to Inf
472 hugeval = hugeval * hugeval
473* Using HUGEVAL may lead to spurious underflows.
474 incr_thresh = real( n ) * eps
475*
476 DO j = 1, nrhs
477 y_prec_state = extra_residual
478 IF ( y_prec_state .EQ. extra_y ) THEN
479 DO i = 1, n
480 y_tail( i ) = 0.0
481 END DO
482 END IF
483
484 dxrat = 0.0
485 dxratmax = 0.0
486 dzrat = 0.0
487 dzratmax = 0.0
488 final_dx_x = hugeval
489 final_dz_z = hugeval
490 prevnormdx = hugeval
491 prev_dz_z = hugeval
492 dz_z = hugeval
493 dx_x = hugeval
494
495 x_state = working_state
496 z_state = unstable_state
497 incr_prec = .false.
498
499 DO cnt = 1, ithresh
500*
501* Compute residual RES = B_s - op(A_s) * Y,
502* op(A) = A, A**T, or A**H depending on TRANS (and type).
503*
504 CALL scopy( n, b( 1, j ), 1, res, 1 )
505 IF ( y_prec_state .EQ. base_residual ) THEN
506 CALL sgemv( trans, n, n, -1.0, a, lda, y( 1, j ), 1,
507 $ 1.0, res, 1 )
508 ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
509 CALL blas_sgemv_x( trans_type, n, n, -1.0, a, lda,
510 $ y( 1, j ), 1, 1.0, res, 1, prec_type )
511 ELSE
512 CALL blas_sgemv2_x( trans_type, n, n, -1.0, a, lda,
513 $ y( 1, j ), y_tail, 1, 1.0, res, 1, prec_type )
514 END IF
515
516! XXX: RES is no longer needed.
517 CALL scopy( n, res, 1, dy, 1 )
518 CALL sgetrs( trans, n, 1, af, ldaf, ipiv, dy, n, info )
519*
520* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
521*
522 normx = 0.0
523 normy = 0.0
524 normdx = 0.0
525 dz_z = 0.0
526 ymin = hugeval
527*
528 DO i = 1, n
529 yk = abs( y( i, j ) )
530 dyk = abs( dy( i ) )
531
532 IF ( yk .NE. 0.0 ) THEN
533 dz_z = max( dz_z, dyk / yk )
534 ELSE IF ( dyk .NE. 0.0 ) THEN
535 dz_z = hugeval
536 END IF
537
538 ymin = min( ymin, yk )
539
540 normy = max( normy, yk )
541
542 IF ( colequ ) THEN
543 normx = max( normx, yk * c( i ) )
544 normdx = max( normdx, dyk * c( i ) )
545 ELSE
546 normx = normy
547 normdx = max( normdx, dyk )
548 END IF
549 END DO
550
551 IF ( normx .NE. 0.0 ) THEN
552 dx_x = normdx / normx
553 ELSE IF ( normdx .EQ. 0.0 ) THEN
554 dx_x = 0.0
555 ELSE
556 dx_x = hugeval
557 END IF
558
559 dxrat = normdx / prevnormdx
560 dzrat = dz_z / prev_dz_z
561*
562* Check termination criteria
563*
564 IF (.NOT.ignore_cwise
565 $ .AND. ymin*rcond .LT. incr_thresh*normy
566 $ .AND. y_prec_state .LT. extra_y)
567 $ incr_prec = .true.
568
569 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
570 $ x_state = working_state
571 IF ( x_state .EQ. working_state ) THEN
572 IF ( dx_x .LE. eps ) THEN
573 x_state = conv_state
574 ELSE IF ( dxrat .GT. rthresh ) THEN
575 IF ( y_prec_state .NE. extra_y ) THEN
576 incr_prec = .true.
577 ELSE
578 x_state = noprog_state
579 END IF
580 ELSE
581 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
582 END IF
583 IF ( x_state .GT. working_state ) final_dx_x = dx_x
584 END IF
585
586 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
587 $ z_state = working_state
588 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
589 $ z_state = working_state
590 IF ( z_state .EQ. working_state ) THEN
591 IF ( dz_z .LE. eps ) THEN
592 z_state = conv_state
593 ELSE IF ( dz_z .GT. dz_ub ) THEN
594 z_state = unstable_state
595 dzratmax = 0.0
596 final_dz_z = hugeval
597 ELSE IF ( dzrat .GT. rthresh ) THEN
598 IF ( y_prec_state .NE. extra_y ) THEN
599 incr_prec = .true.
600 ELSE
601 z_state = noprog_state
602 END IF
603 ELSE
604 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
605 END IF
606 IF ( z_state .GT. working_state ) final_dz_z = dz_z
607 END IF
608*
609* Exit if both normwise and componentwise stopped working,
610* but if componentwise is unstable, let it go at least two
611* iterations.
612*
613 IF ( x_state.NE.working_state ) THEN
614 IF ( ignore_cwise) GOTO 666
615 IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
616 $ GOTO 666
617 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 ) GOTO 666
618 END IF
619
620 IF ( incr_prec ) THEN
621 incr_prec = .false.
622 y_prec_state = y_prec_state + 1
623 DO i = 1, n
624 y_tail( i ) = 0.0
625 END DO
626 END IF
627
628 prevnormdx = normdx
629 prev_dz_z = dz_z
630*
631* Update solution.
632*
633 IF ( y_prec_state .LT. extra_y ) THEN
634 CALL saxpy( n, 1.0, dy, 1, y( 1, j ), 1 )
635 ELSE
636 CALL sla_wwaddw( n, y( 1, j ), y_tail, dy )
637 END IF
638
639 END DO
640* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
641 666 CONTINUE
642*
643* Set final_* when cnt hits ithresh.
644*
645 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
646 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
647*
648* Compute error bounds
649*
650 IF (n_norms .GE. 1) THEN
651 errs_n( j, la_linrx_err_i ) =
652 $ final_dx_x / (1 - dxratmax)
653 END IF
654 IF ( n_norms .GE. 2 ) THEN
655 errs_c( j, la_linrx_err_i ) =
656 $ final_dz_z / (1 - dzratmax)
657 END IF
658*
659* Compute componentwise relative backward error from formula
660* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
661* where abs(Z) is the componentwise absolute value of the matrix
662* or vector Z.
663*
664* Compute residual RES = B_s - op(A_s) * Y,
665* op(A) = A, A**T, or A**H depending on TRANS (and type).
666*
667 CALL scopy( n, b( 1, j ), 1, res, 1 )
668 CALL sgemv( trans, n, n, -1.0, a, lda, y(1,j), 1, 1.0, res,
669 $ 1 )
670
671 DO i = 1, n
672 ayb( i ) = abs( b( i, j ) )
673 END DO
674*
675* Compute abs(op(A_s))*abs(Y) + abs(B_s).
676*
677 CALL sla_geamv ( trans_type, n, n, 1.0,
678 $ a, lda, y(1, j), 1, 1.0, ayb, 1 )
679
680 CALL sla_lin_berr ( n, n, 1, res, ayb, berr_out( j ) )
681*
682* End of loop for each RHS.
683*
684 END DO
685*
686 RETURN
687*
688* End of SLA_GERFSX_EXTENDED
689*
690 END
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
subroutine sgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
SGETRS
Definition sgetrs.f:119
subroutine sla_geamv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SLA_GEAMV computes a matrix-vector product using a general matrix to calculate error bounds.
Definition sla_geamv.f:175
subroutine sla_gerfsx_extended(prec_type, trans_type, n, nrhs, a, lda, af, ldaf, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, errs_n, errs_c, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)
SLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general matric...
subroutine sla_lin_berr(n, nz, nrhs, res, ayb, berr)
SLA_LIN_BERR computes a component-wise relative backward error.
character *1 function chla_transtype(trans)
CHLA_TRANSTYPE
subroutine sla_wwaddw(n, x, y, w)
SLA_WWADDW adds a vector into a doubled-single vector.
Definition sla_wwaddw.f:79
real function slamch(cmach)
SLAMCH
Definition slamch.f:68