# Near-field table build routing Near-field DuffyRadial tables are built by a batched OpenCL kernel. The scalar per-entry builder that stands behind it is orders of magnitude slower and converges differently at the same requested quadrature orders, so which one ran is both a performance fact and an accuracy fact. Four routings are recorded: `batched` : the device kernel ran. `scalar` : the scalar builder was chosen **up front**, because the batched path was never eligible: the kernel does not support it, or — the case that catches people — `get_table` was called with neither a `queue` nor a `cl_ctx`, both of which are optional in the public manager API. This route logs `[duffy:builder] mode=scalar` at `INFO` and raises no warning. It is not a failure, and it is not loud; if a build is unexpectedly slow, check that a queue reached the manager. `scalar-adaptive` : the scalar builder ran because the build configuration asked for adaptive quadrature. `scalar-fallback` : the batched build was attempted and **raised**. This one is never silent. ## The fallback is loud A `scalar-fallback` — and only that routing — announces itself at the moment it happens: a logged `WARNING`, a `[duffy:builder] mode=scalar-fallback` line, and a `RuntimeWarning` carrying the kernel class, dimension, exception type and reason. The other three routings only log at `INFO`, so a slow build that was never a batched attempt looks exactly like a fast one in a default log. ## Every routing is recorded - The builder records `table.build_routing` — one of the four above for any table built since routing recording exists; a payload cached before that loads as `unknown`, a legacy state that says nothing about how it was built — and, for a fallback, `table.build_fallback_reason`. Both are persisted with the cached payload, so a warm, cache-loaded table still reports how it was originally built (`volumential.opcounters.direct_build_routing`). - `volumential.opcounters.direct_build_routing` is what measurement code reads to emit the routing as a `direct_build_routing` CSV column, which is why it belongs in the record of any table-timing run ({doc}`../benchmarks/index`). A run that does not record it cannot be checked after the fact: read the routing off the table object instead, or run under strict mode. ## Strict mode ```bash export VOLUMENTIAL_DUFFY_NO_FALLBACK=1 ``` turns the fallback into a `RuntimeError` instead. Campaign runs use this so that a table which dropped to the scalar builder cannot end up recorded as a batched build by a run nobody was watching. Any value other than unset, `0`, `false`, `no` or `off` enables strict mode. Strict mode also applies on the **load** path, where the builder never runs. A cached table is refused with an `UnverifiedBuildRoutingError` when its recorded routing is - `scalar-fallback` — it is exactly the differently converged data the switch exists to refuse; - `unknown` — a payload written before routing was recorded, so the build cannot be vouched for; or - anything outside the recognized set (`batched`, `scalar`, `scalar-adaptive`, `scalar-fallback`, `unknown`) — a damaged payload whose routing reads `scalar-fallbac` says nothing about which builder ran, so it is corrupt provenance rather than verified provenance. In every case the error names the table and the remedy: rebuild it with `force_recompute=True`, or unset the switch to accept the cached data. Without this check, a strict campaign whose cache had already been warmed would load and use exactly the data the switch exists to refuse — and, because the routing is faithfully restored from the payload, would report it correctly while doing so. Externally assembled tables (`build_method = ExternalAssembly`, for example a registered windowed RKE assembly; see {doc}`../design-notes/windowed-channels`) are exempt: they were never a DuffyRadial build, which is why the assemblers clear the routing, and their provenance is carried by the build method, the provenance kind and the payload checksum the load path verifies. Strictness is an environment switch rather than a `DuffyBuildConfig` field because the build config is hashed into the table-cache fingerprint, and an operational strictness policy should not invalidate cached numerical data. ## Complex exponentials in the generated quadrature kernel The fused Duffy quadrature kernel rewrites `exp(re + i*im)` into `exp(re) * (cos(im) + i*sin(im))` before code generation *wherever it can prove both halves are real doubles* — see "Why it is guarded" below, and note that an exponent it cannot prove keeps its `cdouble_exp`. The rewrite applies to the standard real-parameter kernels, which then reach the device as real `exp`/`cos`/`sin` calls; a caller-supplied `sumpy_knl` whose phase is not provably real — `HelmholtzKernel(dim, allow_evanescent=True)` is the case the suite pins — deliberately keeps `cdouble_exp` and does not get the speedup below. `pyopencl` implements `cdouble_exp` with the OpenCL `sincos(x, &cosx)` out-parameter builtin, which on the PoCL 7.0 / LLVM 19.1.7 CPU driver costs about 200 ns per call against about 1.6 ns for a separate `sin`/`cos` pair; since the quadrature evaluates the kernel at every Duffy node, that one builtin made the 3D Helmholtz direct table build roughly ten times slower than the otherwise identical Yukawa build. The rewrite is `exp(a+b) = exp(a)exp(b)` with Euler's formula over an exact structural split of the exponent, so it is valid for genuinely complex exponents (the damped `exp((-a + i b) r)` form included) and leaves real exponents untouched. ### Where the 200 ns come from `sincos` itself is not slow. On PoCL's host-CPU device, fp64 `sincos` comes from the bundled kernel library, whose range reduction and polynomial call `fma()` 18 to 25 times per element, and `fma()` there is the correctly rounded fused operation. On an x86-64 CPU **without** an FMA unit — AVX-only, pre-Haswell — the backend cannot lower that to a hardware instruction and emits a call to the C library's *software* `fma`, at about 6.2 ns per element: 18 of them predict 112 ns, against 104 ns measured for `sincos` in a self-contained reproducer. `sin` and `cos` never enter the kernel library at all — they are clang builtins, vectorized and lowered to `libmvec` — and in that same reproducer the pair costs about 0.7 ns, a ratio of about **146x**. Those two numbers are the reproducer's, not the ones quoted above: the ~200 ns and ~1.6 ns figures were measured through the fused Duffy program, where each call sits in the surrounding quadrature kernel, so the absolute costs differ and the ratio there is about 130x. Two harnesses, one conclusion — do not mix a numerator from one with a denominator from the other. On the other CPU OpenCL runtime on the same machine the ratio is about 1x, and on a GPU about 1x as well (both NVIDIA's runtime and PoCL's CUDA device compile `sincos` from the vendor's device library, and `fma` is a single instruction there). So this is specific to PoCL's host-CPU device on targets with no FMA unit, and it is not confined to `sincos`: every fp64 builtin that PoCL implements using `fma()` — `remquo`, `remainder`, `acospi`, `asinpi`, `atanpi`, `atan2pi`, `acosh`, `asinh`, `atanh`, and the `sincos`/`log`/`exp` helpers — pays it on such a host, as does an explicit fp64 `fma()` call in kernel code, which is the same software libcall (about 6 ns per call measured). Ordinary fp64 builtins and arithmetic that do not route through `fma()` (`fabs`, `floor`, `sqrt`, plain `*` and `+`) are unaffected, so fp64 results that touch neither the listed builtins nor an explicit `fma()` need no reclassification. The rewrite above removes the exposure for the standard real-parameter kernels, and only for those — an evanescent Helmholtz or other unprovable phase keeps its `cdouble_exp` and stays exposed on such a host. The general consequence is that a run's metadata has to record the CPU class, and that seconds from a host without hardware FMA are not comparable with seconds from one that has it ({doc}`../benchmarks/index` is the full list of what a run has to record). Tracked in [#138](https://github.com/xywei/volumential/issues/138). ### Why it is guarded The rewrite is exact in value but not in conditioning once the *phase* `im` can itself be complex: for `z = x + i y`, `cos z` and `sin z` both grow like `exp(|y|)/2` while `exp(i z)` decays like `exp(-y)`, so Euler's formula turns a decaying exponential into a cancelling difference of two large terms. It also replaces one `cdouble_exp`, which promotes its whole argument to double, with bare real calls whose precision `loopy` infers from the expression. So the rewrite happens only for a phase, and a magnitude, that are *provably* real doubles, checked node by node. Every leaf must be a real-typed constant at least as wide as a double, a variable the kernel has not left unproven, an arithmetic combination of those, or a call to a function that is real for real arguments; anything else — an unrecognised node type, a `hankel1` call, a `complex128(0j)` that promotes the operation around it, a post-CSE `CommonSubexpression` wrapping any of those — keeps its `cdouble_exp`. An expression made only of constants is refused whatever their Python types, because nothing in it fixes the emitted precision: `loopy` writes the constant real half of `exp(-200 + 1j*k)` as `exp((float) (-200.0f))`, which underflows where `cdouble_exp` kept the finite `exp(-200)`. A kernel argument counts as proven only when its declared dtype is a real floating type at least as wide as a double. Arguments a caller supplies through `extra_kernel_kwarg_types` are checked by the same rule, since they are not in `integral_knl.get_args()`. Complex (the wave number of `HelmholtzKernel(dim, allow_evanescent=True)`), narrow, integer and undeclared dtypes are all unproven. That is deliberately blunt: `loopy`'s constant-dtype inference cannot be reproduced from the expression tree — an integer argument alone narrows the result of a floating builtin, and even a plain `3.0` beside an integer is emitted as `3.0f` — so the guard does not try to model it. No sumpy kernel this table builds has such an argument (Helmholtz's `k` and Yukawa's `lam` are both `float64`), so the rule costs nothing in practice and gives a guarantee instead of an approximation. The global scaling constant is rewritten under the same guard, since it is evaluated inside both quadrature loops, and the split walks the `/1` wrapper `SympyToPymbolicMapper` leaves around it. ### Measured effect At 3D, `q = 3`, source box level 2: the Helmholtz per-(entry × node) cost drops from about 74 ns to about 8 ns, matching the real-valued Yukawa kernel, with table entries agreeing to `3e-16` relative.