xLARFGP: avoid overflows

christoph-conrads · christoph-conrads · commit 64897c447470 · 2023-11-10T20:08:50.000+01:00
Avoid overflows when XNORM << ABS(ALPHA) << 1. fixes #934
diff --git a/SRC/clarfgp.f b/SRC/clarfgp.f
@@ -122,7 +122,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
 *     ..
 *     .. Local Scalars ..
       INTEGER            J, KNT
-      REAL               ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
+      REAL               ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM
       COMPLEX            SAVEALPHA
 *     ..
 *     .. External Functions ..
@@ -143,11 +143,12 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
          RETURN
       END IF
 *
+      EPS = SLAMCH( 'Precision' )
       XNORM = SCNRM2( N-1, X, INCX )
       ALPHR = REAL( ALPHA )
       ALPHI = AIMAG( ALPHA )
 *
-      IF( XNORM.EQ.ZERO ) THEN
+      IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
 *
 *        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
 *
diff --git a/SRC/dlarfgp.f b/SRC/dlarfgp.f
@@ -122,7 +122,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
 *     ..
 *     .. Local Scalars ..
       INTEGER            J, KNT
-      DOUBLE PRECISION   BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
+      DOUBLE PRECISION   BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
@@ -141,11 +141,12 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
          RETURN
       END IF
 *
+      EPS = DLAMCH( 'Precision' )
       XNORM = DNRM2( N-1, X, INCX )
 *
-      IF( XNORM.EQ.ZERO ) THEN
+      IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
 *
-*        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0
+*        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0.
 *
          IF( ALPHA.GE.ZERO ) THEN
 *           When TAU.eq.ZERO, the vector is special-cased to be
diff --git a/SRC/slarfgp.f b/SRC/slarfgp.f
@@ -122,7 +122,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
 *     ..
 *     .. Local Scalars ..
       INTEGER            J, KNT
-      REAL               BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
+      REAL               BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM
 *     ..
 *     .. External Functions ..
       REAL               SLAMCH, SLAPY2, SNRM2
@@ -141,9 +141,10 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
          RETURN
       END IF
 *
+      EPS = SLAMCH( 'Precision' )
       XNORM = SNRM2( N-1, X, INCX )
 *
-      IF( XNORM.EQ.ZERO ) THEN
+      IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
 *
 *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0.
 *
diff --git a/SRC/zlarfgp.f b/SRC/zlarfgp.f
@@ -122,7 +122,7 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
 *     ..
 *     .. Local Scalars ..
       INTEGER            J, KNT
-      DOUBLE PRECISION   ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
+      DOUBLE PRECISION   ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM
       COMPLEX*16         SAVEALPHA
 *     ..
 *     .. External Functions ..
@@ -143,11 +143,12 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
          RETURN
       END IF
 *
+      EPS = DLAMCH( 'Precision' )
       XNORM = DZNRM2( N-1, X, INCX )
       ALPHR = DBLE( ALPHA )
       ALPHI = DIMAG( ALPHA )
 *
-      IF( XNORM.EQ.ZERO ) THEN
+      IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
 *
 *        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
 *

Original file line number	Diff line number	Diff line change
`@@ -122,7 +122,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )`
`122`	`122`	`* ..`
`123`	`123`	`* .. Local Scalars ..`
`124`	`124`	`INTEGER J, KNT`
`125`		`- REAL ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM`
	`125`	`+ REAL ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM`
`126`	`126`	`COMPLEX SAVEALPHA`
`127`	`127`	`* ..`
`128`	`128`	`* .. External Functions ..`
`@@ -143,11 +143,12 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )`
`143`	`143`	`RETURN`
`144`	`144`	`END IF`
`145`	`145`	`*`
	`146`	`+ EPS = SLAMCH( 'Precision' )`
