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
Description
xLARFGP may cause an overflow in line 195. It is my expectation that no overflow takes place.
The operations executed are (
ALPHAis assigned to several times, signs are ignored):The reciprocal of
ALPHA''is computed in line 223 and this computation may overflow when0 < ALPHA ≪ 1and0 < XNORM ≪ ALPHA. For example,ALPHA = ε²andXNORM = ε⁴.The C code demonstrates the problem (compile it with
gcc -Wextra -Wall -std=c99 -pedantic slarfgp.c -llapack -lm):Checklist