Return INFO = 1 from xLALSD and xBDSDC instead of stopping in XERBLA when the bidiagonal contains a NaN

[rmlarsen]

Disclaimer: This PR was prepared using Claude Code.

Summary

xLALSD (the solver behind xGELSD) and xBDSDC scale the bidiagonal by its largest entry before the divide and conquer phase, guarded only by ORGNRM.EQ.ZERO. xLANST propagates a NaN, and xLASCL rejects CFROM = NaN as an illegal argument, so the call ends in XERBLA, and the reference XERBLA stops the program. A NaN anywhere in A reaches the bidiagonal of xGELSD, and so does an infinity, because the driver's own scaling multiplies by zero and leaves a NaN where the infinity was. This PR tests the norm with xISNAN at those sites and returns INFO = 1, the value xLASDA already uses when a singular value fails to converge. Nothing changes for finite input.

Description

With a NaN at A(2,2) of a 6-by-4 least-squares problem, one process per call:

driver	master
DGELS, DGELST, DGETSLS, DGELSY	returns, INFO = 0, NaN solution
DGELSS	returns, INFO = 3 (DBDSQR did not converge)
DGELSD, ZGELSD	process stops: ** On entry to DLASCL parameter number 4 had an illegal value

The backtrace is DGELSD -> DLALSD -> DLASCL -> XERBLA, from the unconditional DLASCL( 'G', 0, 0, ORGNRM, ONE, ... ) after ORGNRM = DLANST( 'M', N, D, E ). An infinity in A ends the same way: ANRM > BIGNUM makes DGELSD call DLASCL( 'G', 0, 0, ANRM, BIGNUM, ... ), which multiplies by zero and turns the infinity into a NaN. xGESDD guards its own input against a NaN (returns INFO = -4, #469), but an infinity passes that guard, becomes a NaN in the same way, and stops the process inside xBDSDC. Callers that replace XERBLA with an exception see an illegal-argument error from DLASCL for what is a data problem.

Whether a NaN is detected in a scan that uses MAX and comparisons is compiler dependent; the sites changed here use xISNAN, which is not.

Fix. At the scaling sites of {s,d,c,z}lalsd.f (the N = 1 branch and the general one) and {s,d}bdsdc.f, add ELSE IF( xISNAN( ... ) ) THEN INFO = 1; RETURN. The INFO descriptions of xLALSD, xBDSDC and xGELSD document the new return. Ten files, six with code changes.

xBDSDC reaches the changed site only on its divide and conquer path (COMPQ = 'I' or 'P' and N > SMLSIZ). Its small-problem path and its singular-values-only path go through xLASDQ as before and are not changed by this PR: for N <= 2 they return INFO = 0 with a NaN result, for larger N a non-convergence count.

Minimal reproducer

program minimal
  implicit none
  integer, parameter :: m = 6, n = 4
  double precision :: a(m, n), b(m), s(n), work(4096), z
  integer :: iwork(256), info, rank, i, j
  do j = 1, n
    do i = 1, m
      a(i, j) = 1 + mod(3*i + 5*j, 7) / 7d0
      if (i == j) a(i, j) = a(i, j) + 3
    end do
  end do
  z = 0
  a(2, 2) = z / z
  b = 1
  call dgelsd(m, n, 1, a, m, b, m, s, -1d0, rank, work, 4096, iwork, info)
  print '(a,i0)', 'DGELSD returned, info = ', info
end program

BEFORE (master):      ** On entry to DLASCL parameter number  4 had an illegal value
                     (program stopped by XERBLA)
AFTER (this branch):  DGELSD returned, info = 1

Validation

NaN and infinity at every position, all precisions, both sides of the divide and conquer threshold

repro/nan_sweep.f90 runs {s,d,c,z}GELSD on shapes (1,1), (2,2), (5,3), (3,5), (25,25), (26,26), (40,25), (25,40) and (80,60) with a NaN or an infinity at every position (a lattice of 40 positions for the larger shapes), and {s,d}BDSDC with a NaN on every diagonal and off-diagonal entry for the same sizes. XERBLA is replaced by a routine that stops the process with status 99, so a completed run proves that no case reached it. The same program with finite input prints the solutions in hex for a bit-for-bit comparison with master.

routine	cases	returned	INFO = 1
SGELSD, DGELSD, CGELSD, ZGELSD	574 each	all	all
SBDSDC, DBDSDC, N > 25	289 each	all	all
SBDSDC, DBDSDC, N <= 25	116 each	all	pre-existing xLASDQ behavior (see above)

Finite input: 27 (routine, shape) cases, bit-identical between master and this branch. Master stops on the first NaN case.

Regression test. xERRLS gets the case: a 2 by 2 matrix with a NaN in one entry, passed to xGELSD, which must return INFO = 1 without calling XERBLA. It belongs with the error-exit tests rather than in xDRVLS, because the driver returns no solution for such a matrix and the point of the test is that it returns at all; a NaN matrix type in xDRVLS would report the fixed behaviour as a failure (INFO .NE. 0 goes to ALAERH) and let the master behaviour through, since its residual ratios are NaN and RESULT( J ) .GE. THRESH is false for a NaN. On the parent commit the four xerrls sections report *** xGELSD on a matrix with a NaN returned INFO = -4 instead of 1 *** preceded by the xLASCL illegal-argument line, and *** xLS routines failed the tests of the error exits ***; with the fix they pass. xGELSD needs a real workspace of about 800 for a 2 by 2 problem, so the arrays in xERRLS grow accordingly, and the four files join the -Onopropagate list that already keeps the NAG compiler from folding the SQRT( -ONE ) that builds the NaN.

Test suite. The full LAPACK test suite passes on this branch: 215 of 215 CTest entries, 5441901 LAPACK tests and 315872 BLAS tests with 0 numerical errors and 0 other errors, the same totals as the parent commit f96546fc9 built and run the same way (GCC 13.3, CMAKE_BUILD_TYPE=Release, BUILD_INDEX64_EXT_API=ON).

Performance. Timed on a 13th Gen Intel(R) Core(TM) i7-13700HX under WSL2 with the reference BLAS, GCC 13.3, -O2. To separate the change from code-placement effects (which move untouched routines by up to 28% between two separately linked static libraries on this machine), the parent library is a shared object shared by both sides, and each benchmark binary carries its own copy of only the changed routines, parent or branch, which interposes over the library's; everything else is byte-identical. Four rounds in alternating order, one core, one process per run, on an idle machine; medians of the per-round medians, with DPOTRF as an untouched control. The benchmark driver and raw output are available on request.

routine	n	parent	this branch	ratio
DGELSD (us/call)	16	11.4	10.8	0.95
	64	179	172	0.97
	256	9870	9964	1.01
	512	78812	82091	1.04
DBDSDC (us/call)	32	64.7	64.9	1.00
	128	1573	1567	1.00
	512	69940	69707	1.00
	1024	518447	513789	0.99
DPOTRF control (us/call)	64 / 256 / 1024			1.03 / 1.04 / 1.00

No measurable difference; the change is one xISNAN test per call.

Not changed here. xBDSQR with singular vectors loops indefinitely on an infinity in D (reachable through xBDSDC with N <= SMLSIZ); xGELSD cannot pass an infinity down because its scaling turns it into a NaN first. That is a separate defect.

Checklist

• The documentation has been updated. (INFO descriptions of xLALSD, xBDSDC and xGELSD.)
• If the PR solves a specific issue, it is set to be closed on merge. (No tracking issue; happy to open one.)

[codecov]

❌ 1 Tests Failed:

Tests completed	Failed	Passed	Skipped
7893	1	7892	0

View the top 1 failed test(s) by shortest run time

openmp-ubuntu-26_04-arm-shared.BLAS_64.blas::dblat1_64.out (DOUBLE PRECISION Level 1 BLAS routines)

Stack Traces | 0s run time

10 numerical error(s), 2529 test(s) run
                                       FAIL

To view more test analytics, go to the Test Analytics Dashboard
_{📋 Got 3 mins? Take this short survey to help us improve Test Analytics.}

[rmlarsen]

Verified on an Apple M4 (macOS, Homebrew gfortran 16.2, Release build with the CI flags). With this branch merged onto current master, the full test suite passes, the new tests fail without the fix, and the reproducer behaves as described above.