C/ZLARFGP must re-scale when ALPHA is not real

[christoph-conrads]

Description

When a vector x with a single nonzero entry x_1 sufficiently small in modulus is passed to CLARFGP or ZLARFGP and if the imaginary part of x_1 is nonzero, then the computed scalar τ is inaccurate. Specifically if either the real of the imaginary part of x_1 (or both) is denormalized, then the the modulus computed by xLARFGP does not possess a small relative error (see SRC/clarfgp.f, line 172 for example).

The problem does not occur with xLARFG (note the missing "P" in the function) name because it scales the vector before computing the modulus of ALPHA when at least one of the following two conditions holds (see SRC/clarfg.f):

• the norm of the vector entries (x_2, x_3, ..., x_n) is nonzero or
• the imaginary part of x_1 is nonzero.

xlarfgP does not check for the second condition.

The issue was found while working on #406. If inaccurate factors τ computed by xlarfgP are passed to xUNGQR, then the matrices calculated by the latter function are far from being unitary (‖U^* U - I‖ is noticeably different from zero).

Checklist

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

[christoph-conrads]

Minimum example (compile with gcc -Wextra -Wall -std=c99 -pedantic -llapack):

#include <complex.h>
#include <stdio.h>

typedef int lapack_int;

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

int main()
{
    lapack_int n = 1;
    float complex x[1] = { 2.073921727e-43f + 3.082856622e-44f * I };
    lapack_int incx = 1;
    float complex tau = -1.0f;

    clarfgp_(&n, x, x + 1, &incx, &tau);

    float complex tau_expected = 1.086842348e-02 - 1.470330722e-01 * I;

    printf("computed tau: %.6e%+.6ej\n", crealf(tau), cimagf(tau));
    printf(
        "expected tau: %.6e%+.6ej\n", crealf(tau_expected), cimagf(tau_expected)
    );
}

Output on my machine:

$ ./a.out 
computed tau: 1.333332e-02-1.466667e-01j
expected tau: 1.086842e-02-1.470331e-01j

The expected value was computed with the help of Python NumPy by the expression

xnorm = np.abs(alpha)
(1 - alpha.real / xnorm) - (alpha.imag / xnorm) * 1j

in double-precision arithmetic.