`146`	`147`	`XNORM = SCNRM2( N-1, X, INCX )`
`147`	`148`	`ALPHR = REAL( ALPHA )`
`148`	`149`	`ALPHI = AIMAG( ALPHA )`
`149`	`150`	`*`
`150`		`- IF( XNORM.EQ.ZERO ) THEN`
	`151`	`+ IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN`
`151`	`152`	`*`
`152`	`153`	`* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.`
`153`	`154`	`*`
Original file line number	Diff line number	Diff line change
`@@ -122,7 +122,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )`
`122`	`122`	`* ..`
`123`	`123`	`* .. Local Scalars ..`
`124`	`124`	`INTEGER J, KNT`
`125`		`- DOUBLE PRECISION BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM`
	`125`	`+ DOUBLE PRECISION BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM`
`126`	`126`	`* ..`
`127`	`127`	`* .. External Functions ..`
`128`	`128`	`DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2`
`@@ -141,11 +141,12 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )`
`141`	`141`	`RETURN`
`142`	`142`	`END IF`
`143`	`143`	`*`
	`144`	`+ EPS = DLAMCH( 'Precision' )`
`144`	`145`	`XNORM = DNRM2( N-1, X, INCX )`
`145`	`146`	`*`
`146`		`- IF( XNORM.EQ.ZERO ) THEN`
	`147`	`+ IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN`
`147`	`148`	`*`
`148`		`-* H = [+/-1, 0; I], sign chosen so ALPHA >= 0`
	`149`	`+* H = [+/-1, 0; I], sign chosen so ALPHA >= 0.`
`149`	`150`	`*`
`150`	`151`	`IF( ALPHA.GE.ZERO ) THEN`
`151`	`152`	`* When TAU.eq.ZERO, the vector is special-cased to be`
Original file line number	Diff line number	Diff line change
`@@ -122,7 +122,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )`
`122`	`122`	`* ..`
`123`	`123`	`* .. Local Scalars ..`
`124`	`124`	`INTEGER J, KNT`
`125`		`- REAL BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM`
	`125`	`+ REAL BETA, BIGNUM, EPS, SAVEALPHA, SMLNUM, XNORM`
`126`	`126`	`* ..`
`127`	`127`	`* .. External Functions ..`
`128`	`128`	`REAL SLAMCH, SLAPY2, SNRM2`
`@@ -141,9 +141,10 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )`
`141`	`141`	`RETURN`
`142`	`142`	`END IF`
`143`	`143`	`*`
	`144`	`+ EPS = SLAMCH( 'Precision' )`
`144`	`145`	`XNORM = SNRM2( N-1, X, INCX )`
`145`	`146`	`*`
`146`		`- IF( XNORM.EQ.ZERO ) THEN`
	`147`	`+ IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN`
`147`	`148`	`*`
`148`	`149`	`* H = [+/-1, 0; I], sign chosen so ALPHA >= 0.`
`149`	`150`	`*`
Original file line number	Diff line number	Diff line change
`@@ -122,7 +122,7 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )`
`122`	`122`	`* ..`
`123`	`123`	`* .. Local Scalars ..`
`124`	`124`	`INTEGER J, KNT`
`125`		`- DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM`
	`125`	`+ DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM`
`126`	`126`	`COMPLEX*16 SAVEALPHA`
`127`	`127`	`* ..`
`128`	`128`	`* .. External Functions ..`
`@@ -143,11 +143,12 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )`
`143`	`143`	`RETURN`
`144`	`144`	`END IF`
`145`	`145`	`*`
	`146`	`+ EPS = DLAMCH( 'Precision' )`
`146`	`147`	`XNORM = DZNRM2( N-1, X, INCX )`
`147`	`148`	`ALPHR = DBLE( ALPHA )`
`148`	`149`	`ALPHI = DIMAG( ALPHA )`
`149`	`150`	`*`
`150`		`- IF( XNORM.EQ.ZERO ) THEN`
	`151`	`+ IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN`
`151`	`152`	`*`
`152`	`153`	`* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.`
`153`	`154`	`*`