Divide by a subnormal pivot in the banded LU factorizations instead of scaling by its reciprocal

[rmlarsen]

Disclaimer: This PR was prepared using Claude Code.

Summary

The banded LU factorizations xGBTF2 and xGBTRF (all eight: S/D/C/Z × unblocked/blocked) form the multipliers of column J as xSCAL( KM, ONE / pivot, ... ). When the pivot is subnormal, ONE / pivot overflows to Inf, the column of L becomes Inf/NaN, and the factorization completes with INFO = 0. xGBSV then returns a NaN solution as success for a system that is perfectly well conditioned. The dense LU (xGETF2, xGETRF2) has guarded this since at least LAPACK 3.2 (it is present in the 2008 import of the trunk) by dividing element-wise when |pivot| < SFMIN; this PR adds the same guard to the banded routines. Nothing changes for a pivot of normal magnitude.

Description

xGETF2 computes the column of multipliers as

IF( ABS(A( J, J )) .GE. SFMIN ) THEN
   CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
ELSE
   DO 20 I = 1, M-J
      A( J+I, J ) = A( J+I, J ) / A( J, J )
20 CONTINUE
END IF

and xGETRF2 the same at the base of its recursion. The banded routines have only the first branch:

CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )      ! dgbtf2.f:251
CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), 1 )    ! dgbtrf.f:327

For a pivot below 2^-1022 (double) the reciprocal is not representable, so every multiplier in the column is Inf, or NaN where the entry is zero; the xGER trailing update then spreads that through the rest of the band. Partial pivoting does not help: it picks the largest entry of the column, and when the whole matrix is small every column is. xGTTRF and xGETC2 are unaffected because they divide. The tridiagonal driver xGTSV and the dense driver xGESV both solve the reproducer below to full accuracy.

Fix. The xGETF2 form at both sites: keep the reciprocal scaling when |pivot| >= SFMIN, divide element-wise otherwise. Eight files. SFMIN = xLAMCH('S') is computed once after the quick return, as in xGETF2; neither banded routine referenced xLAMCH before.

Minimal reproducer

A banded, diagonally dominant matrix with every entry near 2^-1030; b = A * ones, so the exact solution is the vector of ones.

program minimal
  implicit none
  integer, parameter :: n = 8, kl = 2, ku = 2, ldab = 2*kl + ku + 1
  double precision :: a(n,n), ab(ldab,n), b(n), x(n), s
  integer :: ipiv(n), info, i, j
  s = 2d0**(-1030)
  a = 0
  do i = 1, n
    a(i,i) = 4*s
    do j = max(1,i-kl), min(n,i+ku)
      if (j /= i) a(i,j) = s
    end do
  end do
  b = sum(a, dim=2)
  do j = 1, n
    do i = max(1,j-ku), min(n,j+kl)
      ab(kl+ku+1+i-j, j) = a(i,j)
    end do
  end do
  x = b
  call dgesv(n, 1, a, n, ipiv, x, n, info)
  print '(a,i0,a,es9.2)', 'DGESV  (dense): info = ', info, '   max|x-1| = ', maxval(abs(x-1))
  x = b
  call dgbsv(n, kl, ku, 1, ab, ldab, ipiv, x, n, info)
  print '(a,i0,a,es9.2)', 'DGBSV  (band) : info = ', info, '   max|x-1| = ', maxval(abs(x-1))
end program

BEFORE (master):
DGESV  (dense): info = 0   max|x-1| =  1.60E-14
DGBSV  (band) : info = 0   max|x-1| =       NaN

AFTER (this branch):
DGESV  (dense): info = 0   max|x-1| =  1.60E-14
DGBSV  (band) : info = 0   max|x-1| =  1.60E-14

Validation

Band versus dense sweep, 2124 cases per build

xGBSV against xGESV on the same non-symmetric, diagonally dominant banded matrix with b = A * ones: four precisions, n in {1, 2, 3, 5, 8, 13, 40, 64, 130}, KL and KU each in {0, 1, 2, 3, 5, 40} (so both the unblocked xGBTF2 path and, for KL = 40 >= NB, the blocked xGBTRF path), and three scalings: normal (s = 1), subnormal (s = 2^-1030 in double, 2^-140 in single), and near the overflow threshold (s = huge/64). Each line records INFO and max|x - 1| for both drivers plus FNV-1a hashes of the banded factor and IPIV.

scale	cases	band finite and < 1e-3, master	band finite and < 1e-3, this branch	dense finite and < 1e-3	band NaN/Inf, master	band NaN/Inf, this branch
normal	708	706	706	706	0	0
near-max	708	702	702	702	0	0
subnormal	708	148	654	654	560	0

On master, 560 of the 708 subnormal cases come back NaN or Inf with INFO = 0; the 148 that survive have KL = 0 or n = 1, where there is no multiplier to scale. On this branch the banded error equals the dense error in every case (the largest ratio err_band / err_dense over all finite pairs is 1.00 in all four precisions), and INFO agrees with the dense driver in 708 of 708. The 54 subnormal cases where both drivers exceed 1e-3 are all single precision: at 2^-140 a REAL carries about 9 significant bits, so that is the input, not the solver.

For the normal and near-max scalings all 1416 lines, hashes included, are byte-identical between master and this branch: the new branch is reached only when |pivot| < SFMIN.

Regression test. The ?GB path gets a matrix type 9: the type 1 matrix scaled in place into the subnormal range, one eighth of the safe minimum, which is below the reciprocal of the overflow threshold in every precision. xLATMS cannot generate such a matrix, since it scales its output to the requested norm, and it returns Inf or NaN entries for a subnormal ANORM on the shapes that go through its Givens chase.

At that scale the stored matrix carries too few bits for a reconstruction residual to mean anything: in single precision an entry near 2^-134 keeps about fifteen, and the ratio of a correct factorization lands near 300. The type therefore checks what the guarded division promises, that the factor is finite, through the max-norm of the factor rather than xGBT01. The condition-number and solve ratios, which would form the reciprocal of a subnormal norm, are skipped for it as they already are for a block size other than the first.

On the parent commit the type fails 798 times per precision, one for each shape and block size, with the factor holding an infinity; with the fix all four precisions pass.

Test suite. The full LAPACK test suite passes on this branch: 215 of 215 CTest entries, 5447193 LAPACK tests and 315872 BLAS tests with 0 numerical errors and 0 other errors, built and run the same way as the parent commit f96546fc9. That includes the ?GB routine and driver families (DGB: 30261 routine tests, 36567 driver tests) and the _64 extended-API variants. Built with GCC 13.3, CMAKE_BUILD_TYPE=Release, BUILD_INDEX64_EXT_API=ON.

Checklist

• The documentation has been updated. (No interface or documented behavior changes; the inline comment "Compute multipliers" is unchanged.)
• If the PR solves a specific issue, it is set to be closed on merge. (No tracking issue; happy to open one.)

[codecov]

Codecov Report

❌ Patch coverage is 95.83333% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 69.36%. Comparing base (a6c6e74) to head (30850dd).
✅ All tests successful. No failed tests found.

Files with missing lines	Patch %	Lines
TESTING/LIN/cchkgb.f	90.90%	1 Missing ⚠️
TESTING/LIN/dchkgb.f	90.90%	1 Missing ⚠️
TESTING/LIN/schkgb.f	90.90%	1 Missing ⚠️
TESTING/LIN/zchkgb.f	90.90%	1 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1379   +/-   ##
=======================================
  Coverage   69.36%   69.36%           
=======================================
  Files        6122     6122           
  Lines      486711   486783   +72     
  Branches    23268    23268           
=======================================
+ Hits       337584   337652   +68     
- Misses     148689   148693    +4     
  Partials      438      438           

Components	Coverage Δ
BLAS	97.94% <ø> (ø)
CBLAS	96.98% <ø> (ø)
LAPACK	82.39% <100.00%> (+<0.01%)	⬆️
LAPACKE	2.17% <ø> (ø)
TMGLIB	55.69% <ø> (ø)
BLAS testing	88.33% <ø> (ø)
CBLAS testing	89.63% <ø> (ø)
LAPACK testing	82.21% <91.66%> (+<0.01%)	⬆️
LAPACKE testing	∅ <ø> (∅)

Files with missing lines	Coverage Δ
SRC/cgbtf2.f	97.95% <100.00%> (+0.23%)	⬆️
SRC/cgbtrf.f	100.00% <100.00%> (ø)
SRC/dgbtf2.f	97.95% <100.00%> (+0.23%)	⬆️
SRC/dgbtrf.f	100.00% <100.00%> (ø)
SRC/sgbtf2.f	97.95% <100.00%> (+0.23%)	⬆️
SRC/sgbtrf.f	100.00% <100.00%> (ø)
SRC/zgbtf2.f	97.95% <100.00%> (+0.23%)	⬆️
SRC/zgbtrf.f	100.00% <100.00%> (ø)
TESTING/LIN/alahd.f	0.00% <ø> (ø)
TESTING/LIN/cchkaa.F	72.47% <100.00%> (ø)
... and 7 more

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a6c6e74...30850dd. Read the comment docs.

[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.