Description

xLARFGP computes a Householder reflector. Given an input vector x of dimension m, xLARFGP computes H = I - β vv^* such that

• H x = ‖x‖ e1 and
• v(1) = 1,

where β is a scalar and v is a vector.

#934 attempted to fix an overflow when 0 ≪ ‖x(2:m)‖ ≪ x(1) ≪ 1 by settings x(2:m) to zero under certain conditions. This does not resolve the issue. For example, let m = 2 and x = [ε^4, 2 ε^5], then an overflow takes place in the following code when 1 / ALPHA is computed:

From the point of view of mathematics, this line should simply compute

• x / (x(1) - ‖x‖) if x(1) < 0,
• -x / [‖x(2:m)‖² · (x(1) + ‖x‖)] otherwise.

(See, e.g., Matrix Computations, 4th edition, Algorithm 5.1.1). Obviously the divisions cannot overflow and now the problem source becomes evident: the computation of the reciprocal of ALPHA. The use of, e.g., SLASCL in the case of INCX = 1 resolves the problem.

Minimal example demonstrating the problem (compile with cc -Wextra -Wall -std=c99 -pedantic -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, 4.0f), 2.0f * powf(eps, 5.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