From 1c9eb577fdf76f521f1d9f08c5a32b2672c9f70a Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 17:33:53 +0000 Subject: [PATCH 01/14] doc: add optimization & bug-hazard review notes Code review of the hot kernels and minimizer drivers: latent red-black race for odd n in rosso_krauth, redundant objective evaluation and reduction/sync overhead in lbfgs, duplicate edge sampling in gradient(), RNG pool contention in dynamics, plus build-config and convention hazards (line-tension mapping, default build type). https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- docs/optimization_notes.md | 153 +++++++++++++++++++++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 docs/optimization_notes.md diff --git a/docs/optimization_notes.md b/docs/optimization_notes.md new file mode 100644 index 0000000..d96fd46 --- /dev/null +++ b/docs/optimization_notes.md @@ -0,0 +1,153 @@ +# Optimization potential & known issues + +Technical notes from a code review of the qEW port (June 2026). Items are +ordered by expected impact. File references point at the current `main`-line +sources. Nothing here changes physics; everything preserves the documented +analytic consistency between `objective()` and `gradient()` and the golden-file +contracts. + +The hot path, per the benchmark (`executables/benchmark`), is the bicubic noise +gather inside `gradient()` / `objective()` — memory-bandwidth-bound (16 +coefficient loads + ~70 FLOPs per `DeviceNoise::sample()`). Most of the wins +below attack either (a) the number of `sample()` calls, or (b) per-iteration +kernel-launch / device-sync overhead in the minimizers. + +## Bugs / correctness hazards + +1. **Red–black sweep races for odd `n`** — `src/rosso_krauth.h:147` + (`if ((i & 1) != color) return;`). On a periodic ring with odd `n`, sites + `0` and `n-1` have the same parity but are nearest neighbours, so the same + colour pass updates both while each reads the other's `h`. All current + configs use even `nx`, so this is latent — but nothing enforces it. + *Fix:* assert `n % 2 == 0` in `rk_block()` (cheap), or fall back to a + 3-colour sweep for odd `n`. + +2. **Inconsistent pinning-length → line-tension map across run drivers** — + `src/qew_runs.h`: `static_sweep()` (line 67) and `depinning_ramp()` + (line 201) use `line_tension = pinning_length^1.5 * amplitude`, while + `dynamic_anneal()` (line 134) uses `(pinning_length / xi)^1.5 * amplitude`. + The depinning `main` comments that its `pinning_length` is already + `lambda_p/xi`, so each driver may individually match its Python source, but + the same struct-field name means two different physical quantities in + different drivers. *Action:* verify each against `qew_*.py` and either + unify the convention or rename the fields (`pinning_length_over_xi`). No + test currently pins this mapping. + +3. **Floating-point loop counter in the force ramp** — + `src/qew_runs.h:214`: `for (real_t f = 0; f <= c.f_max; f += c.f_step)`. + Exact for the default `f_step = 1.0`, but for a non-representable step + (e.g. `0.1`) accumulated rounding can drop or shift the final step and + perturb the HDF5 key strings (`driving_force=%.4g`). *Fix:* integer index, + `f = k * f_step`. + +4. **No default `CMAKE_BUILD_TYPE`** — root `CMakeLists.txt` never sets a + fallback, so configuring without `-DCMAKE_BUILD_TYPE=Release` silently + produces an unoptimized (`-O0`) build of a bandwidth-bound HPC code. *Fix:* + default to `Release` when unset (warn). Cosmetic: the project is still + named `my-fantastic-simulation-code` (line 5), and `hello.cpp` / + `kokkos_smoketest.cpp` scaffolding is still compiled into the production + `lib` target. + +## Optimization potential + +### Minimizer (L-BFGS, `src/lbfgs.h`) + +5. **Redundant objective evaluation every iteration** — line 200: + `phi0 = objective(h, ...)` is recomputed at the start of each line search, + but `h` is exactly the point accepted by the previous iteration's line + search, whose `phi` is already known. Carrying the accepted value forward + removes one full objective (one `sample()` per site + reduction) per + iteration — when the first trial step is accepted (the common, + preconditioned case) this *halves* the objective-evaluation cost. + +6. **Fuse the trial step into the objective reduce** — each Armijo trial does + a `parallel_for` (write `h = x0 + a*d`, line 212) plus a `parallel_reduce` + (objective). Evaluating the objective directly at `x0(i) + a*d(i)` inside + the reduce, committing `h` only on acceptance, saves one kernel launch and + `n` stores per trial and removes the `deep_copy(h, x0)` undo path + (line 224) and the `deep_copy(x0, h)` save (line 199). + +7. **Fewer reductions / device syncs per iteration** — + - `sy`, `ss`, `yy` (lines 245–247) are three separate `parallel_reduce`s + over the same two vectors; fuse into one reduce with a 3-component + reducer. The `lbfgs_sy` `parallel_for` (line 240) can be folded in too. + - The two-loop recursion issues `2m` blocking dot products per iteration, + each a device→host sync; for the production `n = 4096` a GPU run is + latency- rather than bandwidth-limited here. The `axpy` can be fused + with the next loop's dot (`q -= a*y` and `s·q` in one kernel). + - History update: `deep_copy(S[row], s_tmp)` / `deep_copy(Y[row], y_tmp)` + (lines 256–257) copy two full vectors per iteration; `std::swap`-ing the + `View`s into the slot is free. (The FFT plans are bound to `q`/`d`, not + `s_tmp`/`y_tmp`, so swapping is safe.) + +### Force kernel (`src/qew_model.h`) + +8. **Each noise edge is sampled twice in `gradient()`** — site `i` samples the + midpoints of edges `(i-1,i)` and `(i,i+1)` (lines 105–106), so every edge + is evaluated by both adjacent sites: `2n` bicubic gathers for `n` distinct + edges. A two-pass scheme — edge kernel writing `dv_dy` (and the arclength + slope term) per edge into a scratch view, then a cheap combine kernel — + does `n` gathers + `~3n` doubles of scratch traffic instead of `2n` + gathers (16 loads each). Expect close to 2× on the dominant kernel. + (`objective()` already samples each edge exactly once; `rk_site_velocity` + re-evaluates at varying `h` and cannot share.) + +9. **`sample()` computes derivatives the caller discards** — + `objective()` uses only `s.v`, yet `DeviceNoise::sample()` + (`src/filtered_noise.h:67`) always accumulates `dv_dgx` *and* `dv_dgy` and + `bspline_weights()` always fills `dw[4]`. `gradient()` never uses `dv_dx` + anywhere in the codebase. Splitting into value-only / value+dy variants + (template tag or separate methods) trims ~⅓–½ of the FLOPs in the 4×4 + inner loop. Gain is modest where memory-bound, but the Armijo line search + is objective-only and benefits directly. + +10. **Integer modulo in the inner loops** — each `sample()` does 8 `%` + operations for the neighbour wrap (`src/filtered_noise.h:79–80`), plus the + stencils' `(i±1) % n`; integer division is expensive on GPUs. After the + positive wrap, `ix ∈ [0, nx]` (the `gx == nx` rounding edge case is + currently absorbed by `ix % nx`), so conditional add/subtract wraps are + sufficient — keep an explicit clamp for that edge case. Power-of-two + grids could use bit masks. + +### Rosso–Krauth (`src/rosso_krauth.h`) + +11. **Redundant mean-`h` reduction per sweep** — `rk_block()` runs two + reductions per sweep (lines 161–167): the advance sum and a fresh mean of + `h`. Since every move is forward, `mean_h = mean_h0 + Σ_sweeps + total_adv / n`; tracking the cumulative advance eliminates the second + reduction. Sweeps run up to `max_sweeps = 2·10⁵`, so this halves the + per-sweep reduction/sync count. + +### Langevin dynamics (`src/dynamics.h`) + +12. **RNG pool acquisition per site per step** — `pool.get_state()` / + `free_state()` inside the `parallel_for` (lines 48–50) acquires a locked + state per site per step. With the benchmark's `nx = 262144` on a GPU the + default-sized `Random_XorShift64_Pool` serializes on lock contention. + Alternatives: size the pool to the launch width, or pre-fill a normal + deviate view per step with `Kokkos::fill_random` and fuse the update. + Also, the `gradient()` call inside `step()` inherits item 8's 2× edge + duplication. + +### Minor / one-time + +13. **FIRE** (`src/fire.h`) — `nv2`/`nf2` are two separate reductions + (lines 103–108) and the v-step (line 79) could fuse with the power + reduction (line 85). Low priority: FIRE is the fallback minimizer. +14. **Noise setup** (`src/filtered_noise.cpp`) — phase generation, mean and + variance are serial host loops over `nx·ny`; `build_coefficients()` stages + through three complex 2-D buffers + host mirrors. One-time cost + (`setup ≈` seconds in the benchmark), only worth touching if many seeds + are ever swept in one run. +15. **`to_host_vec`** (`src/qew_runs.h:28`) — extra element-wise copy from the + mirror into a `std::vector`; could `memcpy` from the mirror when layouts + match. Negligible at current snapshot rates. + +## Suggested order of attack + +Items 5–7 (L-BFGS overhead) and 8 (edge-duplication) are independent, +low-risk, and cover the production `qew_static` runtime; the +gradient-consistency and golden-file tests pin their correctness. Item 11 is +a five-line change covering `qew_depinning`. Item 1 should be fixed (or +asserted away) before anyone runs an odd-`n` line. Items 2 and 4 are +correctness-of-use issues rather than speed, and are cheap insurance. From 3fefe3f4f9a5f3ba81d8defa70fb943a057eeb30 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 17:51:58 +0000 Subject: [PATCH 02/14] doc: phased implementation plan for the optimization-notes findings https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- docs/fix_plan.md | 160 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 docs/fix_plan.md diff --git a/docs/fix_plan.md b/docs/fix_plan.md new file mode 100644 index 0000000..06916c7 --- /dev/null +++ b/docs/fix_plan.md @@ -0,0 +1,160 @@ +# Implementation plan for `optimization_notes.md` + +Companion to [`optimization_notes.md`](optimization_notes.md); item numbers +below refer to that document. Work is grouped into five phases that are +independently land-able, ordered so that correctness fixes go first and each +perf change is benchmarked in isolation. Every phase ends with the full +`ctest` suite on a Serial **and** an OpenMP build (the threaded run doubles as +the race check), plus a before/after run of `executables/benchmark` for the +perf phases. + +Ground rules: + +- The gradient-consistency test and the golden-file tests are the correctness + gates; they must pass unchanged (same tolerances) in every phase. +- Tests that assert `n_func_evals` counts will be updated *deliberately* in + Phase 2, in the same commit as the behaviour change, with the new expected + count derived by hand. +- One commit (or small commit series) per phase, so a regression bisects to a + phase. + +## Phase 0 — Safety & build hygiene (items 1, 3, 4) + +Small, riskless, do first. + +1. **Odd-`n` guard** (`rosso_krauth.h`): in `rk_block()`, before launching + kernels, `throw std::invalid_argument` (host code) when `n % 2 != 0`, + citing the red-black neighbour conflict. A 3-colour sweep for odd `n` is + *not* worth the complexity until someone needs odd lines. Add a death/throw + test. +2. **Integer ramp counter** (`qew_runs.h::depinning_ramp`): iterate + `for (int k = 0; k * f_step <= f_max; ++k)` with `f = k * f_step`. The + default `f_step = 1.0` produces bit-identical forces, so no test or HDF5 + key changes. +3. **CMake**: default `CMAKE_BUILD_TYPE=Release` when unset (with a STATUS + message); rename the project to `qew`. Leave `hello.cpp` / + `kokkos_smoketest.cpp` in place for now — they are referenced by tests and + removing them is cosmetic churn; fold into a later cleanup if desired. + +## Phase 1 — Decision point: line-tension convention (item 2) + +Blocked on information outside the repo: the `../Overleaf/qEW/qew_*.py` +sources. **Needs Lars to confirm** which convention each Python driver uses. +Then, in one commit: + +- Rename the field in whichever config struct holds the ratio to + `pinning_length_over_xi` (or convert at the call site so all three drivers + use one formula). +- Add a tiny unit test per driver pinning the `pinning_length → + line_tension` mapping, so the convention can never silently diverge again. + +Nothing later depends on this; it can land any time. + +## Phase 2 — L-BFGS overhead (items 5–7) + +All inside `lbfgs.h`, in three commits, benchmarking `qew_static`-shaped runs +(preconditioned, n = 4096) after each: + +1. **Carry the objective forward** (item 5). Hoist one + `phi = objective(h, ...)` above the iteration loop; on line-search + acceptance set `phi = phi_accepted`; delete the per-iteration `phi0` + recompute. The line-search-failure `continue` path restores `h = x0`, so + the carried value stays valid there too. Arithmetic is unchanged → + golden L-BFGS test passes bit-for-bit; only `n_func_evals` assertions in + `test_lbfgs.cpp` need recounting. +2. **Fuse trial step into the objective** (item 6). Add a local helper + `trial_objective(x0, d, alpha, p, noise)` — the `objective()` reduce with + `h(i)` replaced by `x0(i) + alpha * d(i)` (and the same for the `ip` + neighbour, which is the identical FP expression, so results stay + bit-identical). Commit `h` with a single `parallel_for` only on + acceptance; the `deep_copy(x0, h)` save and `deep_copy(h, x0)` undo both + disappear. Keep `objective()` itself untouched (public API, used by + tests/benchmark). +3. **Reduction fusion + history swap** (item 7). + - One `parallel_reduce` with three reducers computing `sy`, `ss`, `yy` + while also writing `s_tmp`/`y_tmp` (replaces one `parallel_for` + three + reduces). + - `std::swap(S[row], s_tmp)` / `std::swap(Y[row], y_tmp)` instead of + `deep_copy` (safe: the FFT plans are bound to `q`/`d` only). + - Fold `q -= a*y` into the subsequent dot where profiling on a GPU build + shows the sync count matters; otherwise skip — this is the only + speculative sub-item, and it slightly obfuscates the two-loop recursion. + +## Phase 3 — Hot-kernel work reduction (items 8, 9, 10) + +The order matters: land 10 (mechanical, bit-identical) first so the bigger +rewrites diff cleanly against it. + +1. **Modulo removal** (item 10, `filtered_noise.h::sample`, + `qew_model.h`, `rosso_krauth.h`): replace `(i±1) % n` stencil wraps with + compare-and-adjust; in `sample()`, clamp the `gx == nx` rounding edge case + explicitly (`if (ix >= nx) { ix -= nx; gx -= nx; }` before computing + `fx`), then build the four neighbour indices with conditional wraps. + Same arithmetic → golden noise-interpolation test must pass bit-for-bit. +2. **Value-only / dy-only sampling** (item 9). Add to `DeviceNoise`: + `sample_value(x, y)` (no `dw`, no derivative accumulators) and keep + `sample()` as-is. In `qew_model.h`, route calls through C++20 + `requires`-detected helpers (`noise_value(noise, x, y)` falls back to + `noise.sample(x, y).v` for the analytic test stubs, so no test functor + has to change). Split `bspline_weights` into weights-only and + weights+derivatives variants. `objective()` switches to the value path. +3. **Edge-pass gradient** (item 8). Introduce + `struct GradientWorkspace { View1D edge; }` (sized lazily to `n`). + `gradient(h, p, noise, grad, ws)` becomes two kernels: + - edge kernel: for edge `(i, i+1)` store + `edge(i) = 0.5 * dv_dy(midpoint)` — plus, for the arclength model, the + normalized slope `dh/sqrt(1+dh²)`, packed into a second view or a 2-D + workspace, since both adjacent sites reuse it; + - combine kernel: `grad(i) = dx * (line_force(from edge slopes or h) + + edge(im) + edge(i) - f)`. + Keep a back-compat `gradient(h, p, noise, grad)` overload that allocates a + workspace internally (used by tests); `lbfgs_minimize`, `fire_minimize`, + and `LangevinDynamics` own a persistent workspace. Floating-point sums + re-associate slightly (`0.5*(sl+sr)` becomes `0.5*sl + 0.5*sr`), so the + gradient-consistency test (tolerance-based) is the gate, and the golden + L-BFGS comparison tolerance may need a documented epsilon bump — verify + it stays within the existing tolerance first. + `rk_site_velocity` is intentionally untouched (single-site re-evaluation + cannot share edges). + +Acceptance: benchmark `grad` ms/eval roughly halves on OpenMP; `obj` +improves measurably from item 9; no test tolerance loosened beyond what is +justified in the commit message. + +## Phase 4 — Rosso–Krauth & dynamics (items 11, 12) + +1. **Drop the per-sweep mean reduction** (item 11): track + `cum_adv += total_adv` and compute `res.mean_h = mean_h0 + cum_adv / n`. + Bit-identical apart from summation order of the mean; the RK tests + (washboard `f_c`, runaway) are tolerance-based and unaffected. +2. **Langevin RNG** (item 12): measure first. On the OpenMP build the pool + is per-thread and fine; the concern is GPU lock contention. Run the + benchmark's `dyn` timing on the CUDA CI runner; only if `dyn − grad` + overhead is significant, switch to a pre-generated normal-deviate view + filled by a dedicated kernel where each thread owns one pool state for the + whole row-chunk (avoids per-element lock churn), fused with the update + step. The equipartition test gates correctness; exact trajectories may + change (different draw order) — that test is statistical, so no fixture + updates expected. + +## Phase 5 — Optional cleanups (items 13–15) + +FIRE reduction fusion, host-side noise-setup loops, `to_host_vec` memcpy. +Do only if profiling ever shows them; otherwise close as won't-fix and note +that in `optimization_notes.md`. + +## Sequencing summary + +| Phase | Contents | Risk | Blocked on | +|---|---|---|---| +| 0 | odd-`n` guard, ramp counter, CMake | none | — | +| 1 | line-tension convention | none | Python sources (Lars) | +| 2 | L-BFGS: carry φ, fused trial, fused reduces, view swap | low | — | +| 3 | modulo removal → sample variants → edge-pass gradient | medium | — | +| 4 | RK mean tracking; Langevin RNG (measure-first) | low | GPU runner for 12 | +| 5 | FIRE / setup / copy cleanups | — | profiling evidence | + +Phases 0, 2, 3, 4 are mutually independent except that Phase 3's edge-pass +should rebase on Phase 2's line-search changes (both touch the `lbfgs.h` +call sites). Suggested landing order: 0 → 2 → 3 → 4, with 1 whenever the +Python check happens. From 4bf6f8e021cd6ecf4440d4dddfcbd4273bbb7333 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 18:12:12 +0000 Subject: [PATCH 03/14] fix: phase 0 safety & build hygiene - rk_block: reject odd line lengths up front; with periodic boundaries sites 0 and n-1 share red-black parity but are nearest neighbours, so one colour pass would update both concurrently. Add a test. - depinning_ramp: integer ramp counter instead of accumulating a floating-point loop variable (bit-identical for the default step). - CMake: default CMAKE_BUILD_TYPE to Release; rename the placeholder project name to qew. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- CMakeLists.txt | 9 ++++++++- src/qew_runs.h | 5 ++++- src/rosso_krauth.h | 7 +++++++ tests/test_rosso_krauth.cpp | 23 +++++++++++++++++++++++ 4 files changed, 42 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b652aa6..5fe997d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,14 @@ cmake_minimum_required(VERSION 3.14) if(POLICY CMP0135) cmake_policy(SET CMP0135 NEW) endif() -project(my-fantastic-simulation-code LANGUAGES C CXX VERSION 0.1) +project(qew LANGUAGES C CXX VERSION 0.1) + +# Bandwidth-bound HPC code: an accidental no-flags build (-O0) is never what +# anyone wants, so default to Release when the user does not choose. +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + set(CMAKE_BUILD_TYPE Release CACHE STRING "Build type" FORCE) + message(STATUS "CMAKE_BUILD_TYPE not set, defaulting to Release") +endif() # Kokkos 5.x requires C++20 or newer. set(CMAKE_CXX_STANDARD 20) diff --git a/src/qew_runs.h b/src/qew_runs.h index 486cc25..9c2308a 100644 --- a/src/qew_runs.h +++ b/src/qew_runs.h @@ -211,7 +211,10 @@ DepinningResult depinning_ramp(const DepinningConfig &c, const Noise &noise) { DepinningResult res; res.line_tension = line_tension; - for (real_t f = 0; f <= c.f_max; f += c.f_step) { + // Integer counter: accumulating f += f_step drifts for non-representable + // steps and can drop the final force (and perturb the HDF5 key strings). + for (int k = 0; static_cast(k) * c.f_step <= c.f_max; ++k) { + const real_t f = static_cast(k) * c.f_step; p.driving_force = f; const RkResult r = rk_block(h, p, noise, c.rp); // restart from previous diff --git a/src/rosso_krauth.h b/src/rosso_krauth.h index 800ad24..04719a6 100644 --- a/src/rosso_krauth.h +++ b/src/rosso_krauth.h @@ -1,6 +1,8 @@ #ifndef QEW_ROSSO_KRAUTH_H #define QEW_ROSSO_KRAUTH_H +#include + #include #include "qew_model.h" @@ -127,6 +129,11 @@ template RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, const RkParams &rp) { const int n = static_cast(h.extent(0)); + // On a periodic ring with odd n, sites 0 and n-1 share parity but are + // nearest neighbours, so one colour pass would update both concurrently. + if (n % 2 != 0) + throw std::invalid_argument( + "rk_block: line length must be even for the red-black sweep"); const real_t dx = p.physical_size / static_cast(n); const real_t lt = p.line_tension; const real_t f = p.driving_force; diff --git a/tests/test_rosso_krauth.cpp b/tests/test_rosso_krauth.cpp index 8d2c828..bbf708b 100644 --- a/tests/test_rosso_krauth.cpp +++ b/tests/test_rosso_krauth.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -288,3 +289,25 @@ TEST(RossoKrauth, BlocksArclengthModel) { Kokkos::Max(max_v)); EXPECT_LT(max_v, 1e-6); } + +// An odd line length must be rejected up front: with periodic boundaries sites +// 0 and n-1 share red-black parity but are nearest neighbours, so one colour +// pass would update both concurrently (a data race, not just wrong physics). +TEST(RossoKrauth, RejectsOddLineLength) { + const int n = 17; + const SinusoidNoise noise; + + Params p; + p.physical_size = n; + p.line_tension = 1.0; + p.driving_force = 0.0; + p.model = Model::Linear; + + RkParams rp; + rp.dstep = 0.01; + rp.max_advance = 1.0; + rp.runaway = 3.0; + + View1D h = flat_line(n, 0.0); + EXPECT_THROW(rk_block(h, p, noise, rp), std::invalid_argument); +} From 73a852ba82b948a6add9ecf610864e699a9573e4 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 18:20:04 +0000 Subject: [PATCH 04/14] perf: lbfgs carries the objective value across iterations phi at the current iterate equals the accepted trial value from the previous line search, so recomputing it each iteration was one full objective evaluation (a noise sample per site + reduction) wasted per step -- in the common first-trial-accepted case, half the objective cost. Only the first iteration evaluates it explicitly. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/lbfgs.h | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/lbfgs.h b/src/lbfgs.h index f5c5979..e0607a3 100644 --- a/src/lbfgs.h +++ b/src/lbfgs.h @@ -146,6 +146,8 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, gradient(h, p, noise, g); LbfgsResult res; int ls_fail_streak = 0; + real_t phi = 0; // objective at the current iterate, carried across + // iterations (each accepted trial evaluates it anyway) for (int it = 0; it < lp.max_iter; ++it) { const real_t gmax = lbfgs_max_abs(n, g); @@ -197,8 +199,11 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, // ---- Armijo backtracking line search ---- Kokkos::deep_copy(x0, h); - const real_t phi0 = objective(h, p, noise); // h == x0 here - ++res.n_func_evals; + if (it == 0) { // later iterations reuse the accepted trial value + phi = objective(h, p, noise); + ++res.n_func_evals; + } + const real_t phi0 = phi; // objective at x0 // alpha0 = 1 for L-BFGS / preconditioned steps (well-scaled direction); // for an unpreconditioned first steepest-descent step, scale by 1/|g|. real_t alpha = static_cast(1); @@ -212,9 +217,10 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_for( "lbfgs_trial", n, KOKKOS_LAMBDA(const int i) { h(i) = x0(i) + al * d(i); }); - const real_t phi = objective(h, p, noise); + const real_t phi_trial = objective(h, p, noise); ++res.n_func_evals; - if (phi <= phi0 + lp.c1 * alpha * dphi0) { + if (phi_trial <= phi0 + lp.c1 * alpha * dphi0) { + phi = phi_trial; ok = true; break; } From b3990c58c36d0517e687a8279556f3a4a0656fa4 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 18:22:22 +0000 Subject: [PATCH 05/14] perf: evaluate line-search trials without materialising the trial point objective() is refactored onto a shared body that reads the line position through a device-callable functor, with a new overload evaluating E at x0 + alpha*d directly inside the reduce. The Armijo loop in lbfgs_minimize no longer writes h per trial, so the saved x0 copy and the failed-search undo copy disappear; the accepted step is committed with a single axpy. One analytic definition of the energy is shared by both paths, keeping objective/gradient consistency intact. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/lbfgs.h | 20 +++++++++--------- src/qew_model.h | 55 +++++++++++++++++++++++++++++++++++++++++-------- 2 files changed, 56 insertions(+), 19 deletions(-) diff --git a/src/lbfgs.h b/src/lbfgs.h index e0607a3..84bdafa 100644 --- a/src/lbfgs.h +++ b/src/lbfgs.h @@ -101,7 +101,7 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, const int m = lp.m; View1D g("lbfgs_g", n), g_new("lbfgs_gnew", n), d("lbfgs_d", n); - View1D q("lbfgs_q", n), x0("lbfgs_x0", n), s_tmp("lbfgs_s", n), y_tmp("lbfgs_y", n); + View1D q("lbfgs_q", n), s_tmp("lbfgs_s", n), y_tmp("lbfgs_y", n); std::vector S(m), Y(m); for (int k = 0; k < m; ++k) { S[k] = View1D("lbfgs_S", n); @@ -198,12 +198,13 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, } // ---- Armijo backtracking line search ---- - Kokkos::deep_copy(x0, h); + // Trials are evaluated at h + alpha*d without writing h; the step is + // committed only on acceptance, so failure needs no saved copy/undo. if (it == 0) { // later iterations reuse the accepted trial value phi = objective(h, p, noise); ++res.n_func_evals; } - const real_t phi0 = phi; // objective at x0 + const real_t phi0 = phi; // objective at the current iterate // alpha0 = 1 for L-BFGS / preconditioned steps (well-scaled direction); // for an unpreconditioned first steepest-descent step, scale by 1/|g|. real_t alpha = static_cast(1); @@ -213,11 +214,7 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, } bool ok = false; for (int ls = 0; ls < lp.max_ls; ++ls) { - const real_t al = alpha; - Kokkos::parallel_for( - "lbfgs_trial", n, - KOKKOS_LAMBDA(const int i) { h(i) = x0(i) + al * d(i); }); - const real_t phi_trial = objective(h, p, noise); + const real_t phi_trial = objective(h, d, alpha, p, noise); ++res.n_func_evals; if (phi_trial <= phi0 + lp.c1 * alpha * dphi0) { phi = phi_trial; @@ -227,7 +224,6 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, alpha *= lp.ls_shrink; } if (!ok) { - Kokkos::deep_copy(h, x0); // undo the failed trial if (slots.empty() || ++ls_fail_streak > 2) { res.iterations = it; res.max_force = gmax; @@ -240,9 +236,13 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, } ls_fail_streak = 0; + // Commit the accepted step. + const real_t al = alpha; + Kokkos::parallel_for( + "lbfgs_step", n, KOKKOS_LAMBDA(const int i) { h(i) += al * d(i); }); + // ---- curvature pair s = alpha d, y = g_new - g ---- gradient(h, p, noise, g_new); - const real_t al = alpha; Kokkos::parallel_for( "lbfgs_sy", n, KOKKOS_LAMBDA(const int i) { s_tmp(i) = al * d(i); diff --git a/src/qew_model.h b/src/qew_model.h index c4b0eb1..62c6674 100644 --- a/src/qew_model.h +++ b/src/qew_model.h @@ -26,12 +26,13 @@ struct Params { Model model = Model::Linear; }; -// Total energy E[h] = sum_i dx * (line_energy_i + V_i - f h_i), where the -// elastic and noise terms are evaluated on the forward edge (i, i+1) exactly as -// in qew_objective. Runs as a parallel_reduce on the default execution space. -template -real_t objective(const View1D &h, const Params &p, const Noise &noise) { - const int n = static_cast(h.extent(0)); +namespace detail { + +// Shared objective body: the line position enters only through the +// device-callable functor `at(i)`, so the plain evaluation and the +// line-search trial evaluation share one analytic definition. +template +real_t objective_eval(int n, const Params &p, const Noise &noise, const At &at) { const real_t dx = p.physical_size / static_cast(n); const real_t lt = p.line_tension; const real_t f = p.driving_force; @@ -42,7 +43,9 @@ real_t objective(const View1D &h, const Params &p, const Noise &noise) { "qew::objective", n, KOKKOS_LAMBDA(const int i, real_t &acc) { const int ip = (i + 1) % n; // periodic forward neighbour - const real_t dh_dx = (h(ip) - h(i)) / dx; + const real_t hi = at(i); + const real_t hip = at(ip); + const real_t dh_dx = (hip - hi) / dx; real_t line_energy; if (model == Model::Linear) { @@ -55,15 +58,49 @@ real_t objective(const View1D &h, const Params &p, const Noise &noise) { // Edge midpoint, matching python xcenter/hcenter (periodic wrap). const real_t xc = static_cast(i + ip) * static_cast(0.5) * dx; - const real_t hc = (h(ip) + h(i)) * static_cast(0.5); + const real_t hc = (hip + hi) * static_cast(0.5); const NoiseSample s = noise.sample(xc, hc); - acc += dx * (line_energy + s.v - f * h(i)); + acc += dx * (line_energy + s.v - f * hi); }, total); return total; } +struct LineAt { + View1D h; + KOKKOS_INLINE_FUNCTION real_t operator()(int i) const { return h(i); } +}; + +struct TrialAt { + View1D x0, d; + real_t alpha; + KOKKOS_INLINE_FUNCTION real_t operator()(int i) const { + return x0(i) + alpha * d(i); + } +}; + +} // namespace detail + +// Total energy E[h] = sum_i dx * (line_energy_i + V_i - f h_i), where the +// elastic and noise terms are evaluated on the forward edge (i, i+1) exactly as +// in qew_objective. Runs as a parallel_reduce on the default execution space. +template +real_t objective(const View1D &h, const Params &p, const Noise &noise) { + return detail::objective_eval(static_cast(h.extent(0)), p, noise, + detail::LineAt{h}); +} + +// E at the line-search trial point x0 + alpha*d, evaluated without +// materialising it -- saves a kernel launch and n stores per trial, and the +// caller commits (or discards) the step afterwards. +template +real_t objective(const View1D &x0, const View1D &d, real_t alpha, + const Params &p, const Noise &noise) { + return detail::objective_eval(static_cast(x0.extent(0)), p, noise, + detail::TrialAt{x0, d, alpha}); +} + // Force grad_i = dE/dh_i. Elastic term: harmonic Laplacian (linear) or // slope-normalized difference (arclength); noise term: half the sum of dV/dy at // the left edge (i-1, i) and right edge (i, i+1). Exact port of qew_gradient, From 3bc6706ea712f1c84545f713734f421d1a04983a Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 18:23:39 +0000 Subject: [PATCH 06/14] perf: fuse lbfgs curvature-pair kernels, swap views instead of copying - s/y pair construction and the sy/ss/yy dot products collapse from one parallel_for plus three parallel_reduces into a single fused reduce (one pass over the data, one device sync instead of four). - The accepted pair is std::swap'ed into the history slot and g/g_new are swapped, replacing three full-vector deep_copies per iteration (safe: the FFT plans are bound to q and d only). https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/lbfgs.h | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/lbfgs.h b/src/lbfgs.h index 84bdafa..4a5bb3d 100644 --- a/src/lbfgs.h +++ b/src/lbfgs.h @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -242,15 +243,21 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, "lbfgs_step", n, KOKKOS_LAMBDA(const int i) { h(i) += al * d(i); }); // ---- curvature pair s = alpha d, y = g_new - g ---- + // One fused kernel writes the pair and reduces all three scalars. gradient(h, p, noise, g_new); - Kokkos::parallel_for( - "lbfgs_sy", n, KOKKOS_LAMBDA(const int i) { - s_tmp(i) = al * d(i); - y_tmp(i) = g_new(i) - g(i); - }); - const real_t sy = lbfgs_dot(n, s_tmp, y_tmp); - const real_t ss = lbfgs_dot(n, s_tmp, s_tmp); - const real_t yy = lbfgs_dot(n, y_tmp, y_tmp); + real_t sy = 0, ss = 0, yy = 0; + Kokkos::parallel_reduce( + "lbfgs_sy", n, + KOKKOS_LAMBDA(const int i, real_t &asy, real_t &ass, real_t &ayy) { + const real_t s = al * d(i); + const real_t y = g_new(i) - g(i); + s_tmp(i) = s; + y_tmp(i) = y; + asy += s * y; + ass += s * s; + ayy += y * y; + }, + sy, ss, yy); if (sy > lp.curvature_eps * std::sqrt(ss * yy)) { // skip safeguard int row; if (static_cast(slots.size()) < m) { @@ -259,15 +266,17 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, row = slots.front(); slots.erase(slots.begin()); } - Kokkos::deep_copy(S[row], s_tmp); - Kokkos::deep_copy(Y[row], y_tmp); + // Swap the pair into the history slot (s_tmp/y_tmp become the + // slot's old buffers, reused as scratch next iteration). + std::swap(S[row], s_tmp); + std::swap(Y[row], y_tmp); rho[row] = static_cast(1) / sy; gamma = sy / yy; slots.push_back(row); } else { ++res.n_curvature_skips; // non-positive curvature: keep H0 pos. def. } - Kokkos::deep_copy(g, g_new); + std::swap(g, g_new); // g <- g_new; old g becomes scratch } res.iterations = lp.max_iter; From 45a9139f2fbfaad40edaca09e15a055fa68fcaea Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 18:26:58 +0000 Subject: [PATCH 07/14] bench: time an end-to-end preconditioned L-BFGS relaxation The per-kernel timings cannot see minimiser-loop overheads (line-search objective evals, two-loop recursion, history updates), which is where the lbfgs optimisations land. Measured on this machine (Serial, n=4096, arclength, FilteredNoise): 5.73s/12893 fevals before the lbfgs changes, 4.77s/10424 fevals after, at an identical 2470 iterations. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- executables/benchmark/main.cpp | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/executables/benchmark/main.cpp b/executables/benchmark/main.cpp index 516e34a..77403e0 100644 --- a/executables/benchmark/main.cpp +++ b/executables/benchmark/main.cpp @@ -10,6 +10,7 @@ #include "dynamics.h" #include "filtered_noise.h" +#include "lbfgs.h" #include "qew_model.h" #include "qew_types.h" @@ -69,6 +70,35 @@ int main(int argc, char *argv[]) { "dyn=%7.3f ms/ea\n", nthreads, t_setup, 1e3 * t_grad / reps, 1e3 * t_obj / reps, 1e3 * t_dyn / reps); + + // End-to-end L-BFGS relaxation shaped like a qew_static run + // (preconditioned, flat start). Tracks the minimiser-loop overheads + // (line-search objective evals, two-loop recursion, history updates) + // that the per-kernel timings above cannot see. + const int snx = 4096, sny = 512; + const real_t sLx = 32.0, sLy = 1.0, sxi = 0.1; + FilteredNoise snoise(snx, sny, sLx, sLy, 1.0, sxi, sxi, /*seed=*/1); + Params sp; + sp.physical_size = sLx; + sp.line_tension = 1.0; // pinning_length = 1 + sp.driving_force = 0.0; + sp.model = Model::Arclength; + View1D hs("hs", snx); + Kokkos::deep_copy(hs, static_cast(0.5) * sLy); + LbfgsParams sopt; + sopt.ftol = 1e-6; + sopt.max_iter = 5000; + sopt.precondition = true; + sopt.precond_shift = + static_cast(0.5) * (sLx / snx) / (sxi * sxi); + timer.reset(); + const LbfgsResult sres = + lbfgs_minimize(hs, sp, snoise.device_noise(), sopt); + const double t_lbfgs = timer.seconds(); + std::printf( + "lbfgs(n=%d, precond): %6.3fs iters=%d feval=%ld %s\n", snx, + t_lbfgs, sres.iterations, sres.n_func_evals, + sres.converged ? "converged" : "NOT converged"); } Kokkos::finalize(); return 0; From 87becda3ed5a2b5b7b90edb4b46c50eb051bcbfe Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 18:27:34 +0000 Subject: [PATCH 08/14] doc: mark phases 0 and 2 done, record measured lbfgs speedup https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- docs/fix_plan.md | 7 +++++++ docs/optimization_notes.md | 7 ++++++- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/docs/fix_plan.md b/docs/fix_plan.md index 06916c7..42b398a 100644 --- a/docs/fix_plan.md +++ b/docs/fix_plan.md @@ -1,5 +1,12 @@ # Implementation plan for `optimization_notes.md` +> **Status (June 2026):** Phases 0 and 2 are implemented and merged into this +> branch. Measured effect of Phase 2 on the benchmark's end-to-end +> preconditioned L-BFGS relaxation (Serial, n = 4096, arclength): 5.73 s / +> 12893 objective evals → 4.77 s / 10424 at an identical 2470 iterations +> (−17% wall time). Phases 1, 3, 4, 5 are open; Phase 1 still needs the +> Python-source check described below. + Companion to [`optimization_notes.md`](optimization_notes.md); item numbers below refer to that document. Work is grouped into five phases that are independently land-able, ordered so that correctness fixes go first and each diff --git a/docs/optimization_notes.md b/docs/optimization_notes.md index d96fd46..77f49e8 100644 --- a/docs/optimization_notes.md +++ b/docs/optimization_notes.md @@ -1,7 +1,12 @@ # Optimization potential & known issues Technical notes from a code review of the qEW port (June 2026). Items are -ordered by expected impact. File references point at the current `main`-line +ordered by expected impact. + +> **Status:** items 1, 3, 4 (bugs) and 5, 6, 7 (L-BFGS overhead) are fixed on +> this branch — see `fix_plan.md` for the phase breakdown and measured +> numbers. Items 2 (convention check), 8–12 (kernel work reduction, RK, +> dynamics) and the minor items remain open. File references point at the current `main`-line sources. Nothing here changes physics; everything preserves the documented analytic consistency between `objective()` and `gradient()` and the golden-file contracts. From f83a4621bfc7c06f709d9ceebfe28046fbc41cba Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 19:01:12 +0000 Subject: [PATCH 09/14] fix: unify the pinning-length convention on the Larkin ratio lambda_p/xi The paper (arXiv:2410.21838, eq. (10)) controls the model through the dimensionless ratio lambda_p/xi = (Gamma/U0)^{2/3}, varied by scanning the line tension at fixed disorder. Previously the same config field name meant that ratio in static_sweep/depinning_ramp but a physical length (divided by xi internally) in dynamic_anneal. All drivers now take pinning_length_over_xi and share Gamma = amplitude * (lambda_p/xi)^{3/2}, pinned by a mapping assertion per driver. Line tensions are numerically unchanged in every run. HDF5 key strings are untouched; the values carry the ratio, which already held for static/depinning output (the dynamic demo's recorded value changes 0.01 -> 0.1 accordingly). https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- docs/fix_plan.md | 16 ++++++++----- docs/optimization_notes.md | 27 +++++++++++----------- executables/benchmark/main.cpp | 2 +- executables/qew_depinning/main.cpp | 6 ++--- executables/qew_dynamic/main.cpp | 4 +++- executables/qew_static/main.cpp | 8 ++++--- src/qew_io.h | 13 ++++++----- src/qew_runs.h | 36 +++++++++++++++++++----------- tests/test_executables.cpp | 29 +++++++++++++++--------- 9 files changed, 85 insertions(+), 56 deletions(-) diff --git a/docs/fix_plan.md b/docs/fix_plan.md index 42b398a..38ef92f 100644 --- a/docs/fix_plan.md +++ b/docs/fix_plan.md @@ -1,11 +1,15 @@ # Implementation plan for `optimization_notes.md` -> **Status (June 2026):** Phases 0 and 2 are implemented and merged into this -> branch. Measured effect of Phase 2 on the benchmark's end-to-end -> preconditioned L-BFGS relaxation (Serial, n = 4096, arclength): 5.73 s / -> 12893 objective evals → 4.77 s / 10424 at an identical 2470 iterations -> (−17% wall time). Phases 1, 3, 4, 5 are open; Phase 1 still needs the -> Python-source check described below. +> **Status (June 2026):** Phases 0, 1 and 2 are implemented on this branch. +> Phase 1 was settled by the paper (arXiv:2410.21838 eq. (10): +> `lambda_p/xi = (Gamma/U0)^{2/3}`, varied by scanning Gamma at fixed +> disorder): all drivers now take the dimensionless ratio as +> `pinning_length_over_xi` and share `Gamma = A (lambda_p/xi)^{3/2}`, with +> per-driver mapping tests; line tensions are numerically unchanged. Measured +> effect of Phase 2 on the benchmark's end-to-end preconditioned L-BFGS +> relaxation (Serial, n = 4096, arclength): 5.73 s / 12893 objective evals → +> 4.77 s / 10424 at an identical 2470 iterations (−17% wall time). Phases 3, +> 4, 5 are open. Companion to [`optimization_notes.md`](optimization_notes.md); item numbers below refer to that document. Work is grouped into five phases that are diff --git a/docs/optimization_notes.md b/docs/optimization_notes.md index 77f49e8..32710d7 100644 --- a/docs/optimization_notes.md +++ b/docs/optimization_notes.md @@ -3,10 +3,10 @@ Technical notes from a code review of the qEW port (June 2026). Items are ordered by expected impact. -> **Status:** items 1, 3, 4 (bugs) and 5, 6, 7 (L-BFGS overhead) are fixed on -> this branch — see `fix_plan.md` for the phase breakdown and measured -> numbers. Items 2 (convention check), 8–12 (kernel work reduction, RK, -> dynamics) and the minor items remain open. File references point at the current `main`-line +> **Status:** items 1, 2, 3, 4 (bugs) and 5, 6, 7 (L-BFGS overhead) are fixed +> on this branch — see `fix_plan.md` for the phase breakdown and measured +> numbers. Items 8–12 (kernel work reduction, RK, dynamics) and the minor +> items remain open. File references point at the current `main`-line sources. Nothing here changes physics; everything preserves the documented analytic consistency between `objective()` and `gradient()` and the golden-file contracts. @@ -28,15 +28,16 @@ kernel-launch / device-sync overhead in the minimizers. 3-colour sweep for odd `n`. 2. **Inconsistent pinning-length → line-tension map across run drivers** — - `src/qew_runs.h`: `static_sweep()` (line 67) and `depinning_ramp()` - (line 201) use `line_tension = pinning_length^1.5 * amplitude`, while - `dynamic_anneal()` (line 134) uses `(pinning_length / xi)^1.5 * amplitude`. - The depinning `main` comments that its `pinning_length` is already - `lambda_p/xi`, so each driver may individually match its Python source, but - the same struct-field name means two different physical quantities in - different drivers. *Action:* verify each against `qew_*.py` and either - unify the convention or rename the fields (`pinning_length_over_xi`). No - test currently pins this mapping. + `src/qew_runs.h`: `static_sweep()` and `depinning_ramp()` used + `line_tension = pinning_length^1.5 * amplitude`, while `dynamic_anneal()` + used `(pinning_length / xi)^1.5 * amplitude` — the same struct-field name + meant two different physical quantities. **Resolved:** the paper + (arXiv:2410.21838, eq. (10)) controls the model through the dimensionless + Larkin ratio `lambda_p/xi = (Gamma/U0)^{2/3}`, so all drivers now take a + `pinning_length_over_xi` field and share `Gamma = A * (lambda_p/xi)^{3/2}`; + each driver's mapping is pinned by a unit test. HDF5 key strings are + unchanged and their values carry the ratio (already true for static and + depinning output; dynamic output previously recorded the physical length). 3. **Floating-point loop counter in the force ramp** — `src/qew_runs.h:214`: `for (real_t f = 0; f <= c.f_max; f += c.f_step)`. diff --git a/executables/benchmark/main.cpp b/executables/benchmark/main.cpp index 77403e0..0a36599 100644 --- a/executables/benchmark/main.cpp +++ b/executables/benchmark/main.cpp @@ -80,7 +80,7 @@ int main(int argc, char *argv[]) { FilteredNoise snoise(snx, sny, sLx, sLy, 1.0, sxi, sxi, /*seed=*/1); Params sp; sp.physical_size = sLx; - sp.line_tension = 1.0; // pinning_length = 1 + sp.line_tension = 1.0; // lambda_p/xi = 1 sp.driving_force = 0.0; sp.model = Model::Arclength; View1D hs("hs", snx); diff --git a/executables/qew_depinning/main.cpp b/executables/qew_depinning/main.cpp index b171b5d..2242056 100644 --- a/executables/qew_depinning/main.cpp +++ b/executables/qew_depinning/main.cpp @@ -31,8 +31,8 @@ int main(int argc, char *argv[]) { cfg.model = Model::Arclength; const char *model_name = "arclength"; cfg.amplitude = 1.0; - cfg.xi = 0.1; // correlation length - cfg.pinning_length = 5.0; // lambda_p/xi (stiff) + cfg.xi = 0.1; // correlation length + cfg.pinning_length_over_xi = 5.0; // lambda_p/xi (stiff) cfg.Lx = 8.0; cfg.Ly = 1.0; cfg.nx = 1024; @@ -52,7 +52,7 @@ int main(int argc, char *argv[]) { const DepinningResult res = depinning_ramp(cfg, noise.device_noise()); std::printf("model=%s lambda_p/xi=%.3g line_tension=%.4g\n", model_name, - static_cast(cfg.pinning_length), + static_cast(cfg.pinning_length_over_xi), static_cast(res.line_tension)); for (const DepinningStep &st : res.steps) { std::printf("f=%7.3f sweeps=%-7d mean_h=%8.4f %s\n", diff --git a/executables/qew_dynamic/main.cpp b/executables/qew_dynamic/main.cpp index 10311e2..4350c6d 100644 --- a/executables/qew_dynamic/main.cpp +++ b/executables/qew_dynamic/main.cpp @@ -30,7 +30,9 @@ int main(int argc, char *argv[]) { const char *model_name = "arclength"; cfg.amplitude = 1.0; cfg.xi = 0.1; // correlation length - cfg.pinning_length = 0.01; + // lambda_p/xi = 0.1 (floppy); previously stored as the physical length + // 0.01 and divided by xi -- the line tension is unchanged. + cfg.pinning_length_over_xi = 0.1; cfg.driving_force = 0.0; cfg.Lx = 2.0; cfg.Ly = 1.0; diff --git a/executables/qew_static/main.cpp b/executables/qew_static/main.cpp index 963b706..159a6de 100644 --- a/executables/qew_static/main.cpp +++ b/executables/qew_static/main.cpp @@ -36,7 +36,8 @@ int main(int argc, char *argv[]) { cfg.nx = 4096; cfg.ny = 512; cfg.driving_force = 0.0; - cfg.pinning_lengths = {0.01, 0.1, 1.0, 10.0, 100.0}; + // lambda_p/xi: floppy (0.01) through stiff (100), cf. the paper's sweep. + cfg.pinning_lengths_over_xi = {0.01, 0.1, 1.0, 10.0, 100.0}; cfg.ftol = 1e-6; // ample for roughness analysis cfg.max_iter = 5000; const std::uint64_t seed = 1; @@ -47,8 +48,9 @@ int main(int argc, char *argv[]) { static_sweep(cfg, noise.device_noise()); for (const StaticProfile &prof : profiles) { - std::printf("lp=%-6.3g lt=%-10.4g iters=%-7d max|F|=%.3e %s\n", - prof.pinning_length, static_cast(prof.line_tension), + std::printf("lp/xi=%-6.3g lt=%-10.4g iters=%-7d max|F|=%.3e %s\n", + prof.pinning_length_over_xi, + static_cast(prof.line_tension), prof.result.iterations, static_cast(prof.result.max_force), prof.result.converged ? "converged" : "NOT converged"); diff --git a/src/qew_io.h b/src/qew_io.h index 4da78d7..bdb7f7b 100644 --- a/src/qew_io.h +++ b/src/qew_io.h @@ -19,13 +19,16 @@ namespace qew { +// NB: the "pinning_length" keys/values carry the dimensionless ratio +// lambda_p/xi, the paper's control parameter (see qew_runs.h). inline void write_static(HighFive::File &file, const StaticConfig &c, const std::vector &profiles) { file.createDataSet("physical_size", static_cast(c.Lx)); for (const StaticProfile &prof : profiles) { - const std::string prefix = "driving_force=" + fmt("%.1g", c.driving_force) + - "/pinning_length=" + fmt("%.3g", prof.pinning_length); - file.createDataSet(prefix + "/pinning_length", prof.pinning_length); + const std::string prefix = + "driving_force=" + fmt("%.1g", c.driving_force) + + "/pinning_length=" + fmt("%.3g", prof.pinning_length_over_xi); + file.createDataSet(prefix + "/pinning_length", prof.pinning_length_over_xi); file.createDataSet(prefix + "/h", prof.h); } } @@ -39,7 +42,7 @@ inline void write_dynamic(HighFive::File &file, const DynamicConfig &c, "/driving_force=" + fmt("%.1g", c.driving_force) + "/line_tension=" + fmt("%.3g", res.line_tension); file.createDataSet(prefix + "/pinning_length", - static_cast(c.pinning_length)); + static_cast(c.pinning_length_over_xi)); file.createDataSet(prefix + "/h", snap.h); } } @@ -52,7 +55,7 @@ inline void write_depinning(HighFive::File &file, const DepinningConfig &c, const std::string prefix = "driving_force=" + fmt("%.4g", st.f); file.createDataSet(prefix + "/h", st.h); file.createDataSet(prefix + "/pinning_length", - static_cast(c.pinning_length)); + static_cast(c.pinning_length_over_xi)); } if (res.f_c_upper_bound >= 0) file.createDataSet("f_c_upper_bound", diff --git a/src/qew_runs.h b/src/qew_runs.h index 9c2308a..1bc5d18 100644 --- a/src/qew_runs.h +++ b/src/qew_runs.h @@ -22,7 +22,16 @@ // the kernels it drives, so tests can pass a cheap analytic potential instead of // paying for a FilteredNoise FFT setup. The loop bodies are a faithful copy of // the original mains -- behaviour (timestep heuristic, anneal schedule, force -// ramp + f_c bracketing, pinning-length -> line-tension map) is unchanged. +// ramp + f_c bracketing) is unchanged. +// +// Pinning-length convention (Le et al., arXiv:2410.21838, eq. (10)): the Larkin +// length obeys lambda_p / xi = (Gamma / U0)^{2/3}, varied by holding the +// disorder (amplitude U0, correlation length xi) fixed and scanning the line +// tension Gamma. All drivers therefore take the DIMENSIONLESS ratio +// lambda_p/xi and invert it as +// Gamma = amplitude * (lambda_p/xi)^{3/2}. +// The ratio is also what the HDF5 "pinning_length" keys/values carry, matching +// the paper's labelling (lambda_p/xi = 0.1 floppy ... 100 stiff). namespace qew { inline std::vector to_host_vec(const View1D &h) { @@ -39,19 +48,20 @@ inline std::vector to_host_vec(const View1D &h) { struct StaticConfig { Model model = Model::Arclength; real_t amplitude = 1.0; - real_t xi = 0.1; // correlation length + real_t xi = 0.1; // correlation length; MUST match the noise's xi real_t Lx = 32.0; real_t Ly = 1.0; int nx = 4096; int ny = 512; // noise y-resolution (used only at noise setup) real_t driving_force = 0.0; - std::vector pinning_lengths = {0.01, 0.1, 1.0, 10.0, 100.0}; + // lambda_p/xi ratios to sweep, floppy (<1) to stiff (>1). + std::vector pinning_lengths_over_xi = {0.01, 0.1, 1.0, 10.0, 100.0}; real_t ftol = 1e-6; int max_iter = 5000; }; struct StaticProfile { - double pinning_length = 0; + double pinning_length_over_xi = 0; real_t line_tension = 0; std::vector h; LbfgsResult result; @@ -60,8 +70,8 @@ struct StaticProfile { template std::vector static_sweep(const StaticConfig &c, const Noise &noise) { std::vector out; - out.reserve(c.pinning_lengths.size()); - for (double lp : c.pinning_lengths) { + out.reserve(c.pinning_lengths_over_xi.size()); + for (double lp : c.pinning_lengths_over_xi) { Params p; p.physical_size = c.Lx; p.line_tension = std::pow(lp, static_cast(1.5)) * c.amplitude; @@ -82,7 +92,7 @@ std::vector static_sweep(const StaticConfig &c, const Noise &nois const LbfgsResult r = lbfgs_minimize(h, p, noise, opt); StaticProfile prof; - prof.pinning_length = lp; + prof.pinning_length_over_xi = lp; prof.line_tension = p.line_tension; prof.h = to_host_vec(h); prof.result = r; @@ -97,8 +107,8 @@ std::vector static_sweep(const StaticConfig &c, const Noise &nois struct DynamicConfig { Model model = Model::Arclength; real_t amplitude = 1.0; - real_t xi = 0.1; - real_t pinning_length = 0.01; + real_t xi = 0.1; // correlation length; MUST match the noise's xi + real_t pinning_length_over_xi = 0.1; // lambda_p/xi (floppy) real_t driving_force = 0.0; real_t Lx = 2.0; real_t Ly = 1.0; @@ -131,7 +141,7 @@ DynamicResult dynamic_anneal(const DynamicConfig &c, const Noise &noise) { const real_t dt = std::min(static_cast(0.1), tau / 10 * std::min(fac, fac * fac)); const real_t line_tension = - std::pow(c.pinning_length / c.xi, static_cast(1.5)) * c.amplitude; + std::pow(c.pinning_length_over_xi, static_cast(1.5)) * c.amplitude; Params p; p.physical_size = c.Lx; @@ -172,8 +182,8 @@ DynamicResult dynamic_anneal(const DynamicConfig &c, const Noise &noise) { struct DepinningConfig { Model model = Model::Arclength; real_t amplitude = 1.0; - real_t xi = 0.1; - real_t pinning_length = 5.0; + real_t xi = 0.1; // correlation length; MUST match the noise's xi + real_t pinning_length_over_xi = 5.0; // lambda_p/xi (stiff) real_t Lx = 8.0; real_t Ly = 1.0; int nx = 1024; @@ -198,7 +208,7 @@ struct DepinningResult { template DepinningResult depinning_ramp(const DepinningConfig &c, const Noise &noise) { const real_t line_tension = - std::pow(c.pinning_length, static_cast(1.5)) * c.amplitude; + std::pow(c.pinning_length_over_xi, static_cast(1.5)) * c.amplitude; Params p; p.physical_size = c.Lx; diff --git a/tests/test_executables.cpp b/tests/test_executables.cpp index 69af939..556b267 100644 --- a/tests/test_executables.cpp +++ b/tests/test_executables.cpp @@ -65,19 +65,21 @@ TEST(Executables, StaticSweepShapeAndMapping) { cfg.Ly = 1.0; cfg.nx = 16; cfg.ny = 8; - cfg.pinning_lengths = {0.1, 1.0, 4.0}; + cfg.pinning_lengths_over_xi = {0.1, 1.0, 4.0}; cfg.ftol = 1e-8; cfg.max_iter = 2000; const auto profiles = static_sweep(cfg, ZeroNoise{}); - ASSERT_EQ(profiles.size(), cfg.pinning_lengths.size()); + ASSERT_EQ(profiles.size(), cfg.pinning_lengths_over_xi.size()); for (std::size_t i = 0; i < profiles.size(); ++i) { const auto &prof = profiles[i]; - EXPECT_EQ(prof.pinning_length, cfg.pinning_lengths[i]); - // line_tension = pinning_length^1.5 * amplitude (the documented map). + EXPECT_EQ(prof.pinning_length_over_xi, cfg.pinning_lengths_over_xi[i]); + // Larkin inversion Gamma = amplitude * (lambda_p/xi)^{3/2} + // (arXiv:2410.21838 eq. (10)); shared by all three drivers. EXPECT_NEAR(prof.line_tension, - std::pow(cfg.pinning_lengths[i], 1.5) * cfg.amplitude, 1e-12); + std::pow(cfg.pinning_lengths_over_xi[i], 1.5) * cfg.amplitude, + 1e-12); ASSERT_EQ(prof.h.size(), static_cast(cfg.nx)); // Zero noise, zero force: the flat start is already the minimum. EXPECT_TRUE(prof.result.converged); @@ -94,7 +96,7 @@ TEST(Executables, DynamicAnnealTimestepAndSnapshots) { cfg.model = Model::Linear; cfg.amplitude = 1.0; cfg.xi = 0.1; - cfg.pinning_length = 0.01; + cfg.pinning_length_over_xi = 0.1; cfg.Lx = 1.0; cfg.Ly = 1.0; cfg.nx = 32; @@ -111,8 +113,9 @@ TEST(Executables, DynamicAnnealTimestepAndSnapshots) { const real_t dt_expected = std::min(static_cast(0.1), cfg.tau / 10 * std::min(fac, fac * fac)); EXPECT_NEAR(res.dt, dt_expected, 1e-12); + // Same Larkin map as the static driver: Gamma = A * (lambda_p/xi)^{3/2}. EXPECT_NEAR(res.line_tension, - std::pow(cfg.pinning_length / cfg.xi, 1.5) * cfg.amplitude, 1e-12); + std::pow(cfg.pinning_length_over_xi, 1.5) * cfg.amplitude, 1e-12); // Snapshots are logarithmically spaced: first near 10*dt, each ~10x the last. ASSERT_GE(res.snapshots.size(), 2u); @@ -137,7 +140,7 @@ TEST(Executables, DepinningRampBracketsCriticalForce) { DepinningConfig cfg; cfg.model = Model::Linear; cfg.amplitude = 1.0; - cfg.pinning_length = 1.0; + cfg.pinning_length_over_xi = 1.0; cfg.Lx = 8.0; cfg.Ly = 1.0; cfg.nx = 8; @@ -153,6 +156,10 @@ TEST(Executables, DepinningRampBracketsCriticalForce) { const DepinningResult res = depinning_ramp(cfg, noise); + // Same Larkin map as the other drivers: Gamma = A * (lambda_p/xi)^{3/2}. + EXPECT_NEAR(res.line_tension, + std::pow(cfg.pinning_length_over_xi, 1.5) * cfg.amplitude, 1e-12); + ASSERT_FALSE(res.steps.empty()); EXPECT_TRUE(res.steps.front().result.blocked); // f = 0 blocks EXPECT_EQ(res.steps.front().f, 0.0); @@ -195,7 +202,7 @@ TEST(Executables, WriteStaticKeysAndValues) { cfg.driving_force = 0.0; std::vector profiles(1); - profiles[0].pinning_length = 0.1; + profiles[0].pinning_length_over_xi = 0.1; profiles[0].line_tension = 0.5; profiles[0].h = {1.0, 2.0, 3.0, 4.0}; @@ -225,7 +232,7 @@ TEST(Executables, WriteDynamicKeysAndValues) { DynamicConfig cfg; cfg.Lx = 2.0; cfg.driving_force = 0.0; - cfg.pinning_length = 0.01; + cfg.pinning_length_over_xi = 0.1; DynamicResult res; res.line_tension = 2.0; @@ -249,7 +256,7 @@ TEST(Executables, WriteDynamicKeysAndValues) { TEST(Executables, WriteDepinningSkipsRunawayAndWritesFc) { DepinningConfig cfg; cfg.Lx = 8.0; - cfg.pinning_length = 5.0; + cfg.pinning_length_over_xi = 5.0; DepinningResult res; res.f_c_upper_bound = 0.8; From 74c82c3d2dc5001edbedf09d2cd1c3b36c25a25c Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 19:08:11 +0000 Subject: [PATCH 10/14] perf: conditional periodic wraps + single-quantity noise samplers - Replace integer % with compare-and-adjust wraps in the bicubic sample (8 mods per gather), the model stencils, and the Rosso-Krauth sweep; the gx == nx rounding edge case is folded back explicitly. Integer division is expensive on GPUs and these are the hot loops. - DeviceNoise grows sample_value() (objective path: no dw weights, no derivative accumulators) and sample_dy() (gradient/RK path: dV/dx was computed everywhere and consumed nowhere). Kernels route through requires-detected helpers, so analytic test noises that only provide sample() keep working unchanged. Arithmetic is unchanged; the golden noise-interpolation and golden L-BFGS cross-checks pass. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/filtered_noise.h | 91 ++++++++++++++++++++++++++++++++++++++------ src/qew_model.h | 35 ++++++++++++----- src/rosso_krauth.h | 12 +++--- 3 files changed, 112 insertions(+), 26 deletions(-) diff --git a/src/filtered_noise.h b/src/filtered_noise.h index 5106b1c..731cb02 100644 --- a/src/filtered_noise.h +++ b/src/filtered_noise.h @@ -29,10 +29,9 @@ namespace qew { // Cubic B-spline interpolation weights w[0..3] for the four neighbours -// {i-1, i, i+1, i+2} at fractional position t in [0,1), and their derivatives -// dw[0..3] with respect to t (i.e. d/d(grid coordinate)). +// {i-1, i, i+1, i+2} at fractional position t in [0,1). KOKKOS_INLINE_FUNCTION -void bspline_weights(real_t t, real_t w[4], real_t dw[4]) { +void bspline_weights_value(real_t t, real_t w[4]) { const real_t t2 = t * t; const real_t t3 = t2 * t; const real_t omt = static_cast(1) - t; @@ -44,7 +43,16 @@ void bspline_weights(real_t t, real_t w[4], real_t dw[4]) { static_cast(3) * t2 - static_cast(3) * t3) / static_cast(6); w[3] = t3 / static_cast(6); +} + +// As above, plus the weight derivatives dw[0..3] with respect to t (i.e. +// d/d(grid coordinate)). +KOKKOS_INLINE_FUNCTION +void bspline_weights(real_t t, real_t w[4], real_t dw[4]) { + bspline_weights_value(t, w); + const real_t t2 = t * t; + const real_t omt = static_cast(1) - t; dw[0] = -omt * omt / static_cast(2); dw[1] = (-static_cast(12) * t + static_cast(9) * t2) / static_cast(6); @@ -63,21 +71,40 @@ struct DeviceNoise { real_t dx = 1; // physical grid spacing in x real_t dy = 1; // physical grid spacing in y + // Physical -> grid coordinates: wrap into [0, n), split into the base cell + // index and fraction, and build the four periodic neighbour indices with + // conditional wraps (integer % is expensive in these hot device loops). KOKKOS_INLINE_FUNCTION - NoiseSample sample(real_t x, real_t y) const { - // Physical -> grid coordinates, wrapped into [0, n) (positive modulo). + void locate(real_t x, real_t y, int ixn[4], int iyn[4], real_t &fx, + real_t &fy) const { real_t gx = x / dx; real_t gy = y / dy; gx -= Kokkos::floor(gx / static_cast(nx)) * static_cast(nx); gy -= Kokkos::floor(gy / static_cast(ny)) * static_cast(ny); - const int ix = static_cast(Kokkos::floor(gx)); - const int iy = static_cast(Kokkos::floor(gy)); - const real_t fx = gx - static_cast(ix); - const real_t fy = gy - static_cast(iy); + int ix = static_cast(Kokkos::floor(gx)); + int iy = static_cast(Kokkos::floor(gy)); + // The wrap can round to exactly n; fold it back before the fraction. + if (ix >= nx) { ix -= nx; gx -= static_cast(nx); } + if (iy >= ny) { iy -= ny; gy -= static_cast(ny); } + fx = gx - static_cast(ix); + fy = gy - static_cast(iy); + + ixn[0] = (ix == 0) ? nx - 1 : ix - 1; + ixn[1] = ix; + ixn[2] = (ix + 1 == nx) ? 0 : ix + 1; + ixn[3] = (ixn[2] + 1 == nx) ? 0 : ixn[2] + 1; + iyn[0] = (iy == 0) ? ny - 1 : iy - 1; + iyn[1] = iy; + iyn[2] = (iy + 1 == ny) ? 0 : iy + 1; + iyn[3] = (iyn[2] + 1 == ny) ? 0 : iyn[2] + 1; + } - const int ixn[4] = {(ix - 1 + nx) % nx, ix % nx, (ix + 1) % nx, (ix + 2) % nx}; - const int iyn[4] = {(iy - 1 + ny) % ny, iy % ny, (iy + 1) % ny, (iy + 2) % ny}; + KOKKOS_INLINE_FUNCTION + NoiseSample sample(real_t x, real_t y) const { + int ixn[4], iyn[4]; + real_t fx, fy; + locate(x, y, ixn, iyn, fx, fy); real_t wx[4], dwx[4], wy[4], dwy[4]; bspline_weights(fx, wx, dwx); @@ -105,6 +132,48 @@ struct DeviceNoise { s.dv_dy = dv_dgy / dy; return s; } + + // Value only -- the objective() hot path (line-search trials) does not + // need the derivatives, so skip the dw weights and both accumulators. + KOKKOS_INLINE_FUNCTION + real_t sample_value(real_t x, real_t y) const { + int ixn[4], iyn[4]; + real_t fx, fy; + locate(x, y, ixn, iyn, fx, fy); + + real_t wx[4], wy[4]; + bspline_weights_value(fx, wx); + bspline_weights_value(fy, wy); + + real_t v = 0; + for (int a = 0; a < 4; ++a) { + real_t row = 0; + for (int b = 0; b < 4; ++b) row += coef(ixn[a], iyn[b]) * wy[b]; + v += wx[a] * row; + } + return v; + } + + // dV/dy only -- gradient() and the Rosso-Krauth site velocity use just the + // y-derivative at edge midpoints; dV/dx is never consumed anywhere. + KOKKOS_INLINE_FUNCTION + real_t sample_dy(real_t x, real_t y) const { + int ixn[4], iyn[4]; + real_t fx, fy; + locate(x, y, ixn, iyn, fx, fy); + + real_t wx[4], wy[4], dwy[4]; + bspline_weights_value(fx, wx); + bspline_weights(fy, wy, dwy); + + real_t dv_dgy = 0; + for (int a = 0; a < 4; ++a) { + real_t row_dy = 0; + for (int b = 0; b < 4; ++b) row_dy += coef(ixn[a], iyn[b]) * dwy[b]; + dv_dgy += wx[a] * row_dy; + } + return dv_dgy / dy; + } }; class FilteredNoise { diff --git a/src/qew_model.h b/src/qew_model.h index 62c6674..a182484 100644 --- a/src/qew_model.h +++ b/src/qew_model.h @@ -28,6 +28,25 @@ struct Params { namespace detail { +// Noise access points: prefer the cheap single-quantity samplers when the +// functor provides them (DeviceNoise does); analytic test noises only supply +// the full sample() and fall back to it. Keeps one call site per kernel. +template +KOKKOS_INLINE_FUNCTION real_t noise_value(const Noise &noise, real_t x, real_t y) { + if constexpr (requires { noise.sample_value(x, y); }) + return noise.sample_value(x, y); + else + return noise.sample(x, y).v; +} + +template +KOKKOS_INLINE_FUNCTION real_t noise_dvdy(const Noise &noise, real_t x, real_t y) { + if constexpr (requires { noise.sample_dy(x, y); }) + return noise.sample_dy(x, y); + else + return noise.sample(x, y).dv_dy; +} + // Shared objective body: the line position enters only through the // device-callable functor `at(i)`, so the plain evaluation and the // line-search trial evaluation share one analytic definition. @@ -42,7 +61,7 @@ real_t objective_eval(int n, const Params &p, const Noise &noise, const At &at) Kokkos::parallel_reduce( "qew::objective", n, KOKKOS_LAMBDA(const int i, real_t &acc) { - const int ip = (i + 1) % n; // periodic forward neighbour + const int ip = (i + 1 == n) ? 0 : i + 1; // periodic forward neighbour const real_t hi = at(i); const real_t hip = at(ip); const real_t dh_dx = (hip - hi) / dx; @@ -59,9 +78,7 @@ real_t objective_eval(int n, const Params &p, const Noise &noise, const At &at) // Edge midpoint, matching python xcenter/hcenter (periodic wrap). const real_t xc = static_cast(i + ip) * static_cast(0.5) * dx; const real_t hc = (hip + hi) * static_cast(0.5); - const NoiseSample s = noise.sample(xc, hc); - - acc += dx * (line_energy + s.v - f * hi); + acc += dx * (line_energy + noise_value(noise, xc, hc) - f * hi); }, total); return total; @@ -116,8 +133,8 @@ void gradient(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_for( "qew::gradient", n, KOKKOS_LAMBDA(const int i) { - const int ip = (i + 1) % n; // right neighbour - const int im = (i - 1 + n) % n; // left neighbour + const int ip = (i + 1 == n) ? 0 : i + 1; // right neighbour + const int im = (i == 0) ? n - 1 : i - 1; // left neighbour real_t line_force; if (model == Model::Linear) { @@ -139,11 +156,11 @@ void gradient(const View1D &h, const Params &p, const Noise &noise, const real_t hcl = (h(im) + h(i)) * static_cast(0.5); const real_t xcr = static_cast(i + ip) * static_cast(0.5) * dx; const real_t hcr = (h(ip) + h(i)) * static_cast(0.5); - const NoiseSample sl = noise.sample(xcl, hcl); - const NoiseSample sr = noise.sample(xcr, hcr); + const real_t dvl = detail::noise_dvdy(noise, xcl, hcl); + const real_t dvr = detail::noise_dvdy(noise, xcr, hcr); grad(i) = dx * (line_force + - static_cast(0.5) * (sl.dv_dy + sr.dv_dy) - f); + static_cast(0.5) * (dvl + dvr) - f); }); } diff --git a/src/rosso_krauth.h b/src/rosso_krauth.h index 04719a6..d17cc95 100644 --- a/src/rosso_krauth.h +++ b/src/rosso_krauth.h @@ -54,8 +54,8 @@ KOKKOS_INLINE_FUNCTION real_t rk_site_velocity(int i, int n, real_t dx, real_t l real_t f, Model model, real_t hi, real_t hl, real_t hr, const Noise &noise) { - const int im = (i - 1 + n) % n; - const int ip = (i + 1) % n; + const int im = (i == 0) ? n - 1 : i - 1; + const int ip = (i + 1 == n) ? 0 : i + 1; real_t line_force; if (model == Model::Linear) { @@ -73,8 +73,8 @@ KOKKOS_INLINE_FUNCTION real_t rk_site_velocity(int i, int n, real_t dx, real_t l const real_t hcl = (hl + hi) * static_cast(0.5); const real_t xcr = static_cast(i + ip) * static_cast(0.5) * dx; const real_t hcr = (hi + hr) * static_cast(0.5); - const real_t dvl = noise.sample(xcl, hcl).dv_dy; - const real_t dvr = noise.sample(xcr, hcr).dv_dy; + const real_t dvl = detail::noise_dvdy(noise, xcl, hcl); + const real_t dvr = detail::noise_dvdy(noise, xcr, hcr); const real_t grad = dx * (line_force + static_cast(0.5) * (dvl + dvr) - f); @@ -153,8 +153,8 @@ RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_for( "rk_sweep", n, KOKKOS_LAMBDA(const int i) { if ((i & 1) != color) return; - const int im = (i - 1 + n) % n; - const int ip = (i + 1) % n; + const int im = (i == 0) ? n - 1 : i - 1; + const int ip = (i + 1 == n) ? 0 : i + 1; bool found = false; const real_t hnew = rk_forward_zero(i, n, dx, lt, f, model, h(i), h(im), From 6ba908c3c3c998745867e491543666b6136be465 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 19:12:43 +0000 Subject: [PATCH 11/14] perf: two-pass edge-sharing gradient Each noise edge midpoint was bicubic-gathered twice, once by each adjacent site. gradient() now runs an edge kernel (one dV/dy gather and, for arclength, one normalized slope per edge) into a GradientWorkspace, followed by a cheap per-site combine. Edge values and the combination order are unchanged (FP addition commutes), so results are identical; the golden tests and the unchanged L-BFGS trajectory (2470 iters / 10424 fevals) confirm it. lbfgs/fire/LangevinDynamics own persistent workspaces; a workspace-free overload keeps the old call signature for tests. Serial benchmark (n=262144): grad 33.5 -> 26.5 ms, dyn 52.4 -> 34.4 ms, end-to-end preconditioned L-BFGS (n=4096) 4.77 -> 3.55 s. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/dynamics.h | 3 +- src/fire.h | 5 ++-- src/lbfgs.h | 5 ++-- src/qew_model.h | 76 ++++++++++++++++++++++++++++++++++++------------- 4 files changed, 64 insertions(+), 25 deletions(-) diff --git a/src/dynamics.h b/src/dynamics.h index a4030f9..f430c7a 100644 --- a/src/dynamics.h +++ b/src/dynamics.h @@ -34,7 +34,7 @@ class LangevinDynamics { template void step(const View1D &h, const Params &p, const Noise &noise, real_t drift_coeff, real_t diff_coeff) { - gradient(h, p, noise, grad_); + gradient(h, p, noise, grad_, ws_); const View1D grad = grad_; Pool pool = pool_; const real_t a = drift_coeff; @@ -55,6 +55,7 @@ class LangevinDynamics { private: View1D grad_; + GradientWorkspace ws_; Pool pool_; }; diff --git a/src/fire.h b/src/fire.h index 9c23063..e5ecb0b 100644 --- a/src/fire.h +++ b/src/fire.h @@ -44,6 +44,7 @@ FireResult fire_minimize(const View1D &h, const Params &p, const Noise &noise, const int n = static_cast(h.extent(0)); View1D grad("fire_grad", n); View1D v("fire_v", n); + GradientWorkspace gws; Kokkos::deep_copy(v, static_cast(0)); real_t dt = (fp.dt_init > 0) ? fp.dt_init : static_cast(0.1) * fp.dt_max; @@ -62,7 +63,7 @@ FireResult fire_minimize(const View1D &h, const Params &p, const Noise &noise, return m; }; - gradient(h, p, noise, grad); + gradient(h, p, noise, grad, gws); FireResult res; for (int it = 0; it < fp.max_iter; ++it) { @@ -120,7 +121,7 @@ FireResult fire_minimize(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_for( "fire_xstep", n, KOKKOS_LAMBDA(const int i) { h(i) += dt * v(i); }); - gradient(h, p, noise, grad); + gradient(h, p, noise, grad, gws); } res.iterations = fp.max_iter; diff --git a/src/lbfgs.h b/src/lbfgs.h index 4a5bb3d..a3b2333 100644 --- a/src/lbfgs.h +++ b/src/lbfgs.h @@ -103,6 +103,7 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, View1D g("lbfgs_g", n), g_new("lbfgs_gnew", n), d("lbfgs_d", n); View1D q("lbfgs_q", n), s_tmp("lbfgs_s", n), y_tmp("lbfgs_y", n); + GradientWorkspace gws; std::vector S(m), Y(m); for (int k = 0; k < m; ++k) { S[k] = View1D("lbfgs_S", n); @@ -144,7 +145,7 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, irfft_plan.emplace(exec, qhat, d, KokkosFFT::Direction::backward, -1); } - gradient(h, p, noise, g); + gradient(h, p, noise, g, gws); LbfgsResult res; int ls_fail_streak = 0; real_t phi = 0; // objective at the current iterate, carried across @@ -244,7 +245,7 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, // ---- curvature pair s = alpha d, y = g_new - g ---- // One fused kernel writes the pair and reduces all three scalars. - gradient(h, p, noise, g_new); + gradient(h, p, noise, g_new, gws); real_t sy = 0, ss = 0, yy = 0; Kokkos::parallel_reduce( "lbfgs_sy", n, diff --git a/src/qew_model.h b/src/qew_model.h index a182484..e9e7257 100644 --- a/src/qew_model.h +++ b/src/qew_model.h @@ -118,19 +118,62 @@ real_t objective(const View1D &x0, const View1D &d, real_t alpha, detail::TrialAt{x0, d, alpha}); } +// Scratch for the two-pass gradient: per-edge quantities computed once and +// shared by the two adjacent sites (the dominant bicubic gather would +// otherwise run twice per edge). Own one of these in loops that call +// gradient() repeatedly so the allocations are reused; the workspace-free +// overload below allocates on the fly. +struct GradientWorkspace { + View1D edge_dvdy; // dV/dy at the midpoint of edge (i, i+1) + View1D edge_slope; // arclength only: dh/sqrt(1+dh^2) on edge (i, i+1) + + void ensure(int n, Model model) { + if (static_cast(edge_dvdy.extent(0)) != n) + edge_dvdy = View1D("qew_edge_dvdy", n); + if (model == Model::Arclength && + static_cast(edge_slope.extent(0)) != n) + edge_slope = View1D("qew_edge_slope", n); + } +}; + // Force grad_i = dE/dh_i. Elastic term: harmonic Laplacian (linear) or // slope-normalized difference (arclength); noise term: half the sum of dV/dy at -// the left edge (i-1, i) and right edge (i, i+1). Exact port of qew_gradient, -// run as a parallel_for. `grad` must be pre-allocated with the same extent as h. +// the left edge (i-1, i) and right edge (i, i+1). Exact port of qew_gradient. +// `grad` must be pre-allocated with the same extent as h. +// +// Two passes: edge kernel (one noise gather + one slope per edge), then a +// cheap per-site combine. Site i's left edge is edge i-1, its right edge is +// edge i; the edge midpoints match python xcenter/hcenter (periodic wrap), so +// the result is identical to the former one-pass version that sampled every +// edge twice. template void gradient(const View1D &h, const Params &p, const Noise &noise, - const View1D &grad) { + const View1D &grad, GradientWorkspace &ws) { const int n = static_cast(h.extent(0)); const real_t dx = p.physical_size / static_cast(n); const real_t lt = p.line_tension; const real_t f = p.driving_force; const Model model = p.model; + ws.ensure(n, model); + const View1D edvdy = ws.edge_dvdy; + const View1D eslope = ws.edge_slope; + const bool arc = (model == Model::Arclength); + + Kokkos::parallel_for( + "qew::gradient_edges", n, KOKKOS_LAMBDA(const int i) { + const int ip = (i + 1 == n) ? 0 : i + 1; + const real_t xc = + static_cast(i + ip) * static_cast(0.5) * dx; + const real_t hc = (h(ip) + h(i)) * static_cast(0.5); + edvdy(i) = detail::noise_dvdy(noise, xc, hc); + if (arc) { + const real_t dh = (h(ip) - h(i)) / dx; + eslope(i) = + dh / Kokkos::sqrt(static_cast(1) + dh * dh); + } + }); + Kokkos::parallel_for( "qew::gradient", n, KOKKOS_LAMBDA(const int i) { const int ip = (i + 1 == n) ? 0 : i + 1; // right neighbour @@ -141,29 +184,22 @@ void gradient(const View1D &h, const Params &p, const Noise &noise, const real_t d2h = h(im) + h(ip) - static_cast(2) * h(i); line_force = -lt * d2h / (dx * dx); } else { - const real_t dhr = (h(ip) - h(i)) / dx; // right slope - const real_t dhl = (h(i) - h(im)) / dx; // left slope - line_force = - lt * - (dhl / Kokkos::sqrt(static_cast(1) + dhl * dhl) - - dhr / Kokkos::sqrt(static_cast(1) + dhr * dhr)) / - dx; + line_force = lt * (eslope(im) - eslope(i)) / dx; } - // dV/dy at the left and right edge midpoints (periodic wrap matches - // python xcenter_left/right and hcenter_left/right). - const real_t xcl = static_cast(im + i) * static_cast(0.5) * dx; - const real_t hcl = (h(im) + h(i)) * static_cast(0.5); - const real_t xcr = static_cast(i + ip) * static_cast(0.5) * dx; - const real_t hcr = (h(ip) + h(i)) * static_cast(0.5); - const real_t dvl = detail::noise_dvdy(noise, xcl, hcl); - const real_t dvr = detail::noise_dvdy(noise, xcr, hcr); - grad(i) = dx * (line_force + - static_cast(0.5) * (dvl + dvr) - f); + static_cast(0.5) * + (edvdy(im) + edvdy(i)) - f); }); } +template +void gradient(const View1D &h, const Params &p, const Noise &noise, + const View1D &grad) { + GradientWorkspace ws; + gradient(h, p, noise, grad, ws); +} + } // namespace qew #endif // QEW_MODEL_H From 5cc8276b7d1ba39c88689028d5b4ba19d1ac276b Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 19:15:45 +0000 Subject: [PATCH 12/14] perf: track Rosso-Krauth mean advance instead of re-reducing h All moves are forward-only, so mean_h = mean_h0 + cumulative_advance/n; the per-sweep mean-h parallel_reduce (one of two reductions per sweep, with sweeps up to 2e5) is replaced by a host-side accumulation of the advance sum that is computed anyway. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/rosso_krauth.h | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/rosso_krauth.h b/src/rosso_krauth.h index d17cc95..87d88c1 100644 --- a/src/rosso_krauth.h +++ b/src/rosso_krauth.h @@ -148,6 +148,7 @@ RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, mean_h0 /= static_cast(n); RkResult res; + real_t cum_adv = 0; // forward-only moves: mean_h = mean_h0 + cum_adv / n for (int sweep = 0; sweep < rp.max_sweeps; ++sweep) { for (int color = 0; color < 2; ++color) { Kokkos::parallel_for( @@ -168,17 +169,14 @@ RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_reduce( "rk_adv_sum", n, KOKKOS_LAMBDA(const int i, real_t &acc) { acc += adv(i); }, total_adv); - real_t mean_h = 0; - Kokkos::parallel_reduce( - "rk_mean", n, - KOKKOS_LAMBDA(const int i, real_t &acc) { acc += h(i); }, mean_h); - mean_h /= static_cast(n); + cum_adv += total_adv; + const real_t mean_adv = cum_adv / static_cast(n); res.sweeps = sweep + 1; res.total_advance = total_adv; - res.mean_h = mean_h; + res.mean_h = mean_h0 + mean_adv; - if (mean_h - mean_h0 > rp.runaway) { // depinned: advanced too far + if (mean_adv > rp.runaway) { // depinned: advanced too far res.runaway = true; res.blocked = false; return res; From c38e37076ecb4a5cb88bdf5effe00aa6f8556fa3 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 19:15:45 +0000 Subject: [PATCH 13/14] perf: chunked RNG-state acquisition in the Langevin step One pool state per contiguous lane (chunk of ~32 sites) instead of a locked get_state/free_state per site, the Kokkos-recommended pattern. Serial CPU timing is unchanged within noise (the kick is Box-Muller math there, with uncontended locks); the motivation is wide-GPU lock serialisation, to be verified on the CUDA runner. The zero-temperature path now skips the RNG branch entirely. Draw order changes, so thermal trajectories differ realization-wise; the equipartition and zero-diffusion determinism tests pass. https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- src/dynamics.h | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/src/dynamics.h b/src/dynamics.h index f430c7a..fde88b9 100644 --- a/src/dynamics.h +++ b/src/dynamics.h @@ -39,18 +39,28 @@ class LangevinDynamics { Pool pool = pool_; const real_t a = drift_coeff; const real_t b = diff_coeff; - const bool thermal = (b > 0); - Kokkos::parallel_for( - "langevin_step", static_cast(h.extent(0)), - KOKKOS_LAMBDA(const int i) { - real_t xi = 0; - if (thermal) { + const int n = static_cast(h.extent(0)); + if (b > 0) { + // One RNG state per lane, each covering a contiguous block: a + // pool acquire per chunk instead of per site (the per-site + // acquire is lock traffic that serialises on wide GPUs). + const int n_lanes = (n < 8192) ? n : 8192; + const int chunk = (n + n_lanes - 1) / n_lanes; + Kokkos::parallel_for( + "langevin_step", n_lanes, KOKKOS_LAMBDA(const int t) { auto gen = pool.get_state(); - xi = static_cast(gen.normal()); + const int lo = t * chunk; + const int hi = (lo + chunk < n) ? lo + chunk : n; + for (int i = lo; i < hi; ++i) + h(i) += -a * grad(i) + + b * static_cast(gen.normal()); pool.free_state(gen); - } - h(i) += -a * grad(i) + b * xi; - }); + }); + } else { // zero temperature: pure gradient descent, no RNG at all + Kokkos::parallel_for( + "langevin_step", n, + KOKKOS_LAMBDA(const int i) { h(i) += -a * grad(i); }); + } } private: From dea8ca61ca72ca06ca6d6097c558a63573dd2beb Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 9 Jun 2026 19:17:15 +0000 Subject: [PATCH 14/14] doc: mark phases 3 and 4 done, record measured speedups https://claude.ai/code/session_01QE8Wd5igjty82gh7GvPuVd --- docs/fix_plan.md | 17 +++++++++++------ docs/optimization_notes.md | 8 ++++---- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/docs/fix_plan.md b/docs/fix_plan.md index 38ef92f..e74d510 100644 --- a/docs/fix_plan.md +++ b/docs/fix_plan.md @@ -1,15 +1,20 @@ # Implementation plan for `optimization_notes.md` -> **Status (June 2026):** Phases 0, 1 and 2 are implemented on this branch. +> **Status (June 2026): all phases 0–4 are implemented on this branch.** > Phase 1 was settled by the paper (arXiv:2410.21838 eq. (10): > `lambda_p/xi = (Gamma/U0)^{2/3}`, varied by scanning Gamma at fixed > disorder): all drivers now take the dimensionless ratio as > `pinning_length_over_xi` and share `Gamma = A (lambda_p/xi)^{3/2}`, with -> per-driver mapping tests; line tensions are numerically unchanged. Measured -> effect of Phase 2 on the benchmark's end-to-end preconditioned L-BFGS -> relaxation (Serial, n = 4096, arclength): 5.73 s / 12893 objective evals → -> 4.77 s / 10424 at an identical 2470 iterations (−17% wall time). Phases 3, -> 4, 5 are open. +> per-driver mapping tests; line tensions are numerically unchanged. +> Measured on the benchmark (Serial): Phase 2 took the end-to-end +> preconditioned L-BFGS relaxation (n = 4096, arclength) from 5.73 s / 12893 +> objective evals to 4.77 s / 10424 at an identical 2470 iterations; Phase 3 +> took it to 3.55 s (38% total) and the hot kernels (n = 262144) from +> grad 33.5 / dyn 52.4 ms to 26.5 / 34.4 ms, with the L-BFGS trajectory +> bit-identical throughout. Phase 4's RK change removes one of two +> reductions per sweep; the Langevin RNG chunking is serial-neutral and its +> GPU benefit still needs the CUDA runner (see below). Phase 5 remains +> won't-fix unless profiling says otherwise. Companion to [`optimization_notes.md`](optimization_notes.md); item numbers below refer to that document. Work is grouped into five phases that are diff --git a/docs/optimization_notes.md b/docs/optimization_notes.md index 32710d7..1136e92 100644 --- a/docs/optimization_notes.md +++ b/docs/optimization_notes.md @@ -3,10 +3,10 @@ Technical notes from a code review of the qEW port (June 2026). Items are ordered by expected impact. -> **Status:** items 1, 2, 3, 4 (bugs) and 5, 6, 7 (L-BFGS overhead) are fixed -> on this branch — see `fix_plan.md` for the phase breakdown and measured -> numbers. Items 8–12 (kernel work reduction, RK, dynamics) and the minor -> items remain open. File references point at the current `main`-line +> **Status:** items 1–12 are addressed on this branch — see `fix_plan.md` +> for the phase breakdown and measured numbers. Item 12's GPU benefit +> (RNG lock relief) still needs verification on the CUDA runner; the minor +> items 13–15 are won't-fix unless profiling says otherwise. File references point at the current `main`-line sources. Nothing here changes physics; everything preserves the documented analytic consistency between `objective()` and `gradient()` and the golden-file contracts.