408 $ nrhs, ab, ldab, afb, ldafb, ipiv,
409 $ colequ, c, b, ldb, y, ldy,
410 $ berr_out, n_norms, err_bnds_norm,
411 $ err_bnds_comp, res, ayb, dy,
412 $ y_tail, rcond, ithresh, rthresh,
413 $ dz_ub, ignore_cwise, info )
421 INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS,
422 $ prec_type, trans_type, n_norms, ithresh
423 LOGICAL COLEQU, IGNORE_CWISE
428 REAL AB( ldab, * ), AFB( ldafb, * ), B( ldb, * ),
429 $ y( ldy, * ), res(*), dy(*), y_tail(*)
430 REAL C( * ), AYB(*), RCOND, BERR_OUT(*),
431 $ err_bnds_norm( nrhs, * ),
432 $ err_bnds_comp( nrhs, * )
439 INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE
440 REAL YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
441 $ dzrat, prevnormdx, prev_dz_z, dxratmax,
442 $ dzratmax, dx_x, dz_z, final_dx_x, final_dz_z,
443 $ eps, hugeval, incr_thresh
447 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
448 $ noprog_state, base_residual, extra_residual,
450 parameter ( unstable_state = 0, working_state = 1,
451 $ conv_state = 2, noprog_state = 3 )
452 parameter ( base_residual = 0, extra_residual = 1,
454 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
455 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
456 INTEGER CMP_ERR_I, PIV_GROWTH_I
457 parameter ( final_nrm_err_i = 1, final_cmp_err_i = 2,
459 parameter ( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
460 parameter ( cmp_rcond_i = 7, cmp_err_i = 8,
462 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
464 parameter ( la_linrx_itref_i = 1,
465 $ la_linrx_ithresh_i = 2 )
466 parameter ( la_linrx_cwise_i = 3 )
467 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
469 parameter ( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
470 parameter ( la_linrx_rcond_i = 3 )
477 CHARACTER CHLA_TRANSTYPE
480 INTRINSIC abs, max, min
484 IF (info.NE.0)
RETURN
485 trans = chla_transtype(trans_type)
486 eps = slamch(
'Epsilon' )
487 hugeval = slamch(
'Overflow' )
489 hugeval = hugeval * hugeval
491 incr_thresh =
REAL( N ) * EPS
495 y_prec_state = extra_residual
496 IF ( y_prec_state .EQ. extra_y )
THEN
513 x_state = working_state
514 z_state = unstable_state
522 CALL scopy( n, b( 1, j ), 1, res, 1 )
523 IF ( y_prec_state .EQ. base_residual )
THEN
524 CALL sgbmv( trans, m, n, kl, ku, -1.0, ab, ldab,
525 $ y( 1, j ), 1, 1.0, res, 1 )
526 ELSE IF ( y_prec_state .EQ. extra_residual )
THEN
527 CALL blas_sgbmv_x( trans_type, n, n, kl, ku,
528 $ -1.0, ab, ldab, y( 1, j ), 1, 1.0, res, 1,
531 CALL blas_sgbmv2_x( trans_type, n, n, kl, ku, -1.0,
532 $ ab, ldab, y( 1, j ), y_tail, 1, 1.0, res, 1,
537 CALL scopy( n, res, 1, dy, 1 )
538 CALL sgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, dy, n,
550 yk = abs( y( i, j ) )
553 IF ( yk .NE. 0.0 )
THEN
554 dz_z = max( dz_z, dyk / yk )
555 ELSE IF ( dyk .NE. 0.0 )
THEN
559 ymin = min( ymin, yk )
561 normy = max( normy, yk )
564 normx = max( normx, yk * c( i ) )
565 normdx = max( normdx, dyk * c( i ) )
568 normdx = max( normdx, dyk )
572 IF ( normx .NE. 0.0 )
THEN
573 dx_x = normdx / normx
574 ELSE IF ( normdx .EQ. 0.0 )
THEN
580 dxrat = normdx / prevnormdx
581 dzrat = dz_z / prev_dz_z
585 IF ( .NOT.ignore_cwise
586 $ .AND. ymin*rcond .LT. incr_thresh*normy
587 $ .AND. y_prec_state .LT. extra_y )
590 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
591 $ x_state = working_state
592 IF ( x_state .EQ. working_state )
THEN
593 IF ( dx_x .LE. eps )
THEN
595 ELSE IF ( dxrat .GT. rthresh )
THEN
596 IF ( y_prec_state .NE. extra_y )
THEN
599 x_state = noprog_state
602 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
604 IF ( x_state .GT. working_state ) final_dx_x = dx_x
607 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
608 $ z_state = working_state
609 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
610 $ z_state = working_state
611 IF ( z_state .EQ. working_state )
THEN
612 IF ( dz_z .LE. eps )
THEN
614 ELSE IF ( dz_z .GT. dz_ub )
THEN
615 z_state = unstable_state
618 ELSE IF ( dzrat .GT. rthresh )
THEN
619 IF ( y_prec_state .NE. extra_y )
THEN
622 z_state = noprog_state
625 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
627 IF ( z_state .GT. working_state ) final_dz_z = dz_z
634 IF ( x_state.NE.working_state )
THEN
635 IF ( ignore_cwise )
GOTO 666
636 IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
638 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 )
GOTO 666
641 IF ( incr_prec )
THEN
643 y_prec_state = y_prec_state + 1
654 IF (y_prec_state .LT. extra_y)
THEN
655 CALL saxpy( n, 1.0, dy, 1, y(1,j), 1 )
666 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
667 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
671 IF ( n_norms .GE. 1 )
THEN
672 err_bnds_norm( j, la_linrx_err_i ) =
673 $ final_dx_x / (1 - dxratmax)
675 IF (n_norms .GE. 2)
THEN
676 err_bnds_comp( j, la_linrx_err_i ) =
677 $ final_dz_z / (1 - dzratmax)
688 CALL scopy( n, b( 1, j ), 1, res, 1 )
689 CALL sgbmv(trans, n, n, kl, ku, -1.0, ab, ldab, y(1,j),
693 ayb( i ) = abs( b( i, j ) )
698 CALL sla_gbamv( trans_type, n, n, kl, ku, 1.0,
699 $ ab, ldab, y(1, j), 1, 1.0, ayb, 1 )
subroutine sla_gbamv(TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X, INCX, BETA, Y, INCY)
SLA_GBAMV performs a matrix-vector operation to calculate error bounds.
subroutine sla_wwaddw(N, X, Y, W)
SLA_WWADDW adds a vector into a doubled-single vector.
subroutine sgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGBMV
subroutine sla_lin_berr(N, NZ, NRHS, RES, AYB, BERR)
SLA_LIN_BERR computes a component-wise relative backward error.
subroutine sla_gbrfsx_extended(PREC_TYPE, TRANS_TYPE, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, 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_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded...
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
real function slamch(CMACH)
SLAMCH
subroutine sgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
SGBTRS
character *1 function chla_transtype(TRANS)
CHLA_TRANSTYPE
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY