A = [a a a], x = [a], where a = DBL_MAX
[0 a a] [0]
[0 0 a] [a]
results in a division of zero by zero. I suppose that the rescaling of
branch 345-348.
// main.c
#include <stdio.h>
#include <float.h>
extern void dlatrs_(char *uplo, char *trans, char *diag, char *normin, int *n,
double *A, int *ldA, double *x, double *scale, double *cnorm, int *info);
int main() {
int n = 3, info = 0;
double A[3*3] = {DBL_MAX, 0.0, 0.0, DBL_MAX, DBL_MAX, 0.0, DBL_MAX, DBL_MAX, DBL_MAX};
double x[3] = {DBL_MAX, 0.0, DBL_MAX}, scale, cnorm[3];
char uplo = 'U', trans = 'N', diag = 'N', normin = 'N';
dlatrs_(&uplo, &trans, &diag, &normin, &n, A, &n, x, &scale, cnorm, &info);
printf("scale = %.6e\n", scale);
return 0;
}
Description
DLATRS with the inputs
returns
scale = -nan. The correct result can be computed withdtrsv,x = [1 -1 1]**T.lapack/SRC/dlatrs.f
Lines 341 to 348 in dcba2e2
TSCALevaluates in 346 toONE / (SMLNUM * Inf) = zerofor the sample matrix andlapack/SRC/dlatrs.f
Line 771 in dcba2e2
results in a division of zero by zero. I suppose that the rescaling of
Aandxand the recomputation ofcnormwent missing in theelsebranch 345-348.Copy & paste reproducer
Checklist