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,
405 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
408 LOGICAL COLEQU, IGNORE_CWISE
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, * )
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
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,
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,
443 parameter ( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
444 parameter ( cmp_rcond_i = 7, cmp_err_i = 8,
446 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_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,
453 parameter ( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
454 parameter ( la_linrx_rcond_i = 3 )
468 INTRINSIC abs, max, min
473 upper = lsame( uplo,
'U' )
474 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
476 ELSE IF( n.LT.0 )
THEN
478 ELSE IF( nrhs.LT.0 )
THEN
480 ELSE IF( lda.LT.max( 1, n ) )
THEN
482 ELSE IF( ldaf.LT.max( 1, n ) )
THEN
484 ELSE IF( ldb.LT.max( 1, n ) )
THEN
486 ELSE IF( ldy.LT.max( 1, n ) )
THEN
490 CALL xerbla(
'SLA_SYRFSX_EXTENDED', -info )
493 eps = slamch(
'Epsilon' )
494 hugeval = slamch(
'Overflow' )
496 hugeval = hugeval * hugeval
498 incr_thresh =
REAL( n )*EPS
500 IF ( lsame( uplo,
'L' ) )
THEN
501 uplo2 = ilauplo(
'L' )
503 uplo2 = ilauplo(
'U' )
507 y_prec_state = extra_residual
508 IF ( y_prec_state .EQ. extra_y )
THEN
525 x_state = working_state
526 z_state = unstable_state
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,
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 )
542 CALL blas_ssymv2_x(uplo2, n, -1.0, a, lda,
543 $ y(1, j), y_tail, 1, 1.0, res, 1, prec_type)
547 CALL scopy( n, res, 1, dy, 1 )
548 CALL ssytrs( uplo, n, 1, af, ldaf, ipiv, dy, n, info )
559 yk = abs( y( i, j ) )
562 IF ( yk .NE. 0.0 )
THEN
563 dz_z = max( dz_z, dyk / yk )
564 ELSE IF ( dyk .NE. 0.0 )
THEN
568 ymin = min( ymin, yk )
570 normy = max( normy, yk )
573 normx = max( normx, yk * c( i ) )
574 normdx = max( normdx, dyk * c( i ) )
577 normdx = max(normdx, dyk)
581 IF ( normx .NE. 0.0 )
THEN
582 dx_x = normdx / normx
583 ELSE IF ( normdx .EQ. 0.0 )
THEN
589 dxrat = normdx / prevnormdx
590 dzrat = dz_z / prev_dz_z
594 IF ( ymin*rcond .LT. incr_thresh*normy
595 $ .AND. y_prec_state .LT. extra_y )
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
603 ELSE IF ( dxrat .GT. rthresh )
THEN
604 IF ( y_prec_state .NE. extra_y )
THEN
607 x_state = noprog_state
610 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
612 IF ( x_state .GT. working_state ) final_dx_x = dx_x
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
622 ELSE IF ( dz_z .GT. dz_ub )
THEN
623 z_state = unstable_state
626 ELSE IF ( dzrat .GT. rthresh )
THEN
627 IF ( y_prec_state .NE. extra_y )
THEN
630 z_state = noprog_state
633 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
635 IF ( z_state .GT. working_state ) final_dz_z = dz_z
638 IF ( x_state.NE.working_state.AND.
639 $ ( ignore_cwise.OR.z_state.NE.working_state ) )
642 IF ( incr_prec )
THEN
644 y_prec_state = y_prec_state + 1
655 IF (y_prec_state .LT. extra_y)
THEN
656 CALL saxpy( n, 1.0, dy, 1, y(1,j), 1 )
667 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
668 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
672 IF ( n_norms .GE. 1 )
THEN
673 err_bnds_norm( j, la_linrx_err_i ) =
674 $ final_dx_x / (1 - dxratmax)
676 IF ( n_norms .GE. 2 )
THEN
677 err_bnds_comp( j, la_linrx_err_i ) =
678 $ final_dz_z / (1 - dzratmax)
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 )
692 ayb( i ) = abs( b( i, j ) )
698 $ a, lda, y(1, j), 1, 1.0, ayb, 1 )
subroutine ssytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SSYTRS
subroutine sla_syamv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bou...
subroutine sla_wwaddw(N, X, Y, W)
SLA_WWADDW adds a vector into a doubled-single vector.
subroutine sla_lin_berr(N, NZ, NRHS, RES, AYB, BERR)
SLA_LIN_BERR computes a component-wise relative backward error.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
integer function ilauplo(UPLO)
ILAUPLO
subroutine ssymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SSYMV
subroutine sla_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)
SLA_SYRFSX_EXTENDED improves the computed solution to a system of linear equations for symmetric inde...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY