Skip to content

Scale down in xLARFG and xLARFGP when |BETA| exceeds half the overflow threshold #1381

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rmlarsen wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Scale down in xLARFG and xLARFGP when |BETA| exceeds half the overflow threshold #1381

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 20 additions & 4 deletions SRC/clarfg.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,8 +118,9 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
REAL ONE, ZERO, HALF
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0,
$ HALF = 0.5E+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -131,7 +132,7 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN
INTRINSIC ABS, AIMAG, CMPLX, HUGE, REAL, SIGN
* ..
* .. External Subroutines ..
EXTERNAL CSCAL, CSSCAL
Expand Down Expand Up @@ -179,16 +180,31 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
XNORM = SCNRM2( N-1, X, INCX )
ALPHA = CMPLX( ALPHR, ALPHI )
BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
KNT = -1
CALL CSSCAL( N-1, SAFMIN, X, INCX )
ALPHI = ALPHI*SAFMIN
ALPHR = ALPHR*SAFMIN
XNORM = SCNRM2( N-1, X, INCX )
ALPHA = CMPLX( ALPHR, ALPHI )
BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
END IF
TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA )
CALL CSCAL( N-1, ALPHA, X, INCX )
*
* If ALPHA is subnormal, it may lose relative accuracy
* Undo the scaling. If ALPHA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SAFMIN
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*RSAFMN
ALPHA = BETA
END IF
*
Expand Down
24 changes: 20 additions & 4 deletions SRC/clarfgp.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,8 +116,9 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
REAL TWO, ONE, ZERO
PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 )
REAL TWO, ONE, ZERO, HALF
PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0,
$ HALF = 0.5E+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -131,7 +132,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
$ CLADIV
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN
INTRINSIC ABS, AIMAG, CMPLX, HUGE, REAL, SIGN
* ..
* .. External Subroutines ..
EXTERNAL CSCAL, CSSCAL
Expand Down Expand Up @@ -193,6 +194,18 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
XNORM = SCNRM2( N-1, X, INCX )
ALPHA = CMPLX( ALPHR, ALPHI )
BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
KNT = -1
CALL CSSCAL( N-1, SMLNUM, X, INCX )
ALPHI = ALPHI*SMLNUM
ALPHR = ALPHR*SMLNUM
XNORM = SCNRM2( N-1, X, INCX )
ALPHA = CMPLX( ALPHR, ALPHI )
BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
END IF
SAVEALPHA = ALPHA
ALPHA = ALPHA + BETA
Expand Down Expand Up @@ -245,11 +258,14 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
*
END IF
*
* If BETA is subnormal, it may lose relative accuracy
* Undo the scaling. If BETA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SMLNUM
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*BIGNUM
ALPHA = BETA
END IF
*
Expand Down
23 changes: 19 additions & 4 deletions SRC/dlarfg.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,8 +118,9 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
DOUBLE PRECISION ONE, ZERO, HALF
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0,
$ HALF = 0.5D+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -130,7 +131,7 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
EXTERNAL DLAMCH, DLAPY2, DNRM2
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SIGN
INTRINSIC ABS, HUGE, SIGN
* ..
* .. External Subroutines ..
EXTERNAL DSCAL
Expand Down Expand Up @@ -173,15 +174,29 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
*
XNORM = DNRM2( N-1, X, INCX )
BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
RSAFMN = ONE / SAFMIN
KNT = -1
CALL DSCAL( N-1, SAFMIN, X, INCX )
ALPHA = ALPHA*SAFMIN
XNORM = DNRM2( N-1, X, INCX )
BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
END IF
TAU = ( BETA-ALPHA ) / BETA
CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
*
* If ALPHA is subnormal, it may lose relative accuracy
* Undo the scaling. If ALPHA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SAFMIN
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*RSAFMN
ALPHA = BETA
END IF
*
Expand Down
23 changes: 19 additions & 4 deletions SRC/dlarfgp.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,8 +116,9 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION TWO, ONE, ZERO
PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
DOUBLE PRECISION TWO, ONE, ZERO, HALF
PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0,
$ HALF = 0.5D+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -128,7 +129,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
EXTERNAL DLAMCH, DLAPY2, DNRM2
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SIGN
INTRINSIC ABS, HUGE, SIGN
* ..
* .. External Subroutines ..
EXTERNAL DSCAL
Expand Down Expand Up @@ -185,6 +186,17 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
*
XNORM = DNRM2( N-1, X, INCX )
BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
BIGNUM = ONE / SMLNUM
KNT = -1
CALL DSCAL( N-1, SMLNUM, X, INCX )
ALPHA = ALPHA*SMLNUM
XNORM = DNRM2( N-1, X, INCX )
BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
END IF
SAVEALPHA = ALPHA
ALPHA = ALPHA + BETA
Expand Down Expand Up @@ -229,11 +241,14 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
*
END IF
*
* If BETA is subnormal, it may lose relative accuracy
* Undo the scaling. If BETA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SMLNUM
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*BIGNUM
ALPHA = BETA
END IF
*
Expand Down
23 changes: 19 additions & 4 deletions SRC/slarfg.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,8 +118,9 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
REAL ONE, ZERO, HALF
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0,
$ HALF = 0.5E+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -130,7 +131,7 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
EXTERNAL SLAMCH, SLAPY2, SNRM2
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SIGN
INTRINSIC ABS, HUGE, SIGN
* ..
* .. External Subroutines ..
EXTERNAL SSCAL
Expand Down Expand Up @@ -173,15 +174,29 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
*
XNORM = SNRM2( N-1, X, INCX )
BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
RSAFMN = ONE / SAFMIN
KNT = -1
CALL SSCAL( N-1, SAFMIN, X, INCX )
ALPHA = ALPHA*SAFMIN
XNORM = SNRM2( N-1, X, INCX )
BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
END IF
TAU = ( BETA-ALPHA ) / BETA
CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
*
* If ALPHA is subnormal, it may lose relative accuracy
* Undo the scaling. If ALPHA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SAFMIN
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*RSAFMN
ALPHA = BETA
END IF
*
Expand Down
23 changes: 19 additions & 4 deletions SRC/slarfgp.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,8 +116,9 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
REAL TWO, ONE, ZERO
PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 )
REAL TWO, ONE, ZERO, HALF
PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0,
$ HALF = 0.5E+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -128,7 +129,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
EXTERNAL SLAMCH, SLAPY2, SNRM2
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SIGN
INTRINSIC ABS, HUGE, SIGN
* ..
* .. External Subroutines ..
EXTERNAL SSCAL
Expand Down Expand Up @@ -185,6 +186,17 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
*
XNORM = SNRM2( N-1, X, INCX )
BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
BIGNUM = ONE / SMLNUM
KNT = -1
CALL SSCAL( N-1, SMLNUM, X, INCX )
ALPHA = ALPHA*SMLNUM
XNORM = SNRM2( N-1, X, INCX )
BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
END IF
SAVEALPHA = ALPHA
ALPHA = ALPHA + BETA
Expand Down Expand Up @@ -229,11 +241,14 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
*
END IF
*
* If BETA is subnormal, it may lose relative accuracy
* Undo the scaling. If BETA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SMLNUM
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*BIGNUM
ALPHA = BETA
END IF
*
Expand Down
24 changes: 20 additions & 4 deletions SRC/zlarfg.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,8 +118,9 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU )
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
DOUBLE PRECISION ONE, ZERO, HALF
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0,
$ HALF = 0.5D+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
Expand All @@ -131,7 +132,7 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU )
EXTERNAL DLAMCH, DLAPY3, DZNRM2, ZLADIV
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN
INTRINSIC ABS, DBLE, DCMPLX, DIMAG, HUGE, SIGN
* ..
* .. External Subroutines ..
EXTERNAL ZDSCAL, ZSCAL
Expand Down Expand Up @@ -179,16 +180,31 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU )
XNORM = DZNRM2( N-1, X, INCX )
ALPHA = DCMPLX( ALPHR, ALPHI )
BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN
*
* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when
* |BETA| > HUGE/2; scale X down and recompute them.
*
KNT = -1
CALL ZDSCAL( N-1, SAFMIN, X, INCX )
ALPHI = ALPHI*SAFMIN
ALPHR = ALPHR*SAFMIN
XNORM = DZNRM2( N-1, X, INCX )
ALPHA = DCMPLX( ALPHR, ALPHI )
BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
END IF
TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA )
CALL ZSCAL( N-1, ALPHA, X, INCX )
*
* If ALPHA is subnormal, it may lose relative accuracy
* Undo the scaling. If ALPHA is subnormal, it may lose relative
* accuracy
*
DO 20 J = 1, KNT
BETA = BETA*SAFMIN
20 CONTINUE
IF( KNT.LT.0 )
$ BETA = BETA*RSAFMN
ALPHA = BETA
END IF
*
Expand Down
Loading
Loading