395 $ lda, af, ldaf, ipiv, colequ, c, b,
396 $ ldb, y, ldy, berr_out, n_norms,
397 $ errs_n, errs_c, res, ayb, dy,
398 $ y_tail, rcond, ithresh, rthresh,
399 $ dz_ub, ignore_cwise, info )
407 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
408 $ trans_type, n_norms
409 LOGICAL COLEQU, IGNORE_CWISE
411 DOUBLE PRECISION RTHRESH, DZ_UB
415 COMPLEX*16 A( lda, * ), AF( ldaf, * ), B( ldb, * ),
416 $ y( ldy, * ), res( * ), dy( * ), y_tail( * )
417 DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
418 $ errs_n( nrhs, * ), errs_c( nrhs, * )
425 INTEGER CNT, I, J, X_STATE, Z_STATE, Y_PREC_STATE
426 DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
427 $ dzrat, prevnormdx, prev_dz_z, dxratmax,
428 $ dzratmax, dx_x, dz_z, final_dx_x, final_dz_z,
429 $ eps, hugeval, incr_thresh
434 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
435 $ noprog_state, base_residual, extra_residual,
437 parameter( unstable_state = 0, working_state = 1,
440 parameter( base_residual = 0, extra_residual = 1,
442 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
443 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
444 INTEGER CMP_ERR_I, PIV_GROWTH_I
445 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
447 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
448 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
450 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
452 parameter( la_linrx_itref_i = 1,
453 $ la_linrx_ithresh_i = 2 )
454 parameter( la_linrx_cwise_i = 3 )
455 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
457 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
458 parameter( la_linrx_rcond_i = 3 )
464 DOUBLE PRECISION DLAMCH
465 CHARACTER CHLA_TRANSTYPE
468 INTRINSIC abs, max, min
471 DOUBLE PRECISION CABS1
474 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
478 IF ( info.NE.0 )
RETURN
479 trans = chla_transtype(trans_type)
480 eps = dlamch(
'Epsilon' )
481 hugeval = dlamch(
'Overflow' )
483 hugeval = hugeval * hugeval
485 incr_thresh = dble( n ) * eps
488 y_prec_state = extra_residual
489 IF ( y_prec_state .EQ. extra_y )
THEN
506 x_state = working_state
507 z_state = unstable_state
515 CALL zcopy( n, b( 1, j ), 1, res, 1 )
516 IF ( y_prec_state .EQ. base_residual )
THEN
517 CALL zgemv( trans, n, n, (-1.0d+0,0.0d+0), a, lda,
518 $ y( 1, j ), 1, (1.0d+0,0.0d+0), res, 1)
519 ELSE IF (y_prec_state .EQ. extra_residual)
THEN
520 CALL blas_zgemv_x( trans_type, n, n, (-1.0d+0,0.0d+0), a,
521 $ lda, y( 1, j ), 1, (1.0d+0,0.0d+0),
522 $ res, 1, prec_type )
524 CALL blas_zgemv2_x( trans_type, n, n, (-1.0d+0,0.0d+0),
525 $ a, lda, y(1, j), y_tail, 1, (1.0d+0,0.0d+0), res, 1,
530 CALL zcopy( n, res, 1, dy, 1 )
531 CALL zgetrs( trans, n, 1, af, ldaf, ipiv, dy, n, info )
542 yk = cabs1( y( i, j ) )
543 dyk = cabs1( dy( i ) )
545 IF ( yk .NE. 0.0d+0 )
THEN
546 dz_z = max( dz_z, dyk / yk )
547 ELSE IF ( dyk .NE. 0.0d+0 )
THEN
551 ymin = min( ymin, yk )
553 normy = max( normy, yk )
556 normx = max( normx, yk * c( i ) )
557 normdx = max( normdx, dyk * c( i ) )
560 normdx = max(normdx, dyk)
564 IF ( normx .NE. 0.0d+0 )
THEN
565 dx_x = normdx / normx
566 ELSE IF ( normdx .EQ. 0.0d+0 )
THEN
572 dxrat = normdx / prevnormdx
573 dzrat = dz_z / prev_dz_z
577 IF (.NOT.ignore_cwise
578 $ .AND. ymin*rcond .LT. incr_thresh*normy
579 $ .AND. y_prec_state .LT. extra_y )
582 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
583 $ x_state = working_state
584 IF ( x_state .EQ. working_state )
THEN
585 IF (dx_x .LE. eps)
THEN
587 ELSE IF ( dxrat .GT. rthresh )
THEN
588 IF ( y_prec_state .NE. extra_y )
THEN
591 x_state = noprog_state
594 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
596 IF ( x_state .GT. working_state ) final_dx_x = dx_x
599 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
600 $ z_state = working_state
601 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
602 $ z_state = working_state
603 IF ( z_state .EQ. working_state )
THEN
604 IF ( dz_z .LE. eps )
THEN
606 ELSE IF ( dz_z .GT. dz_ub )
THEN
607 z_state = unstable_state
610 ELSE IF ( dzrat .GT. rthresh )
THEN
611 IF ( y_prec_state .NE. extra_y )
THEN
614 z_state = noprog_state
617 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
619 IF ( z_state .GT. working_state ) final_dz_z = dz_z
626 IF ( x_state.NE.working_state )
THEN
627 IF ( ignore_cwise )
GOTO 666
628 IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
630 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 )
GOTO 666
633 IF ( incr_prec )
THEN
635 y_prec_state = y_prec_state + 1
646 IF ( y_prec_state .LT. extra_y )
THEN
647 CALL zaxpy( n, (1.0d+0,0.0d+0), dy, 1, y(1,j), 1 )
658 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
659 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
663 IF (n_norms .GE. 1)
THEN
664 errs_n( j, la_linrx_err_i ) = final_dx_x / (1 - dxratmax)
667 IF ( n_norms .GE. 2 )
THEN
668 errs_c( j, la_linrx_err_i ) = final_dz_z / (1 - dzratmax)
679 CALL zcopy( n, b( 1, j ), 1, res, 1 )
680 CALL zgemv( trans, n, n, (-1.0d+0,0.0d+0), a, lda, y(1,j), 1,
681 $ (1.0d+0,0.0d+0), res, 1 )
684 ayb( i ) = cabs1( b( i, j ) )
689 CALL zla_geamv ( trans_type, n, n, 1.0d+0,
690 $ a, lda, y(1, j), 1, 1.0d+0, ayb, 1 )
subroutine zla_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)
ZLA_GERFSX_EXTENDED
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zla_lin_berr(N, NZ, NRHS, RES, AYB, BERR)
ZLA_LIN_BERR computes a component-wise relative backward error.
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zla_wwaddw(N, X, Y, W)
ZLA_WWADDW adds a vector into a doubled-single vector.
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
double precision function dlamch(CMACH)
DLAMCH
character *1 function chla_transtype(TRANS)
CHLA_TRANSTYPE
subroutine zla_geamv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZLA_GEAMV computes a matrix-vector product using a general matrix to calculate error bounds...