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.

lapack/SRC/dlatrs.f

Lines 341 to 348 in dcba2e2

	IMAX = IDAMAX( N, CNORM, 1 )
	TMAX = CNORM( IMAX )
	IF( TMAX.LE.BIGNUM ) THEN
	TSCAL = ONE
	ELSE
	TSCAL = ONE / ( SMLNUM*TMAX )
	CALL DSCAL( N, TSCAL, CNORM, 1 )
	END IF

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

lapack/SRC/dlatrs.f

Line 771 in dcba2e2

SCALE = SCALE / TSCAL

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

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

[angsch]

The patch angsch@cf9b925 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.