slarfgp/dlarfgp: avoid overflow when computing 1/ALPHA

In the general case (xnorm > eps*|alpha|, beta >= 0), the
Householder reflector formula computes

  ALPHA = -XNORM^2 / (alpha + beta)

which can be subnormal even when beta itself is safely above
SMLNUM.  When this subnormal ALPHA was used in the subsequent

  CALL SSCAL( N-1, ONE / ALPHA, X, INCX )

the reciprocal overflowed (e.g. ALPHA = 2^-137 → ONE/ALPHA = 2^137
which exceeds SP max ~3.4e38).

Fix: guard the SSCAL with ABS(ALPHA) < SMLNUM.  When ALPHA is very
small, scale X by ONE/SMLNUM then by SMLNUM/ALPHA — both safe
since ONE/SMLNUM = BIGNUM is below the overflow threshold and
SMLNUM/ALPHA ≤ MAX_EXPONENT.

The complex variants (clarfgp/zlarfgp) are already safe because
they compute 1/ALPHA via CLADIV/ZLADIV, which internally handles
overflow/underflow without intermediate overflow.

Fixes #938

SRC/dlarfgp.f

@@ -220,7 +220,12 @@
 *
 *           This is the general case.
 *
-            CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
+            IF( ABS( ALPHA ).LT.SMLNUM ) THEN
+               CALL DSCAL( N-1, ONE / SMLNUM, X, INCX )
+               CALL DSCAL( N-1, SMLNUM / ALPHA, X, INCX )
+            ELSE
+               CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
+            END IF
 *
          END IF
 *

SRC/slarfgp.f

@@ -220,7 +220,12 @@
 *
 *           This is the general case.
 *
-            CALL SSCAL( N-1, ONE / ALPHA, X, INCX )
+            IF( ABS( ALPHA ).LT.SMLNUM ) THEN
+               CALL SSCAL( N-1, ONE / SMLNUM, X, INCX )
+               CALL SSCAL( N-1, SMLNUM / ALPHA, X, INCX )
+            ELSE
+               CALL SSCAL( N-1, ONE / ALPHA, X, INCX )
+            END IF
 *
          END IF
 *