Conversation
SCASUM spelled out ABS(REAL(..)) + ABS(AIMAG(..)) inline while its double counterpart called the helper, which left SCABS1 without a single caller in the library and so untested. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1418 +/- ##
=======================================
Coverage 69.36% 69.36%
=======================================
Files 6122 6122
Lines 486711 486711
Branches 23268 23268
=======================================
+ Hits 337584 337587 +3
+ Misses 148689 148686 -3
Partials 438 438
... and 1 file with indirect coverage changes Continue to review full report in Codecov by Harness.
|
|
I think this is the documented behavior and hence breaking contract and possibly backwards compatinility. |
|
True, it would change the results. I ran some tests to see if there is some consensus among the vendors, but it seems there isn't: n = 1000, 2000 random vectors:
So I'm not sure whether this is documented behavior, other than it being like that for over 20 years. I'd still like the two routines to operate similarly, so we could replace |
|
I happen to translate all Fortran code in SciPy and I am feeling like this also relates to the overflow concerns because I saw this definition a bit too often to attribute to mistakes. I wonder if there is any historical reasoning why this 1- norm is adopted instead of 2- for complex number magnitude since norms are equivalent if the value is not important but ordering is I can "sense" the rationale behind it. |
|
Closing in favor of #1421 |
2 participants
Description
DZASUMaccumulates withSTEMP = STEMP + DCABS1(ZX(I)).SCASUM, its single-precision counterpart, spells the same quantity out inline asSTEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))).SCABS1computes exactlyABS(REAL(Z)) + ABS(AIMAG(Z)), so the two files were doing the same arithmetic by different routes for no stated reason.The consequence is that
SCABS1has no caller anywhere in the Fortran library. The only references to it on master areBLAS/SRC/DEPRECATED/icamax.f, which no build file lists, and the CBLAS wrapperscblas_scabs1.candscabs1sub.f, which no test calls.DCABS1is reached throughDZASUM. This PR makesSCASUMcall the helper, which brings the two precisions into line and givesSCABS1a caller.Coverage
Measured locally with gfortran 13.3.0:
BLAS/SRCbuilt-O0 -g --coverage -fno-inline, the four?blat1drivers linked against it and run, thengcov -b -n.scabs1.fdcabs1.fscasum.fitself is unmoved at 100% of 12 lines and 71.43% of 14 branches, which is whatdzasum.falso scores.This is not bit-for-bit identical
The inline form evaluates
(STEMP + |re|) + |im|, two roundings into the accumulator per element. The helper form evaluatesSTEMP + (|re| + |im|), one. Over 20000 random vectors of length 1000, 19376 gave a different single-precision result.The new grouping is slightly the more accurate of the two, which is what the grouping predicts: adding the two same-magnitude parts to each other first, then once into a larger running total, rounds less than folding each part in separately. Against a
DZASUMreference over 5000 random vectors of length 1000:SCABS1(after)For reference, single-precision epsilon is 1.192e-07.
SCASUMnow also rounds the same wayDZASUMalways has, so the two precisions agree on more than just the expression.This matters because
SCASUMis not only a user-facing routine:CLAEIN,CLATRS,CLATPS,CLATBS,CLATDF,CTREVCandCTREVC3call it, as do several EIG test helpers. The suite results below cover that.Cost
On gfortran 13.3.0, x86-64, n = 100000, timing
SCASUMthree ways:-O2-O2 -fltoSCABS1, this PRTwo separate effects are at work and they pull opposite ways.
The call itself costs about 12%, because without link-time optimization it cannot be inlined across translation units. That is precisely what
DZASUMhas always paid: substituting the inline expression back intodzasum.fand timing the pair the same way gives 1.13x. On that axis this change bringsSCASUMto parity withDZASUMrather than introducing a new cost.The regrouping is worth a factor of two, and it is the grouping rather than the call or the inlining that earns it.
STEMP + (|re| + |im|)puts one dependent addition on the critical path whereSTEMP + |re| + |im|puts two, which gfortran emits as oneaddssinto the accumulator instead of two. Parenthesising the inline expression, with no call at all, reaches the same 0.060s at plain-O2and produces bit-identical results to this PR. So the 2x seen under LTO is the grouping becoming visible once the call is inlined away, not a benefit of LTO as such.The net is that this PR is about 12% slower than master at
-O2and about twice as fast with LTO.DZASUMsits in the same place for the same reason. Only gfortran on x86-64 was timed.