From 85de3e20bc7a2f4334ea54526431dda3e717426f Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sun, 11 Oct 2026 19:28:58 -0500 Subject: [PATCH 1/2] Fix AB13CD MAXIT exit and GAMMA**2 FMA divergence - loop never hit INFO=2 at MAXIT; returned gammal/2 with INFO=0 - gamma*gamma - x contracted to fnmsub; reference rounds GAMMA**2 first - wrapper: init fpeak (garbage on INFO!=0) - 4 seeded Fortran differential cases; 15k random sweep now bitwise Co-Authored-By: Claude Opus 5.5 --- python/wrappers/py_ab.c | 2 +- src/AB/ab13cd.c | 17 ++++++------ .../python/data/ab13cd_fortran_reference.npz | Bin 2478 -> 7510 bytes tests/python/test_ab13cd_fortran_reference.py | 26 ++++++++++++------ 4 files changed, 28 insertions(+), 17 deletions(-) diff --git a/python/wrappers/py_ab.c b/python/wrappers/py_ab.c index 3cb4d56..f7b9366 100644 --- a/python/wrappers/py_ab.c +++ b/python/wrappers/py_ab.c @@ -8227,7 +8227,7 @@ PyObject* py_ab13cd(PyObject* self, PyObject* args) { f64 *c_data = (f64*)PyArray_DATA(c_array); f64 *d_data = (f64*)PyArray_DATA(d_array); - f64 fpeak; + f64 fpeak = 0.0; i32 info; f64 hnorm = ab13cd(n, m, np, a_data, lda, b_data, ldb, c_data, ldc, d_data, ldd, diff --git a/src/AB/ab13cd.c b/src/AB/ab13cd.c index 6a0d8e5..7515dcf 100644 --- a/src/AB/ab13cd.c +++ b/src/AB/ab13cd.c @@ -246,15 +246,21 @@ f64 ab13cd(i32 n, i32 m, i32 np, const f64 *a, i32 lda, const f64 *b, i32 ldb, f64 gamma = zero; f64 gammau = zero; - while (iter < MAXIT) { + for (;;) { iter++; + if (iter > MAXIT) { + *info = 2; + return zero; + } gamma = (one + two * tol) * gammal; + // Separate statement blocks FMA contraction; gfortran -O3 reference hoists GAMMA**2 and rounds it. + f64 gamma2 = gamma * gamma; for (i32 j = 0; j < m; j++) { for (i32 i = 0; i <= j; i++) { dwork[iw4_ + i + j * m] = -dwork[iw2 + i + j * m]; } - dwork[iw4_ + j + j * m] = gamma * gamma - dwork[iw2 + j + j * m]; + dwork[iw4_ + j + j * m] = gamma2 - dwork[iw2 + j + j * m]; } SLC_DLACPY("Full", &m, &n, dwork, &m, &dwork[iw5_], &m); @@ -276,7 +282,7 @@ f64 ab13cd(i32 n, i32 m, i32 np, const f64 *a, i32 lda, const f64 *b, i32 ldb, for (i32 i = 0; i <= j; i++) { dwork[iw7 + i + j * np] = -dwork[iw3_ + i + j * np]; } - dwork[iw7 + j + j * np] = gamma * gamma - dwork[iw3_ + j + j * np]; + dwork[iw7 + j + j * np] = gamma2 - dwork[iw3_ + j + j * np]; } SLC_DLACPY("Full", &np, &n, c, &ldc, &dwork[iw8], &np); @@ -384,11 +390,6 @@ f64 ab13cd(i32 n, i32 m, i32 np, const f64 *a, i32 lda, const f64 *b, i32 ldb, } } - if (iter > MAXIT) { - *info = 2; - return zero; - } - f64 hnorm = (gammal + gammau) / two; dwork[0] = (f64)lwamax; diff --git a/tests/python/data/ab13cd_fortran_reference.npz b/tests/python/data/ab13cd_fortran_reference.npz index bed25cfe443cf108f1f56315f973ffd3a5f22097..b75e1be27f9525cfed1e750d2cca036e80fae0df 100644 GIT binary patch literal 7510 zcmd5=3sh4_8cqNOBUQjh@c~A=Y*2%hhiDO5nPS|hr?C- zsO$GuB(jwLmim@2NYYOW?Iom-E{#3cyZL@~P4{f2P5d%$Rq7FqjeybGJI1 z*s{td36pI$iMVsKYO~S-HESV@$5~-GlSqvR2Kf1h21bY^;vl(HA_|g8q%U-iHzu~M zAXy3dRxCHDnqjqR11k~=I?}`6C+v1Oua9x4qGj$GU-j>yplY0;`JgZsHXj#;83ez77+N;suOUUBU&g*uQvwwgmld7XF58na_t!J^hS|3eB()eM z-VzWh3JCU0s>Mdv>l2%<)Vctc`49EklYCh~G%dy?Kx+S_U?D~7X*jz(RLX`^k-%<`f)2ZHo! zwDuApsSj|Cq27DfwU-58NJgaAB!Oaigq27ZB=r-&zG92ndy|)D(h0+Bx+ABXwgSV1 z)KP~vbU^!^wV9`P_tJJ^h3!A=2gImR8{UH0vf7H|iNf}#dlcW66#5{u#{9? z6=%xtVl3#*=qC;e5?hHBvan!DK(OBnX=u;^h4fS$hT%!O7eNDmUJc>&gy{D z$=^(s|M3JKQ+%X-`=9OLsN_$tr4O!wlOLL!?8(+6ZXMe7=<|b(AhgBBbjkX15OHe# z6JxC&AX(QTNj}|9Pi?{%W`oh8QDKsr4s3@*N52b!; zq`}M2())5Sc(1KC9a0+8n@?BJHskVre|lp$#xJxfk;_6AzI}D(3wfrgTRGScCt-AX z+6wu$Y~{JEr+#h5SNeqa$KTri{cgq9I@piGipu7goZMN?zd;O#$&) z6McS(Yy!#kF1yOhpFvl1pJNL+dW2!V%gpR|+JV-z1HsbRcA$FSZ|ajf7vaI)uBM-I znm}&ud9$mkTrk(uyY1xD4%nF2Tp77kpQ!${x=LVl12(KK%Uf}^9;`V&vNXrG4*uL- zm1(i|H^?(t?;W=BG2Hyf?SN=k1Bfl(`{RRqTw-z3<$aSapThXvie`CN5m31X<(4Hp z26dOLOzL;sg$pWbBBtNC0J0o~;W{l3z?X}TyKPI>p_SQ}SEM)&k!j7zb5DKpV9PSo zEGtny&%GpZ-C5ImFuqJ#I9pN=XEtRP1{XF%(>P=AIa=xP5bvhmsGL0Te2ve$dn_+P zhd(~cvR~8%Q>rI#c8ckQ_av)xT=%NL_H)st`VG}^L(!}X!tMbC8!TK~-|eMIoLAHo zJ_J*mLku&7EgPoDLW#8|HeSg7?i~)7*fZjGqB|@A;|m_HGEp^vT!mz=;?!-hujJCF z9k&|b8spLXrQ!Esyul}N8+^IM6Q@sd_N;3MJEzO+-}|Nk@K60-=Nfq$?($Tm$Mz6# zikYe3=hK>KOcnV1&KQCz&3Tx{`DN2`B`zHLemfx4BbU__I;iV52NL-gg39s+d;@hT>fVen9LYsnA;_6zA@Q1=Zz zG0>!jDbZmL!^jEhbfu9Y(4>YrK7u(6BQL0n#JrlXM>2P&Xuv^g)yAzEY3s&d3GoG++=^ciT}^cSa7-tDf0?6p!l8h#kG4@#2k00C*_5 zRiu0W_cR?mdY$5*r}`74Y5t7p(Myoo|0{i(KO=nf>cY~X{=WYkG=E0?=%sUo1dS$D zS~g=i%u6xcqq9yUWDjZXVU3FRm(ln6-{cCv*R=i`R^!yGJ@>4hYzZ}kg068cEcmMzZ delta 33 jcmca+wN7|L5&Py}Vj4^$tZWQGzyX9BObiSK>>wThlY<4g diff --git a/tests/python/test_ab13cd_fortran_reference.py b/tests/python/test_ab13cd_fortran_reference.py index 02ebf13..bdf175c 100644 --- a/tests/python/test_ab13cd_fortran_reference.py +++ b/tests/python/test_ab13cd_fortran_reference.py @@ -29,16 +29,24 @@ def _ab13cd_driver(a, b, c, d) -> str: {fortran_matrix("a", a)}{fortran_matrix("b", b)}{fortran_matrix("c", c)}{fortran_matrix("d", d)} hnorm = AB13CD(n, m, p, a, n, b, n, c, p, d, p, 1.0d-9, iwork, dwork, & ldwork, cwork, lcwork, bwork, info) - print '(I0,1X,ES26.17E3)', info, hnorm + print '(I0,1X,ES26.17E3,1X,ES26.17E3)', info, hnorm, dwork(2) end program main """ DATA = Path(__file__).parent / "data" / "ab13cd_fortran_reference.npz" -# Pre-DLAPY2 C (hypot) returned exactly half the norm for half_norm_stall and -# was off by 1.5% for lightly_damped_mode. -CASES = ["half_norm_stall", "lightly_damped_mode"] +# half_norm_stall, lightly_damped_mode: hypot (vs DLAPY2) changed the initial +# OMEGA; half_norm_stall then hit the MAXIT exit (GAMMAU=0 -> hnorm/2). +# maxit_info2: Fortran exhausts MAXIT -> INFO=2; C fell out of the loop with +# GAMMAU=0 and INFO=0. +# dlapy2_omega: hypot picked a different initial OMEGA. +# gamma2_rounding, gamma2_dpotrf_info: C contracted GAMMA**2 - X into an FMA; +# Fortran rounds GAMMA**2 first. Changed the Hamiltonian / DPOTRF outcome. +CASES = [ + "half_norm_stall", "lightly_damped_mode", "maxit_info2", "dlapy2_omega", + "gamma2_rounding", "gamma2_dpotrf_info", +] @pytest.mark.parametrize("name", CASES) @@ -47,14 +55,16 @@ def test_ab13cd_hnorm_matches_fortran(name, tmp_path): with np.load(DATA) as data: a, b, c, d = (data[f"{name}_{k}"] for k in "abcd") - info_f, hnorm_f = run_fortran_driver(_ab13cd_driver(a, b, c, d), tmp_path).split() + info_f, hnorm_f, fpeak_f = run_fortran_driver(_ab13cd_driver(a, b, c, d), tmp_path).split() n, m = b.shape p = c.shape[0] - hnorm, _, info = ab13cd( + hnorm, fpeak, info = ab13cd( n, m, p, np.asfortranarray(a), np.asfortranarray(b), np.asfortranarray(c), np.asfortranarray(d), 1e-9, ) - assert info == int(info_f) == 0 - np.testing.assert_allclose(hnorm, float(hnorm_f), rtol=1e-10, atol=0) + assert info == int(info_f) + if info == 0: + np.testing.assert_allclose(hnorm, float(hnorm_f), rtol=1e-10, atol=0) + np.testing.assert_allclose(fpeak, float(fpeak_f), rtol=1e-10, atol=0) From 4e25905c0071fb48f3e8489afe776a6af9061640 Mon Sep 17 00:00:00 2001 From: James Joseph Date: Sun, 11 Oct 2026 19:30:33 -0500 Subject: [PATCH 2/2] ergo: file HHJNJ5 (AB13DD FMA divergence) Co-Authored-By: Claude Opus 5.5 --- .ergo/journal.jsonl | 1 + .ergo/plans.jsonl | 1 + 2 files changed, 2 insertions(+) create mode 100644 .ergo/journal.jsonl diff --git a/.ergo/journal.jsonl b/.ergo/journal.jsonl new file mode 100644 index 0000000..69cebaf --- /dev/null +++ b/.ergo/journal.jsonl @@ -0,0 +1 @@ +{"version":1,"task_id":"HHJNJ5","kind":"created","at":"2026-10-12T00:23:24.976736Z"} diff --git a/.ergo/plans.jsonl b/.ergo/plans.jsonl index fc55958..7bad173 100644 --- a/.ergo/plans.jsonl +++ b/.ergo/plans.jsonl @@ -360,3 +360,4 @@ {"type":"result","ts":"2026-05-11T00:15:38.967907Z","data":{"task_id":"HKA3XF","summary":"Final bb03ad/bb04ad report shows no row more than 1.5x slower than Fortran.","path":"benchmark_c_vs_fortran_report.md","sha256_at_attach":"be1164009b3b0151c7c9b584b48c0f98a9ec1d6cb36fac98f76a659735246aba","mtime_at_attach":"2026-05-11T00:15:20.784333433Z","git_commit_at_attach":"cca77d92b979c0d4d3a5fbe330a576cdeffcd834","ts":"2026-05-11T00:15:38.967907Z"}} {"type":"body","ts":"2026-05-11T00:15:38.97072Z","data":{"id":"HKA3XF","body":"## Goal\n- Run the final benchmark comparison and decide whether the C11 Lyapunov generator speed issue is closed.\n- Collect validation evidence for implementation review.\n\n## Result\n- Reinstalled the package after optimization and reran targeted BB03AD/BB04AD behavior tests.\n- Regenerated the combined `bb03ad bb04ad` C-vs-Fortran report after the BB04AD optimization and report-section work.\n- The generated `## Fortran Faster Rows` section shows only BB04AD ex4.4 rows at 1.11x, 1.07x, and 1.06x Fortran advantage, so no row remains more than 1.5x slower than Fortran and no follow-up task is required by this gate.\n\n## Validation Gates\n- `uv pip install .`\n- `.venv/bin/pytest tests/python/test_bb03ad.py tests/python/test_bb04ad.py -v` (30 passed)\n- `.venv/bin/meson compile -C build-bench-c`\n- `.venv/bin/python scripts/benchmark_c_vs_fortran.py --c-build-dir build-bench-c --output benchmark_c_vs_fortran_report.md -r bb03ad bb04ad`\n- `git status --short` checked; only `.ergo/plans.jsonl` and generated `benchmark_c_vs_fortran_report.md` remain outside code commits.\n","ts":"2026-05-11T00:15:38.97072Z"}} {"type":"state","ts":"2026-05-11T00:15:38.97072Z","data":{"id":"HKA3XF","state":"done","ts":"2026-05-11T00:15:38.97072Z"}} +{"type":"transaction","version":1,"events":[{"type":"new_task","ts":"2026-10-12T00:23:24.976736Z","data":{"id":"HHJNJ5","uuid":"3810e153-8d5f-459f-8cd2-b66cd2969df2","epic_id":"","state":"todo","title":"AB13DD: FMA contraction diverges from Fortran reference","body":"Repro: clang -O2 (fp-contract=on) emits fmadd/fnmsub at src/AB/ab13dd.c lines ~213,405,428,454,581-587,635-639,984-985,1079-1080 (objdump -d -l on ab13dd.o); gfortran -O3 AB13DD.f.o has only 8 fused ops. Same class as AB13CD bug: C `gamma*gamma - x` contracted while reference hoists/rounds GAMMA**2 (AB13CD fixed on claude/intelligent-elgamal-797628 via separate `gamma2` statement; random sweep 3000 cases went 41 -\u003e 0 mismatches).\nNext: map each C fused site to the Fortran object's fused/unfused ops; split statements where reference rounds; random differential sweep of ab13dd vs Fortran (driver pattern: tests/python/test_ab13cd_fortran_reference.py).\nFiles: src/AB/ab13dd.c, SLICOT-Reference/src/AB13DD.f\n","created_at":"2026-10-12T00:23:24.976736Z"}}]}