xLARFGP unnecessarily overflows

[christoph-conrads]

Description

xLARFGP may cause an overflow in line 195. It is my expectation that no overflow takes place.

The operations executed are (ALPHA is assigned to several times, signs are ignored):

ALPHA ≔ X(1)
XNORM ≔ ∥X(1:)∥₂
BETA ≔ ∥X∥₂
ALPHA' ≔ ALPHA + BETA
ALPHA''  ≔ XNORM * (XNORM / ALPHA')

The reciprocal of ALPHA'' is computed in line 223 and this computation may overflow when 0 < ALPHA ≪ 1 and 0 < XNORM ≪ ALPHA. For example, ALPHA = ε² and XNORM = ε⁴.

The C code demonstrates the problem (compile it with gcc -Wextra -Wall -std=c99 -pedantic slarfgp.c -llapack -lm):

#include <math.h>
#include <stddef.h>
#include <stdio.h>

typedef int lapack_int;

void slarfgp_(
    lapack_int* n,
    float* alpha,
    float* x,
    lapack_int* incx,
    float* tau);

#define N 2

int main()
{
    float nan = NAN;
    float eps = ldexpf(1.0f, -23);
    lapack_int n = N;
    float x[N] = { powf(eps, 2.0f), powf(eps, 4.0f) };
    lapack_int incx = 1;
    float tau = nan;

    printf("input=[%.2e; %.2e]\n", x[0], x[1]);

    slarfgp_(&n, x + 0, x + 1, &incx, &tau);

    printf("output=[%.2e; %.2e]\n", x[0], x[1]);
    printf("tau=%.2e\n", tau);

	if(isinf(x[1])) {
		fprintf(stderr, "x[1] is infinite!\n");
	}
}

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue