Avoid overflow in the 2x2 pivot inverse of the symmetric-indefinite factorizations - #1377
Avoid overflow in the 2x2 pivot inverse of the symmetric-indefinite factorizations#1377rmlarsen wants to merge 4 commits into
Conversation
…actorizations
The Bunch-Kaufman 2x2 pivot block D_k = [[d11, d21],[d21, d22]] is inverted
by dividing through by the off-diagonal d21 so that the scaled determinant
DENOM = d11*d22/d21^2 - 1 stays O(1). The factor entries are
l_j = (1/DENOM) * ((d11/d21)*w0_j - w1_j) / d21 (and its d22 partner),
a per-entry division by d21 and a scaling by T = 1/DENOM. The reference code
does these two operations in an order that overflows at the extremes of the
exponent range, silently in the classic path and with a false INFO in the
rook/RK path:
* The classic unblocked/packed routines (xSYTF2, xHETF2, xSPTRF, xHPTRF) and
the classic panel routines (xLASYF, xLAHEF) hoist the reciprocal as
D21 = T/d21 and then multiply. T/d21 overflows to +-Inf when d21 is
subnormal, so a well-conditioned system with a subnormal off-diagonal is
factored into Inf/NaN and solved to NaN with INFO = 0.
* The rook/RK routines (xSYTF2_ROOK, xHETF2_ROOK, xSYTF2_RK, xHETF2_RK) form
T*(d22*w0_j - w1_j) and divide by d21 afterwards. The intermediate
T*(...) overflows when the block is near the overflow threshold, even though
the factor entry is O(1); the resulting Inf trips the isnan pivot guard and
reports INFO > 0 (a false singularity) for a nonsingular matrix.
Divide each entry by d21 first and scale by T afterwards, in every routine:
l_j = T * (((d11/d21)*w0_j - w1_j) / d21).
The quotient equals DENOM*l_j, which is within a factor 1 + alpha^2 (< 2) of
the result, so it overflows only where l_j itself is within that factor of the
overflow threshold, while the subnormal case is a plain correctly-rounded
division. This is the order xSYTRS/xSYTRS_3 already use to apply D^-1 in the
solve. The hoisted 1/d21 becomes a per-entry division; the extra cost is one
division per factor entry of a 2x2 pivot.
The full LAPACK linear-equation test suite (xlintst{s,d,c,z}) passes with no
new failures across the SY/SR/SK/SA/SP and HE/HR/HK/HA families.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1377 +/- ##
==========================================
+ Coverage 69.01% 69.36% +0.34%
==========================================
Files 6122 6122
Lines 486123 486335 +212
Branches 23286 23268 -18
==========================================
+ Hits 335514 337328 +1814
+ Misses 150420 148569 -1851
- Partials 189 438 +249
... and 146 files with indirect coverage changes Continue to review full report in Codecov by Harness.
|
|
@langou I am going to amend this by one commit to fix a potential overflow regression in gfortran's intrinsic complex division, and add a regression test. Sorry for the churn. |
|
There is a numerical robustness limitation in gfortran’s default complex division. I reproduced it at both -O0 and -O3, without fast-math: with M = HUGE(1d0), (M*(17/32 + 17i/32)) / (M*(-1/2 - i/2)) returns -Inf instead of -1.0625. LAPACK's ZLADIV/CLADIV return the correct results, but produce slower code. Investigating an alternative approach... |
Use scaled complex division for extreme operands. Retain reciprocal multiplication only when the pivot, reciprocal, and numerator bounds exclude intermediate overflow. Add regression coverage for all four precisions, both triangles, packed storage, and blocked panels, including a large numerator with a moderate pivot. Register the tests in both CMake and Make builds.
This reverts commit 8a8021e.
|
@langou working around the issue in gfortran's complex division adds a lot of code complexity and the fix arguably belongs in gfortran. Reverted to the original version you reviewed and added a comment in the PR description. |
| DOUBLE PRECISION ZERO | ||
| PARAMETER ( ZERO = 0.0D+0 ) | ||
| DOUBLE PRECISION FOUR | ||
| PARAMETER ( FOUR = 4.0D+0 ) |
There was a problem hiding this comment.
Ah my favorite kind of named constant: The obvious one.
The matrix generator scales what it produces to the requested norm and xLATB4 asks for at most 1/(SFMIN/EPS), so no type of the SY path holds a 2x2 pivot block whose off-diagonal entry is subnormal, which is what the hoisted reciprocal T/d21 overflows on. Type 11 of xCHKSY takes the type 2 matrix, scales two adjacent rows and columns into the subnormal range, and gives them a 2 by 2 block with d11 = d22 = s and d21 = 4 s: the pivot test then chooses that block, and it sits at the end the factorization starts from, the last two indices for 'U' and the first two for 'L'. The reconstruction residual of xSYTRF is the test; on the parent commit it is a NaN for 24 of the shapes per precision, and a NaN is not .GE. THRESH, so the comparison that prints a failure now tests for one. Such a matrix is at the edge of the representable range, so only the factorization is meaningful: its inverse overflows whatever the pivot inversion does, and the reciprocal of its norm is not a condition estimate. The type therefore sets TRFCON, which is how the path already skips those tests, and skips the condition estimate as well. The complex full-storage path takes the same code with a complex d21, where the division computes |d21|^2 and returns before the reciprocal overflows, so the type is added for the real precisions only; the sweep in the pull request covers the packed and Hermitian routines. The test files declare the new xSCAL calls EXTERNAL: the extended-API build renames only the routines a file declares, so without the declaration xlintsts_64 and xlintstd_64 failed to link against the 64-bit BLAS. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
0791d6e to
a59e306
Compare
|
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. One difference from x86: on master, the DSYSV_RK case of the reproducer returns NaN with INFO = 0 on this machine instead of the false INFO = 3. With the fix both cases give the exact solution. |
2 participants
Disclaimer: This PR was prepared using Claude Code / Fable 5.1.
Description
The Bunch–Kaufman symmetric-indefinite factorizations invert a 2×2 pivot block
by dividing through by the off-diagonal$(-1-\alpha^2,\ \alpha^2-1)$ with $\alpha = (1+\sqrt{17})/8$ . A factor entry is
d21, so the scaled determinantdenom = d11*d22/d21^2 - 1staysO(1). The pivot test keepsdenomin the open intervala division by
d21and a scaling byT = 1/denom. The reference code does these two steps in an order that overflows at the extremes of the exponent range — silently in one family, with a falseINFOin the other:xSYTF2,xHETF2,xSPTRF,xHPTRF) and classic panel (xLASYF,xLAHEF) hoist the reciprocal asD21 = T/d21and then multiply.T/d21overflows to±Infwhend21is subnormal, so a well-conditioned system with a subnormal off-diagonal is factored intoInf/NaNand solved toNaNwithINFO = 0(a silently wrong answer).xSYTF2_ROOK,xHETF2_ROOK,xSYTF2_RK,xHETF2_RK) formT*(d22*w0 - w1)and divide byd21afterwards. The intermediateT*(...)overflows when the block is near the overflow threshold, although the factor entry isO(1); theInftrips theisnanpivot guard and reportsINFO > 0, a false singularity, for a nonsingular matrix.Each family is broken at exactly one end of the range and correct at the other, so neither is caught by the other's inputs.
Fix. Divide each entry by
d21first and scale byTafterwards, uniformly in all 30 routines:The quotient equals
denom*l_j, within a factor1 + alpha^2 < 2ofl_j, so it overflows only wherel_jitself is within that factor of the threshold, while the subnormal case becomes a single correctly-rounded division. This is the orderxSYTRS/xSYTRS_3already use to apply the inverse ofDin the solve. The hoisted1/d21becomes a per-entry division: the added cost is one division per factor entry of a 2×2 pivot.Known gfortran limitation. The new order can expose internal overflow in gfortran 13.3’s complex division even when the quotient is small, at both
-O0and-O3without fast-math (reproducer).CLADIV/ZLADIVavoid this but add overhead; that workaround is not included in this PR.Minimal reproducer
Both systems have exact solution
e3.A 3×3 block embedded at the top (
L) or bottom (U) of a 130×130 identity exercises the blocked panel path (xLASYF/xLAHEF) and shows the same before/after for every driver in the family (?SYSV,?SYSV_ROOK,?SYSV_RK,?SPSV,?HESV,?HESV_ROOK,?HESV_RK,?HPSV, real/complex, Hermitian and complex-symmetric). Atmaster,ZHPSV/ZSPSVon the subnormal block corrupt the heap (double free or corruption); the fix removes that too.Regression test. No type of the
?SYpath could hold a subnormal 2 by 2 pivot block:xLATMSscales what it generates to the normxLATB4asks for, and that is at most1/(SFMIN/EPS). Type 11 takes the type 2 matrix, scales two adjacent rows and columns into the subnormal range, and gives them a block withd11 = d22 = sandd21 = 4 s, so the pivot test chooses it; the block sits at the end the factorization starts from, the last two indices for'U'and the first two for'L'. The reconstruction residual ofxSYTRFis the test, and on the parent commit it is a NaN for 24 shapes per precision. A NaN is not.GE. THRESH, so the comparison that prints a failure now tests for one as well.Such a matrix sits at the edge of the representable range, so only the factorization is meaningful there: its inverse overflows whatever the pivot inversion does, and the reciprocal of its norm is not a condition estimate. The type sets
TRFCON, which is how the path already skips those tests, and skips the condition estimate too. The complex full-storage routines take the same code with a complexd21, where the division forms|d21|^2and returns before the reciprocal overflows, so the type is added for the real precisions; the packed and Hermitian routines are covered by the sweep below.Validation. The full LAPACK linear-equation test suite (
xlintst{s,d,c,z}on{s,d,c,z}test.in) passes with no new failures — every SY/SR/SK/SA/S2/SP and HE/HR/HK/HA family at threshold, e.g.DSY2918,DSR1840,DSK1840,DSP2476 tests. Built and run with GCC 13.3.Checklist
D**(-1)derivation comments inxLASYF/xLAHEFand the inline comments now describe the divide-then-scale order.)