LATRS can reach 0/0

[angsch]

Description

DLATRS with the inputs

A = [a a a], x = [a], where a = DBL_MAX
    [0 a a]      [0]
    [0 0 a]      [a]

returns scale = -nan. The correct result can be computed with dtrsv, x = [1 -1 1]**T.

https://github.com/Reference-LAPACK/lapack/blob/dcba2e274af902f82ff554c47fc3178dbf9ac3e0/SRC/dlatrs.f#L341-L348

TSCAL evaluates in 346 to ONE / (SMLNUM * Inf) = zero for the sample matrix and

https://github.com/Reference-LAPACK/lapack/blob/dcba2e274af902f82ff554c47fc3178dbf9ac3e0/SRC/dlatrs.f#L771 results in a division of zero by zero. I suppose that the rescaling of A and x and the recomputation of cnorm went missing in the else branch 345-348.

Copy & paste reproducer

// 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;
}

Checklist

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

[angsch]

The patch https://github.com/angsch/lapack/commit/cf9b925741b8286a5402ef333b772f1a50a93bd6 should fill the gap. I am confident that LATRS can now really solve a strict superset of TRSV without any avoidable infinities and NaNs. Could someone please have a look? Then I would do the other data types as well.

For the example above, scale = 0.5, x = [0.5 -0.5 0.5]**T is a valid representation of the correct solution as of 1/scale * x = [1 -1 1]**T.