xLARFGP: avoid overflows

Avoid overflows when XNORM << ABS(ALPHA) << 1. fixes #934
Reference-LAPACK · Nov 10, 2023 · 64897c4 · 64897c4
1 parent 5b0bc5e
commit 64897c4
Show file tree

Hide file tree

Showing 4 changed files with 13 additions and 9 deletions.
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.
 *