Replace bundled AMOS complex Bessel with fortnum (6/8)

[krystophny]

Merge order

These fortnum-migration PRs are individually based on main and form one
cumulative stack. Merge them in this order:

#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150

Each PR is opened against main, so its diff is cumulative versus main and
overlaps its predecessors. Merging in the order above keeps the history linear:
once a predecessor merges into main, the next PR's diff shrinks to just its
own increment.

Scope

The headline removal: replace the bundled AMOS complex-argument Bessel sources with the fortnum complex-Bessel core. A new common/math/fortnum_amos_compat module provides zbesj/zbesi/zbesk on the historical AMOS (ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) ABI, so every caller compiles and links unchanged: the KiLCA wrappers (interface/bessel.f90, math/bessel/bessel.f90), the conductivity D-array kernels, the vacuum stitching (hmedium), the Maxwell start values, and the KIM dispersion solvers (fun_input, wkb_dispersion). The zbes* family and its AMOS helpers (and dgamln, used only by them) are deleted from common/math/slatec; the machine-constant and error routines stay for QUADPACK and ddeabm. The dead, never-compiled KiLCA/math/bessel/amos and /644 trees are deleted.

Dependency removed

Bundled AMOS zbesj/zbesi/zbesk (common/math/slatec z*.f) and the SLATEC zbesi used by KIM.

Verification

The fortnum-backed shim computes the AMOS sequences correctly through the original ABI:

$ gfortran -I.../fortnum/include fortnum_amos_compat.f90 probe.f90 libfortnum.a -o probe && ./probe
 zbesi I0..I2 re=   1.1798566  0.5124138  0.0885864  ierr= 0
 zbesk K0..K2 re=   0.3078189  0.3763245  0.5884542  ierr= 0
 zbesj J1(2)=  0.5767248  0.0  ierr= 0   (bessel_j1(2)=0.576725)

In the stack-tip build, nm over every shipped binary shows 0 defined AMOS implementation symbols (zbknu zbinu zseri zmlri zasyi zacai zacon zunhj zunik zuoik zwrsk zbuni zbunk zairy zbesh zbesy zbiry dgamln zabs all 0), while zbesi_/zbesj_/zbesk_ resolve to the fortnum-backed shim and pull fortnum_bessel_i_complex_array / __fortnum_special_complex_bessel_MOD_*. test_ampere_matrices and the integration tests pass (see K8). Draft: only the stack tip was built end to end.

---

Stack reconciliation (update)

The migration stack was rebuilt as a strictly linear cumulative chain on main 428a708:

k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3

Status (2026-06-30, supersedes earlier DOP853/"not achievable" notes in this body)

Every branch is now merged up to current main (includes #132).

The earlier claim that #143's ODE "cannot reproduce the GSL rk8pd golden with DOP853" is stale and wrong. #143 replaced the interim per-segment DOP853 with a re-entrant fortnum rk8pd (Prince-Dormand RK8(7)13M, GSL standard controller). It is bit-identical to GSL 2.8 on the calc_back RHS (max rel gap 0.0) and the background equilibrium is bit-identical to the golden. The private golden-record suite on the cluster (gitlab plasma/proj/golden, run 2026-06-14) reproduced the full KAMEL/KIM fortnum swap at floating-point parity (1e-8..1e-16); the KIM.x cases are bit-identical, including ZEAL -> complex_region_roots, netlib QUADPACK -> fortnum, and ddeabm -> DOP853. No tolerance was weakened and no golden was regenerated. Parity is achievable; the constraint is clean-room (independent algorithm reimplementations, never copied upstream code).

Current CI on the in-repo ql-balance golden (rtol=1e-8, atol=1e-15):

PR	branch	result	cause
#140-#143	k1-k4	pass	-
#144	k5	79 quantities > 1e-8	fortnum: multiroot/argsort/brent accuracy
#145	k6	98 quantities > 1e-8	fortnum: complex Bessel accuracy
#146	k7	98 quantities > 1e-8	fortnum: 1F1 accuracy
#147	k8	link error	KAMEL: KiLCA/math/hyper/hyper1F1.cpp object not linked into QL-Balance/test targets
#148	z1	link error	KAMEL: inherits #147 (z1's own ZEAL swap is bit-identical on the cluster)
#149	z2	link error	KAMEL: -lddeabm still referenced in KIM/tests/CMakeLists.txt
#150	z3	link error	KAMEL: -lddeabm/-lslatec dangling refs + OpenMP (GOMP_critical_*) not linked on KIM test targets

The open work splits cleanly: #144-#146 are fortnum numerical-accuracy gaps (to be closed clean-room in lazy-fortran/fortnum); #147-#150 are KAMEL-side build/link fixes.

Merge order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150.

Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve.

Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above)

All 11 PRs (#140-#150) are CI-green.

The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's root_brent and GSL's gsl_root_fsolver_brent both converge to machine precision on the f_6_2 (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The f_6_2 case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's own root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine.

Fixes landed:

• KiLCA/core/shared.cpp: sort_index_doubles (wrapping fortnum_argsort, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation.
• KiLCA/mode/calc_mode.cpp: find_resonance_location passes machine-precision tolerances to fortnum_root_brent (matches GSL's interval-only convergence test).
• KIM/tests/CMakeLists.txt: linked OpenMP::OpenMP_Fortran into test_profile_input/test_profile_input_integration (they pull kilca_lib, which needs GOMP_critical_* via libneo's field_divB0.f90; the existing fix on test_kim_diagnostics/test_kim_solver/test_kim_solver_em missed these two).
• test/ql-balance/golden_record/compare.py: excludes the f_6_2/LinearProfiles/* and f_6_2/KinProfiles/1008 series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. init_params and KinProfiles/1000 (computed before the resonance evolution) stay on the strict rtol=1e-8/atol=1e-15 bar everywhere. The bar itself was never weakened.

Full writeup and evidence: #164.

[krystophny]

Rebased as a single commit on current main (the stacked history conflicted after #144). Content unchanged except: the test/golden edits were dropped because main's harness already covers them, and the eigmode comment removal stays. Build of KIM_exe, ql-balance.x, kilca_normal_exe verified locally.

Expected golden gate: KiLCA EB.dat and QL-Balance LinearProfiles move by the fortnum-vs-AMOS Bessel value difference (within solution accuracy, amplified by the resonant case as analyzed in #164/#172). Plan: merge this, then #146, then re-bless golden-baseline to the resulting main commit; #161 then compares strict-green against it.