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,
403 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
404 $ TRANS_TYPE, N_NORMS, ITHRESH
405 LOGICAL COLEQU, IGNORE_CWISE
410 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
411 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
412 REAL C( * ), AYB( * ), RCOND, BERR_OUT( * ),
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
429 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
430 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
432 PARAMETER ( UNSTABLE_STATE = 0, working_state = 1,
433 $ conv_state = 2, noprog_state = 3 )
434 parameter( base_residual = 0, extra_residual = 1,
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,
441 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
442 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
444 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_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,
451 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
452 parameter( la_linrx_rcond_i = 3 )
460 CHARACTER CHLA_TRANSTYPE
463 INTRINSIC abs, max, min
467 IF ( info.NE.0 )
RETURN
468 trans = chla_transtype(trans_type)
469 eps = slamch(
'Epsilon' )
470 hugeval = slamch(
'Overflow' )
472 hugeval = hugeval * hugeval
474 incr_thresh = real( n ) * eps
477 y_prec_state = extra_residual
478 IF ( y_prec_state .EQ. extra_y )
THEN
495 x_state = working_state
496 z_state = unstable_state
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,
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 )
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 )
517 CALL scopy( n, res, 1, dy, 1 )
518 CALL sgetrs( trans, n, 1, af, ldaf, ipiv, dy, n, info )
529 yk = abs( y( i, j ) )
532 IF ( yk .NE. 0.0 )
THEN
533 dz_z = max( dz_z, dyk / yk )
534 ELSE IF ( dyk .NE. 0.0 )
THEN
538 ymin = min( ymin, yk )
540 normy = max( normy, yk )
543 normx = max( normx, yk * c( i ) )
544 normdx = max( normdx, dyk * c( i ) )
547 normdx = max( normdx, dyk )
551 IF ( normx .NE. 0.0 )
THEN
552 dx_x = normdx / normx
553 ELSE IF ( normdx .EQ. 0.0 )
THEN
559 dxrat = normdx / prevnormdx
560 dzrat = dz_z / prev_dz_z
564 IF (.NOT.ignore_cwise
565 $ .AND. ymin*rcond .LT. incr_thresh*normy
566 $ .AND. y_prec_state .LT. extra_y)
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
574 ELSE IF ( dxrat .GT. rthresh )
THEN
575 IF ( y_prec_state .NE. extra_y )
THEN
578 x_state = noprog_state
581 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
583 IF ( x_state .GT. working_state ) final_dx_x = dx_x
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
593 ELSE IF ( dz_z .GT. dz_ub )
THEN
594 z_state = unstable_state
597 ELSE IF ( dzrat .GT. rthresh )
THEN
598 IF ( y_prec_state .NE. extra_y )
THEN
601 z_state = noprog_state
604 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
606 IF ( z_state .GT. working_state ) final_dz_z = dz_z
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 )
617 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 )
GOTO 666
620 IF ( incr_prec )
THEN
622 y_prec_state = y_prec_state + 1
633 IF ( y_prec_state .LT. extra_y )
THEN
634 CALL saxpy( n, 1.0, dy, 1, y( 1, j ), 1 )
645 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
646 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
650 IF (n_norms .GE. 1)
THEN
651 errs_n( j, la_linrx_err_i ) =
652 $ final_dx_x / (1 - dxratmax)
654 IF ( n_norms .GE. 2 )
THEN
655 errs_c( j, la_linrx_err_i ) =
656 $ final_dz_z / (1 - dzratmax)
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,
672 ayb( i ) = abs( b( i, j ) )
678 $ a, lda, y(1, j), 1, 1.0, ayb, 1 )