From f060233e535175b2de796e094c47fa3f695083ee Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 5 Jul 2026 08:09:04 -0400 Subject: [PATCH 1/8] bench: Float64 + two-site DMRG2 parity on both sides of suites 1-2 ITensorMPS exposes only two-site dmrg, so the single-site suite-1 variant is dropped and both suites now compare MPSKit DMRG2 (trscheme = truncrank(chi)) against ITensor's default dmrg. Both sides now run Float64 (MPSKit's ComplexF64 default was a 2-4x BLAS handicap on this real-symmetric Hamiltonian). Sweep budget cut from the 30-sweep placeholder to 10, pending the chi = 256 pilot. Trajectory recorders on both sides additionally capture per-sweep GC seconds and allocated bytes (investigation diagnostics). Co-Authored-By: Claude Fable 5 --- benchmark/README.md | 35 ++++++++------ benchmark/comparisons/itensor/README.md | 46 +++++++++++-------- benchmark/comparisons/itensor/common.jl | 34 ++++++++++---- .../itensor/suite1_dmrg_trivial.jl | 22 +++++---- .../comparisons/itensor/suite2_dmrg_u1.jl | 20 +++++--- benchmark/suites/common.jl | 31 +++++++++++-- benchmark/suites/suite1_dmrg_trivial.jl | 35 +++++++++----- benchmark/suites/suite2_dmrg_u1.jl | 24 +++++----- 8 files changed, 162 insertions(+), 85 deletions(-) diff --git a/benchmark/README.md b/benchmark/README.md index 69230398f..c468257de 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -8,31 +8,38 @@ This directory holds two things: ## Comparative harness ``` -julia --project=benchmark benchmark/run.jl --smoke # both suites, small sizes, minutes -julia --project=benchmark benchmark/run.jl --suite=1 # full suite 1 (long) +julia --project=benchmark benchmark/run.jl --smoke # suites 1-2, small sizes, minutes +julia --project=benchmark benchmark/run.jl --smoke --suite=5,7 +julia --project=benchmark benchmark/run.jl --suite=1 # full suite 1 (long; full runs belong on Rusty, see slurm/) julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --smoke -julia --project=benchmark benchmark/plot_results.jl # plots from results/*.json +julia --project=benchmark benchmark/plot_results.jl # plots + comparison tables from results/*.json +julia --project=benchmark benchmark/profile_sweep.jl --chi=256 # diagnostic: where one DMRG2 sweep spends its time ``` -Suites (defined in `docs/IMPROVEMENT_PLAN.md` §4.2): +Suites (defined in `docs/IMPROVEMENT_PLAN.md` §4.2; numbering follows the plan): -- **Suite 1** — finite DMRG time-to-accuracy: spin-1 Heisenberg chain, N = 100, no symmetry, χ ∈ {64, 128, 256, 512, 1024}; per-sweep (energy, walltime) trajectories. +- **Suite 1** — finite two-site DMRG time-to-accuracy: spin-1 Heisenberg chain, N = 100, no symmetry, χ ∈ {64, 128, 256, 512, 1024}; per-sweep (energy, walltime, GC, allocations) trajectories. - **Suite 2** — same with U(1) (Sz conservation). +- **Suite 5** — TDVP throughput: global quench from the Néel state, same chain; two-site growth phase then single-site measure phase, per-step (⟨Sz⟩, χ, walltime) trajectories, wall seconds per unit physical time at χ ∈ {64, 128, 256}. +- **Suite 7** — thread scaling: the suite-1 workload at χ = 512, one process per (JULIA_NUM_THREADS × BLAS-threads) grid point (see `slurm/threads_*.sbatch` for the grid driver). -Results land in `results/` as one JSON per run, carrying full metadata (Julia and package versions, thread counts, BLAS vendor/threads, hostname, timestamp) plus the trajectory. +Results land in `results/` as one JSON per run, carrying full metadata (Julia and package versions, thread counts, BLAS vendor/threads, hostname, timestamp) plus the trajectory. `slurm/` holds the Rusty submission scripts (manual submission; see its README). `ANALYSIS.md` is the investigation log — skeleton until real cluster results exist. + +## Methodology decisions (2026-07-05, maintainer) + +- **Investigation first**: this round is about understanding where MPSKit is faster or slower and why; populating `docs/src/benchmarks.md` is a later round. +- **Two-site DMRG on both sides**: ITensorMPS exposes only two-site `dmrg`, so both suites use MPSKit `DMRG2` with `trscheme = truncrank(χ)` against it. No single-site comparison exists. +- **Float64 on both sides**: the Hamiltonian is real-symmetric; matching real arithmetic removes a 2-4x BLAS confound (MPSKit's constructor default is ComplexF64, passed explicitly as Float64 here). +- **Sweep budget 10** (was a 30-sweep placeholder); to be confirmed by a χ = 256 pilot so every χ point plateaus. Per-sweep trajectories are recorded, so time-to-accuracy extraction does not need long post-convergence tails. +- **Fresh random start per χ** (seeded identically): symmetric across libraries and per-χ independent; χ-ramp warm-start protocols are future work. +- **Full runs happen on the Rusty cluster** via the scripts in `slurm/`, submitted manually by the maintainer; local runs are smoke/plumbing only. ## Methodology guardrails (non-negotiable before publishing any number) -- Compare **time-to-accuracy**, never time-per-iteration: sweep semantics differ across libraries (MPSKit runs single-site DMRG here; ITensorMPS runs its recommended two-site variant). -- Identical protocol knobs: same Hamiltonian, χ schedule, truncation policy (fixed χ, cutoff 0, no noise/expansion), and identical `JULIA_NUM_THREADS` / BLAS backend and thread counts on the same machine. +- Compare **time-to-accuracy**, never time-per-iteration: sweep semantics and per-sweep costs differ across libraries even at matched two-site updates. +- Identical protocol knobs: same Hamiltonian, χ schedule, truncation policy (rank-χ, cutoff 0, no noise/expansion), element type, and identical `JULIA_NUM_THREADS` / BLAS backend and thread counts on the same machine. - **Sanity gate**: at matched χ, the converged energies of the two libraries must agree to ~1e-8 before any timing from that pair is meaningful. Smoke runs (4 sweeps, tiny χ) are for plumbing only — they need agree only in leading digits. - Competitor scripts follow the competitor's *official documented idiom* (citations in the source comments); where exact parity is impossible, the choice made favors the competitor and is documented in `comparisons/itensor/README.md`. - Publish everything or publish nothing: scripts, environment manifests, raw result JSONs, and plotting code accompany any published table. Losses are reported as prominently as wins. -## Known open items before full runs - -- The U(1) virtual-sector allocation in `suites/suite2_dmrg_u1.jl` is a flagged judgment call (`# REVIEW`); the ITensor smoke run reached a lower energy at matched χ, so the MPSKit-side sector split must be reviewed by the maintainer before suite-2 timings mean anything. -- Full-mode `nsweeps = 30` is a placeholder; tune so every χ point converges past the 1e-8 gate. -- Each χ point currently starts from a fresh random state (no ramp-up schedule) on both sides — symmetric, but not how practitioners run DMRG; revisit together. - TeNPy comparisons are out of scope for this round (`comparisons/tenpy/`). diff --git a/benchmark/comparisons/itensor/README.md b/benchmark/comparisons/itensor/README.md index ac6d7627d..61725d9ce 100644 --- a/benchmark/comparisons/itensor/README.md +++ b/benchmark/comparisons/itensor/README.md @@ -17,7 +17,7 @@ meaningful.** Time-to-accuracy is only comparable if both libraries are solving variational problem to the *same* energy. Before trusting a single wall-time number: 1. Run the MPSKit side and this side at the same N and χ schedule. -2. Compare the converged (final-sweep, largest-χ, full 30-sweep) energies. They must match +2. Compare the converged (final-sweep, largest-χ, full-run) energies. They must match to ~1e-8. If they do not, the Hamiltonian, the symmetry sector, or the χ convention is mismatched and **no timing from that run may be published.** @@ -36,7 +36,7 @@ From the repository root: JULIA_NUM_THREADS=1 julia --project=benchmark/comparisons/itensor \ benchmark/comparisons/itensor/run.jl --smoke -# full run (N=100, χ∈[64,128,256,512,1024], 30 sweeps) +# full run (N=100, χ∈[64,128,256,512,1024], 10 sweeps) JULIA_NUM_THREADS=1 julia --project=benchmark/comparisons/itensor \ benchmark/comparisons/itensor/run.jl @@ -71,10 +71,10 @@ the two result files' metadata before comparing, and document the hardware. | # | Aspect | MPSKit side | ITensor side | Exact parity? | Favors | |---|--------|-------------|--------------|---------------|--------| | 1 | Model | `heisenberg_XXX(…; J=1, spin=1)` = J·Σ Sᵢ·Sⱼ | official `OpSum` with `"Sz"`,`"S+"`,`"S-"` (DMRG tutorial), same H | **exact** | neither | -| 2 | Fixed χ | random full-χ `FiniteMPS`, `alg_expand=nothing`, non-truncating gauge | `random_mps(...; linkdims=χ)` + `maxdim=mindim=χ`, `cutoff=0`, `noise=0` | **exact** | neither | -| 3 | Truncation cutoff | none (fixed χ) | `cutoff=0.0` | **exact** | neither (0 = most favorable to ITensor anyway: never discards weight) | -| 4 | Element type | `ComplexF64` (MPSKit default) | `Float64` (ITensor default / documented idiom) | no | **ITensor** (real is faster; forcing complex would only slow ITensor and isn't its idiom) | -| 5 | DMRG update | suite 1: **single-site** (`DMRG`); suite 2: **two-site** (`DMRG2`, `trscheme=truncrank(χ)`) | **two-site** (ITensor's default `dmrg`, the recommended idiom) | suite 1: no; suite 2: **matched** | suite 1: **ITensor** (two-site is more robust against local minima); suite 2: neither | +| 2 | Fixed χ | `DMRG2` with `trscheme=truncrank(χ)`: every two-site update truncates back to ≤ χ | `random_mps(...; linkdims=χ)` + `maxdim=mindim=χ`, `cutoff=0`, `noise=0` | **matched** | neither | +| 3 | Truncation cutoff | none (`truncrank` is rank-only) | `cutoff=0.0` | **exact** | neither (0 = most favorable to ITensor anyway: never discards weight) | +| 4 | Element type | `Float64` (explicit) | `Float64` (ITensor default / documented idiom) | **exact** | neither | +| 5 | DMRG update | **two-site** (`DMRG2`, `trscheme=truncrank(χ)`), both suites | **two-site** (ITensor's default `dmrg`, the recommended idiom; no single-site `dmrg` exists) | **matched** | neither | | 6 | Subspace expansion / noise | none | `noise=0.0` (matches MPSKit) | **exact** | conservative for ITensor — see note | | 7 | Sweeps | `maxiter=nsweeps`, `tol=0` (runs all sweeps) | `nsweeps` identical, no early stop in observer | **exact** | neither | | 8 | U(1) sector | Sz = 0 sector; small sector-diverse seed, **`DMRG2` distributes χ across sectors itself** | `random_mps(sites, state; linkdims=χ)` with alternating Néel state; **ITensor distributes χ across sectors itself** | **matched** (both libraries choose the split) | neither | @@ -85,19 +85,25 @@ the two result files' metadata before comparing, and document the hardware. - **Element type (row 4).** ITensorMPS's `random_mps` is `Float64` by default and that is the correct, documented choice for this real-symmetric Hamiltonian. MPSKit's constructors - default to `ComplexF64`. Real arithmetic is strictly cheaper, so this **favors ITensor**. - We keep it: forcing ITensor into complex would slow it down and depart from its idiom. - -- **Single- vs two-site DMRG (row 5).** Suite 1's MPSKit protocol uses single-site DMRG; - ITensor's default `dmrg` is two-site, which is also what its DMRG tutorial and FAQ - recommend, so for suite 1 this **favors ITensor** (two-site is generally *more robust* - against fixed-χ local minima). Suite 2 instead uses MPSKit's two-site `DMRG2`, matching - the update scheme on both sides — this replaced an earlier hand-picked static U(1) - sector split that was demonstrably suboptimal (−26.4027 vs −26.8188 at the χ = 8 smoke - point). The two updates have different - per-sweep cost and semantics — precisely why §4.3 mandates comparing **time-to-accuracy, - not time-per-sweep**. `nsweeps` is matched on both sides only as a loop bound, not as a - claim that a sweep costs the same. + default to `ComplexF64`, so the MPSKit side passes `Float64` explicitly. Earlier + revisions ran MPSKit at its complex default — a 2-4x BLAS handicap that muddied the + algorithmic comparison; matched real arithmetic on both sides was a maintainer decision + (2026-07-05). + +- **Update scheme (row 5).** Both suites now use two-site DMRG on both sides: ITensorMPS + exposes only the two-site `dmrg` (its tutorial/FAQ-recommended idiom), so no single-site + comparison exists, and the earlier single-site suite-1 variant was dropped (maintainer + decision, 2026-07-05). Suite 2's `DMRG2` also replaced an earlier hand-picked static + U(1) sector split that was demonstrably suboptimal (−26.4027 vs −26.8188 at the χ = 8 + smoke point). Sweep costs still differ in implementation detail — + §4.3 mandates comparing **time-to-accuracy, not time-per-sweep**; `nsweeps` is matched + on both sides only as a loop bound. + +- **Initial state (rows 2/8).** Both sides start from a random state seeded identically + per χ. Trivial suite: both start at full χ. U(1) suite: ITensor's `random_mps` starts + at full χ with its own per-sector split, MPSKit grows from a small sector-diverse seed + within the first sweeps (`DMRG2` redistributes χ per update). Each library follows its + natural protocol; neither split is hand-picked. - **Noise / subspace expansion (row 6).** ITensor's DMRG FAQ recommends a small, decreasing `noise` schedule to help two-site DMRG escape local minima at fixed χ. We set `noise=0.0` @@ -128,7 +134,7 @@ no-noise, two-site run **plateaued in a local minimum** by sweep 2 (trajectory: −26.8188. This is exactly the fixed-χ local-minimum behavior ITensor's FAQ recommends `noise` for (row 6), amplified by only 4 sweeps. It is a smoke-scale convergence artifact, **not** a Hamiltonian/sector mismatch — note the U(1) χ=8 run reaches −26.8188 (matching the -MPSKit reference) from its Néel-seeded start. The full 30-sweep run at large χ is where the +MPSKit reference) from its Néel-seeded start. The full run at large χ is where the 1e-8 sanity gate must hold; verify it there before publishing any timing. ## Do not publish numbers from an unverified run diff --git a/benchmark/comparisons/itensor/common.jl b/benchmark/comparisons/itensor/common.jl index d9f402306..2b9dbe142 100644 --- a/benchmark/comparisons/itensor/common.jl +++ b/benchmark/comparisons/itensor/common.jl @@ -43,18 +43,29 @@ end # (see `benchmark/suites/common.jl`, `dmrg_trajectory`). mutable struct TrajectoryObserver <: AbstractObserver t0::UInt64 + gc_prev::Base.GC_Num energies::Vector{Float64} walltimes::Vector{Float64} - TrajectoryObserver() = new(time_ns(), Float64[], Float64[]) + gctimes::Vector{Float64} + allocd::Vector{Int} + TrajectoryObserver() = new(time_ns(), Base.gc_num(), Float64[], Float64[], Float64[], Int[]) end # `measure!` is passed `energy`, `sweep`, `bond`, `sweep_is_done`, `half_sweep`, ... # per the Observer docs. `energy` is real for a real-symmetric Hamiltonian; `real(...)` -# is defensive and matches the MPSKit side, which also stores `real(...)`. +# is defensive and matches the MPSKit side, which also stores `real(...)`. Per-sweep GC +# seconds and allocated bytes are recorded the same way as the MPSKit side +# (`Base.GC_Diff` between sweeps) — a diagnostic for the investigation, not a published +# metric. function ITensorMPS.measure!(o::TrajectoryObserver; kwargs...) if get(kwargs, :sweep_is_done, false) push!(o.energies, real(kwargs[:energy])) push!(o.walltimes, (time_ns() - o.t0) / 1.0e9) + gc_now = Base.gc_num() + diff = Base.GC_Diff(gc_now, o.gc_prev) + push!(o.gctimes, diff.total_time / 1.0e9) + push!(o.allocd, diff.allocd) + o.gc_prev = gc_now end return nothing end @@ -66,18 +77,18 @@ Run ITensorMPS `dmrg` for exactly `nsweeps` sweeps at *fixed* bond dimension `χ the energy and elapsed wall time after every sweep. Returns `(; psi, energy, energies, walltimes)`. -Fixed-χ protocol (parity with the MPSKit run, which uses a random full-χ `FiniteMPS`, -`alg_expand = nothing`, and a non-truncating gauge so χ never changes): +Fixed-χ protocol (parity with the MPSKit run, which uses two-site `DMRG2` with +`trscheme = truncrank(χ)` so every two-site update truncates back to at most χ): * `maxdim = fill(χ, nsweeps)` pins the ceiling to χ on every sweep. * `mindim = fill(χ, nsweeps)` pins the floor to χ, so the bond dimension stays at χ from the first sweep (the initial `psi0` is already a random full-χ MPS) instead of growing adaptively. This is what makes the ITensor trajectory a genuine fixed-χ trajectory comparable to MPSKit's. Both floors are capped by the local Hilbert-space dimension near the chain ends exactly as MPSKit's are. - * `cutoff = 0.0`: no discarded-weight truncation, matching MPSKit's no-cutoff fixed-χ - run. This is also the setting most favorable to ITensor — it never throws away weight - to save time. (Docs: cutoff is "a float ... specifying the truncation error cutoff", - https://docs.itensor.org/ITensorMPS/stable/DMRG.html.) + * `cutoff = 0.0`: no discarded-weight truncation, matching MPSKit's `truncrank(χ)` + (rank-only, no weight threshold). This is also the setting most favorable to ITensor — + it never throws away weight to save time. (Docs: cutoff is "a float ... specifying the + truncation error cutoff", https://docs.itensor.org/ITensorMPS/stable/DMRG.html.) * `noise = 0.0`: no noise term, matching MPSKit (no subspace expansion). The clock (`t0`) is reset immediately before the `dmrg` call so the recorded wall times @@ -89,6 +100,7 @@ function dmrg_trajectory(H, psi0, χ::Int; nsweeps::Int, outputlevel::Int = 0) maxdim = fill(χ, nsweeps) mindim = fill(χ, nsweeps) observer.t0 = time_ns() + observer.gc_prev = Base.gc_num() # dmrg(H, psi0; nsweeps, maxdim, cutoff, observer, ...) per # https://docs.itensor.org/ITensorMPS/stable/DMRG.html energy, psi = dmrg( @@ -96,7 +108,11 @@ function dmrg_trajectory(H, psi0, χ::Int; nsweeps::Int, outputlevel::Int = 0) nsweeps = nsweeps, maxdim = maxdim, mindim = mindim, cutoff = 0.0, noise = 0.0, outputlevel = outputlevel, observer = observer, ) - return (; psi, energy, energies = observer.energies, walltimes = observer.walltimes) + return (; + psi, energy, + energies = observer.energies, walltimes = observer.walltimes, + gctimes = observer.gctimes, allocd = observer.allocd, + ) end """ diff --git a/benchmark/comparisons/itensor/suite1_dmrg_trivial.jl b/benchmark/comparisons/itensor/suite1_dmrg_trivial.jl index a93cda722..a6954c43b 100644 --- a/benchmark/comparisons/itensor/suite1_dmrg_trivial.jl +++ b/benchmark/comparisons/itensor/suite1_dmrg_trivial.jl @@ -10,14 +10,16 @@ # H = MPO(os, sites) # psi0 = random_mps(sites; linkdims = χ) # This builds H = J * Σ_j S_i·S_j, identical to MPSKit's -# `heisenberg_XXX(ComplexF64, Trivial, FiniteChain(N); J, spin = 1)`. +# `heisenberg_XXX(Float64, Trivial, FiniteChain(N); J, spin = 1)`. # -# PARITY NOTE (element type): ITensor uses its native *real* (Float64) arithmetic here, -# which is the documented idiom (`random_mps` "by default has element type Float64") and -# is the correct choice for this real-symmetric Hamiltonian. MPSKit runs its default -# ComplexF64. Real arithmetic is strictly FASTER, so this favors ITensor — forcing complex -# would only slow ITensor down and is not its idiom, so we keep real. Flagged here so the -# choice is explicit. +# PARITY NOTE (element type): both sides use Float64. `random_mps` "by default has element +# type Float64" (documented idiom, correct for this real-symmetric Hamiltonian), and the +# MPSKit side passes Float64 explicitly for exact parity (maintainer decision, 2026-07-05 — +# previously MPSKit ran its ComplexF64 default, a 2-4x BLAS handicap). +# +# PARITY NOTE (update scheme): both sides use two-site DMRG — ITensor's default `dmrg` +# here, `DMRG2(trscheme = truncrank(χ))` on the MPSKit side. ITensorMPS exposes no +# single-site `dmrg`, so this is the only matched comparison. module ITensorSuite1DMRGTrivial using ITensors @@ -81,6 +83,8 @@ function run(; "chi_actual" => chi_actual, "energies" => result.energies, "walltimes" => result.walltimes, + "gctimes" => result.gctimes, + "allocd_bytes" => result.allocd, # ITensor has no direct analogue of MPSKit's Galerkin (subspace) error; # kept as null so the schema matches and plotting code stays shared. "final_galerkin_error" => nothing, @@ -92,9 +96,11 @@ function run(; data = collect_metadata() data["suite"] = "1-dmrg-trivial" - data["description"] = "finite DMRG time-to-accuracy, spin-1 Heisenberg chain, no symmetry (ITensorMPS)" + data["description"] = "finite two-site DMRG time-to-accuracy, spin-1 Heisenberg chain, no symmetry (ITensorMPS)" data["model"] = "heisenberg_XXX" data["symmetry"] = "Trivial" + data["algorithm"] = "dmrg (two-site), maxdim = mindim = chi, cutoff = 0, noise = 0" + data["eltype"] = "Float64" data["N"] = N data["J"] = J data["spin"] = 1 diff --git a/benchmark/comparisons/itensor/suite2_dmrg_u1.jl b/benchmark/comparisons/itensor/suite2_dmrg_u1.jl index 966b11488..535b61501 100644 --- a/benchmark/comparisons/itensor/suite2_dmrg_u1.jl +++ b/benchmark/comparisons/itensor/suite2_dmrg_u1.jl @@ -14,12 +14,14 @@ # uses the alternating "Up"/"Dn" state array for S=1 sites, so this is the documented idiom. # # PARITY / SECTOR NOTE: the AF spin-1 Heisenberg ground state lives in total Sz = 0, which -# the alternating state selects (`flux(psi0) == QN("Sz", 0)`; logged below). Unlike the -# MPSKit side — which must *hand-pick* how to spread χ across U(1) charge sectors (its -# `u1_virtualspace`/`qmax` REVIEW note) — ITensor's `random_mps(sites, state; linkdims)` -# distributes the bond dimension across sectors automatically for the given flux. So the -# per-sector split here is ITensor's own default, chosen by the library, not by us; this is -# the fair, idiomatic choice and cannot be called a strawman. +# the alternating state selects (`flux(psi0) == QN("Sz", 0)`; logged below). Neither side +# hand-picks the per-sector split of χ: ITensor's `random_mps(sites, state; linkdims)` +# distributes the bond dimension across sectors automatically for the given flux, and the +# MPSKit side's two-site `DMRG2` redistributes χ across sectors at every update starting +# from a small sector-diverse seed. Each library chooses its own split — the fair, +# idiomatic setting on both sides. (Residual asymmetry: ITensor starts at full χ, MPSKit +# grows from the seed within the first sweeps; documented in the README, symmetric in +# spirit since both follow their library's natural protocol.) module ITensorSuite2DMRGU1 using ITensors @@ -95,6 +97,8 @@ function run(; "chi_actual" => chi_actual, "energies" => result.energies, "walltimes" => result.walltimes, + "gctimes" => result.gctimes, + "allocd_bytes" => result.allocd, "final_galerkin_error" => nothing, ) ) @@ -104,9 +108,11 @@ function run(; data = collect_metadata() data["suite"] = "2-dmrg-u1" - data["description"] = "finite DMRG time-to-accuracy, spin-1 Heisenberg chain, U(1) symmetry (ITensorMPS)" + data["description"] = "finite two-site DMRG time-to-accuracy, spin-1 Heisenberg chain, U(1) symmetry (ITensorMPS)" data["model"] = "heisenberg_XXX" data["symmetry"] = "U1 (conserve_sz)" + data["algorithm"] = "dmrg (two-site), maxdim = mindim = chi, cutoff = 0, noise = 0" + data["eltype"] = "Float64" data["N"] = N data["J"] = J data["spin"] = 1 diff --git a/benchmark/suites/common.jl b/benchmark/suites/common.jl index b6ac58b22..0046a79c3 100644 --- a/benchmark/suites/common.jl +++ b/benchmark/suites/common.jl @@ -1,8 +1,9 @@ # Shared utilities for the MPSKit competitive-benchmark suites (docs/IMPROVEMENT_PLAN.md §4). # -# These suites measure *time-to-accuracy* for `find_groundstate` with `DMRG`: for a fixed -# bond dimension χ, run a fixed number of sweeps and record the energy and elapsed wall -# time after every sweep via the algorithm's `finalize` callback (see +# These suites measure *time-to-accuracy* for `find_groundstate` with two-site `DMRG2` +# (the update scheme matched to ITensorMPS, which exposes only two-site `dmrg`): for a +# bond-dimension cap χ, run a fixed number of sweeps and record the energy and elapsed +# wall time after every sweep via the algorithm's `finalize` callback (see # `src/algorithms/groundstate/dmrg.jl`, field `finalize`, signature # `finalize(iter, ψ, H, envs) -> (ψ, envs)`). Accuracy is reported downstream (in # `plot_results.jl`) relative to the lowest energy reached across the whole batch of runs @@ -76,21 +77,33 @@ Galerkin-error convergence check practically never fires early — see `find_groundstate!(::AbstractFiniteMPS, H, ::DMRG, envs)`), recording the energy and elapsed wall time after every sweep via the `finalize` hook. +Not used by the comparison suites (ITensorMPS exposes only two-site `dmrg`, so there is +no single-site counterpart to compare against); kept for MPSKit-internal diagnostics +such as `benchmark/profile_sweep.jl`. + Returns `(; ψ, envs, ϵ, energies, walltimes)` where `energies[i]`/`walltimes[i]` are the real part of the energy and the elapsed wall time (seconds) after sweep `i`. """ function dmrg_trajectory(ψ₀, H; nsweeps::Int, tol::Real = 0.0, verbosity::Int = 0) energies = Float64[] walltimes = Float64[] + gctimes = Float64[] + allocd = Int[] t0 = time_ns() + gc_prev = Ref(Base.gc_num()) finalize = function (iter, ψ, H′, envs) push!(energies, real(expectation_value(ψ, H′, envs))) push!(walltimes, (time_ns() - t0) / 1.0e9) + gc_now = Base.gc_num() + diff = Base.GC_Diff(gc_now, gc_prev[]) + push!(gctimes, diff.total_time / 1.0e9) # GC seconds spent during this sweep + push!(allocd, diff.allocd) # bytes allocated during this sweep + gc_prev[] = gc_now return ψ, envs end alg = DMRG(; tol = tol, maxiter = nsweeps, verbosity = verbosity, finalize = finalize) ψ, envs, ϵ = find_groundstate(ψ₀, H, alg) - return (; ψ, envs, ϵ, energies, walltimes) + return (; ψ, envs, ϵ, energies, walltimes, gctimes, allocd) end """ @@ -108,10 +121,18 @@ Returns `(; ψ, envs, ϵ, energies, walltimes)`, same fields as [`dmrg_trajector function dmrg2_trajectory(ψ₀, H; nsweeps::Int, χ::Int, tol::Real = 0.0, verbosity::Int = 0) energies = Float64[] walltimes = Float64[] + gctimes = Float64[] + allocd = Int[] t0 = time_ns() + gc_prev = Ref(Base.gc_num()) finalize = function (iter, ψ, H′, envs) push!(energies, real(expectation_value(ψ, H′, envs))) push!(walltimes, (time_ns() - t0) / 1.0e9) + gc_now = Base.gc_num() + diff = Base.GC_Diff(gc_now, gc_prev[]) + push!(gctimes, diff.total_time / 1.0e9) # GC seconds spent during this sweep + push!(allocd, diff.allocd) # bytes allocated during this sweep + gc_prev[] = gc_now return ψ, envs end alg = DMRG2(; @@ -119,7 +140,7 @@ function dmrg2_trajectory(ψ₀, H; nsweeps::Int, χ::Int, tol::Real = 0.0, verb trscheme = truncrank(χ), finalize = finalize ) ψ, envs, ϵ = find_groundstate(ψ₀, H, alg) - return (; ψ, envs, ϵ, energies, walltimes) + return (; ψ, envs, ϵ, energies, walltimes, gctimes, allocd) end """ diff --git a/benchmark/suites/suite1_dmrg_trivial.jl b/benchmark/suites/suite1_dmrg_trivial.jl index 039b7bc5a..f52cf7d7b 100644 --- a/benchmark/suites/suite1_dmrg_trivial.jl +++ b/benchmark/suites/suite1_dmrg_trivial.jl @@ -1,11 +1,20 @@ # Suite 1 (docs/IMPROVEMENT_PLAN.md §4.2, item 1): finite DMRG time-to-accuracy for the # spin-1 Heisenberg chain, no symmetry. # -# Protocol: for each χ in a schedule, build a random `FiniteMPS` with (plain, ungraded) -# virtual space ℂ^χ and run single-site DMRG for a fixed number of sweeps, recording the -# energy and elapsed wall time after every sweep (see `BenchCommon.dmrg_trajectory`). -# `alg_expand = nothing` (the default) and a non-truncating gauge, so χ never changes -# during the run and each schedule point genuinely probes that bond dimension. +# Protocol: for each χ in a schedule, build a random full-χ `FiniteMPS` with (plain, +# ungraded) virtual space ℂ^χ and run two-site DMRG (`DMRG2`, `trscheme = truncrank(χ)`) +# for a fixed number of sweeps, recording the energy and elapsed wall time after every +# sweep (see `BenchCommon.dmrg2_trajectory`). Every two-site update truncates back to at +# most χ states, so each schedule point genuinely probes that bond dimension. +# +# Two-site is the matched update scheme: ITensorMPS exposes only two-site `dmrg`, so a +# single-site comparison does not exist and the earlier single-site variant of this suite +# was dropped (maintainer decision, 2026-07-05). `BenchCommon.dmrg_trajectory` (single +# site) is kept for MPSKit-internal diagnostics only. +# +# Element type is Float64 on BOTH sides — the Hamiltonian is real-symmetric, real +# arithmetic is each library's best case, and matching it removes a 2-4x BLAS confound +# (maintainer decision, 2026-07-05). module Suite1DMRGTrivial using MPSKit @@ -28,21 +37,21 @@ function run(; seed::Int = 1234, J::Real = 1.0, spin::Real = 1, resultsdir::AbstractString = BenchCommon.results_dir() ) - H = heisenberg_XXX(ComplexF64, Trivial, FiniteChain(N); J = J, spin = spin) + H = heisenberg_XXX(Float64, Trivial, FiniteChain(N); J = J, spin = spin) pspaces = physicalspace(H) # warmup: run the full pipeline once at a tiny size so JIT compilation does not # pollute the first timed trajectory (methodology guardrail §4.3: wall times must # reflect the algorithm, not the compiler) - let Hw = heisenberg_XXX(ComplexF64, Trivial, FiniteChain(6); J = J, spin = spin) - dmrg_trajectory(FiniteMPS(physicalspace(Hw), ℂ^4), Hw; nsweeps = 2) + let Hw = heisenberg_XXX(Float64, Trivial, FiniteChain(6); J = J, spin = spin) + dmrg2_trajectory(FiniteMPS(Float64, physicalspace(Hw), ℂ^4), Hw; nsweeps = 2, χ = 4) end trials = Vector{Dict{String, Any}}() for χ in chis Random.seed!(seed) - ψ₀ = FiniteMPS(pspaces, ℂ^χ) - elapsed = @elapsed result = dmrg_trajectory(ψ₀, H; nsweeps = nsweeps) + ψ₀ = FiniteMPS(Float64, pspaces, ℂ^χ) + elapsed = @elapsed result = dmrg2_trajectory(ψ₀, H; nsweeps = nsweeps, χ = χ) chi_actual = maximum(dim(left_virtualspace(result.ψ, n)) for n in 2:N) @info "suite 1: χ = $χ done" chi_actual final_energy = last(result.energies) total_time = elapsed final_galerkin_error = result.ϵ @@ -53,6 +62,8 @@ function run(; "chi_actual" => chi_actual, "energies" => result.energies, "walltimes" => result.walltimes, + "gctimes" => result.gctimes, + "allocd_bytes" => result.allocd, "final_galerkin_error" => result.ϵ, ) ) @@ -62,9 +73,11 @@ function run(; data = collect_metadata() data["suite"] = "1-dmrg-trivial" - data["description"] = "finite DMRG time-to-accuracy, spin-1 Heisenberg chain, no symmetry" + data["description"] = "finite two-site DMRG time-to-accuracy, spin-1 Heisenberg chain, no symmetry" data["model"] = "heisenberg_XXX" data["symmetry"] = "Trivial" + data["algorithm"] = "DMRG2(trscheme = truncrank(chi))" + data["eltype"] = "Float64" data["N"] = N data["J"] = J data["spin"] = spin diff --git a/benchmark/suites/suite2_dmrg_u1.jl b/benchmark/suites/suite2_dmrg_u1.jl index 18ff943a5..8a4d5a76a 100644 --- a/benchmark/suites/suite2_dmrg_u1.jl +++ b/benchmark/suites/suite2_dmrg_u1.jl @@ -1,7 +1,7 @@ -# Suite 2 (docs/IMPROVEMENT_PLAN.md §4.2, item 2): same recording protocol as suite 1, -# with U(1) (Sz-conservation) symmetry enforced on the spin-1 Heisenberg chain, and -# two-site DMRG (`DMRG2`) instead of single-site so the bond dimension is distributed -# across symmetry sectors automatically (see the sector-allocation note below). +# Suite 2 (docs/IMPROVEMENT_PLAN.md §4.2, item 2): same protocol and algorithm as suite 1 +# (two-site `DMRG2`), with U(1) (Sz-conservation) symmetry enforced on the spin-1 +# Heisenberg chain. The two-site update also distributes the bond dimension across +# symmetry sectors automatically (see the sector-allocation note below). # # API verified against the installed MPSKitModels source # (`~/.julia/packages/MPSKitModels/*/src/models/hamiltonians.jl` and @@ -23,9 +23,8 @@ using Dates include(joinpath(@__DIR__, "common.jl")) using .BenchCommon -# Sector allocation: unlike suite 1 (fixed full-χ space, single-site DMRG), this suite -# uses two-site DMRG (`DMRG2` with `trscheme = truncrank(χ)`), which redistributes the -# bond dimension across U(1) sectors automatically at every two-site update. This mirrors +# Sector allocation: two-site DMRG (`DMRG2` with `trscheme = truncrank(χ)`) redistributes +# the bond dimension across U(1) sectors automatically at every two-site update. This mirrors # ITensor's behavior, whose two-site `dmrg` likewise chooses the per-sector block sizes # itself; a hand-picked static split (tried first) was demonstrably suboptimal # (−26.4027 vs ITensor's −26.8188 at χ = 8 on the N = 20 smoke check). The initial state @@ -44,20 +43,20 @@ function run(; seed::Int = 1234, J::Real = 1.0, spin::Real = 1, resultsdir::AbstractString = BenchCommon.results_dir() ) - H = heisenberg_XXX(ComplexF64, U1Irrep, FiniteChain(N); J = J, spin = spin) + H = heisenberg_XXX(Float64, U1Irrep, FiniteChain(N); J = J, spin = spin) pspaces = physicalspace(H) # warmup: run the full pipeline once at a tiny size so JIT compilation does not # pollute the first timed trajectory (methodology guardrail §4.3: wall times must # reflect the algorithm, not the compiler) - let Hw = heisenberg_XXX(ComplexF64, U1Irrep, FiniteChain(6); J = J, spin = spin) - dmrg2_trajectory(FiniteMPS(physicalspace(Hw), u1_seedspace()), Hw; nsweeps = 2, χ = 4) + let Hw = heisenberg_XXX(Float64, U1Irrep, FiniteChain(6); J = J, spin = spin) + dmrg2_trajectory(FiniteMPS(Float64, physicalspace(Hw), u1_seedspace()), Hw; nsweeps = 2, χ = 4) end trials = Vector{Dict{String, Any}}() for χ in chis Random.seed!(seed) - ψ₀ = FiniteMPS(pspaces, u1_seedspace()) + ψ₀ = FiniteMPS(Float64, pspaces, u1_seedspace()) elapsed = @elapsed result = dmrg2_trajectory(ψ₀, H; nsweeps = nsweeps, χ = χ) chi_actual = maximum(dim(left_virtualspace(result.ψ, n)) for n in 2:N) @@ -69,6 +68,8 @@ function run(; "chi_actual" => chi_actual, "energies" => result.energies, "walltimes" => result.walltimes, + "gctimes" => result.gctimes, + "allocd_bytes" => result.allocd, "final_galerkin_error" => result.ϵ, ) ) @@ -82,6 +83,7 @@ function run(; data["model"] = "heisenberg_XXX" data["symmetry"] = "U1Irrep" data["algorithm"] = "DMRG2(trscheme = truncrank(chi))" + data["eltype"] = "Float64" data["N"] = N data["J"] = J data["spin"] = spin From 11838bc620526e16b9caccf2a72a40e68d5d9a78 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 5 Jul 2026 08:26:46 -0400 Subject: [PATCH 2/8] bench: TDVP-throughput (suite 5) and thread-scaling (suite 7) harnesses MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Suite 5: global quench from the Néel state on the spin-1 Heisenberg chain, two-site growth phase then single-site measure phase, mirrored on both libraries; per-step (Sz, chi, walltime) trajectories and wall-seconds per unit physical time. Krylov exponentiation tolerance matched at 1e-10 on both sides. Smoke gate: mid-chain Sz agrees across libraries to ~1e-13. Suite 7: suite-1 workload at one (N, chi) point, one process per (JULIA_NUM_THREADS x BLAS-threads) grid point, driven externally. Both run.jl entry points gain --suite=5,7 and --blas-threads=n (default 1). Co-Authored-By: Claude Fable 5 --- benchmark/comparisons/itensor/README.md | 50 ++--- benchmark/comparisons/itensor/run.jl | 50 ++++- benchmark/comparisons/itensor/suite5_tdvp.jl | 183 ++++++++++++++++++ .../comparisons/itensor/suite7_threads.jl | 104 ++++++++++ benchmark/run.jl | 54 +++++- benchmark/suites/suite5_tdvp.jl | 175 +++++++++++++++++ benchmark/suites/suite7_threads.jl | 96 +++++++++ 7 files changed, 674 insertions(+), 38 deletions(-) create mode 100644 benchmark/comparisons/itensor/suite5_tdvp.jl create mode 100644 benchmark/comparisons/itensor/suite7_threads.jl create mode 100644 benchmark/suites/suite5_tdvp.jl create mode 100644 benchmark/suites/suite7_threads.jl diff --git a/benchmark/comparisons/itensor/README.md b/benchmark/comparisons/itensor/README.md index 61725d9ce..4b5e46cc2 100644 --- a/benchmark/comparisons/itensor/README.md +++ b/benchmark/comparisons/itensor/README.md @@ -1,10 +1,13 @@ -# ITensorMPS.jl parity benchmarks (suites 1 & 2) +# ITensorMPS.jl parity benchmarks (suites 1, 2, 5, 7) The ITensorMPS.jl side of the MPSKit competitive benchmark described in `docs/IMPROVEMENT_PLAN.md` §4. These scripts reproduce, as closely as the two libraries allow, the exact protocol of the MPSKit side (`benchmark/suites/`, `benchmark/run.jl`): finite DMRG **time-to-accuracy** for the spin-1 Heisenberg chain, no symmetry (suite 1) and -with U(1) Sz conservation (suite 2). +with U(1) Sz conservation (suite 2); **TDVP throughput** for a global quench (suite 5, see +the parity notes in `suite5_tdvp.jl`, including the matched Krylov-exponentiation +tolerance); and **thread scaling** of the suite-1 workload (suite 7, one process per +thread-grid point — see `benchmark/slurm/threads_*.sbatch`). Everything here uses the **official, documented ITensorMPS idiom** so that a reader from the ITensor community cannot call it a strawman. Every nontrivial API call carries a comment @@ -22,8 +25,8 @@ variational problem to the *same* energy. Before trusting a single wall-time num mismatched and **no timing from that run may be published.** At the *smoke* scale (N = 20, χ ∈ {8, 16}, only 4 sweeps) exact agreement is **not** -expected — 4 sweeps is far from convergence, and fixed-χ variational optima differ slightly -across gauges, truncation policy, and single- vs two-site updates. Leading-digit agreement +expected — 4 sweeps is far from convergence, and fixed-χ variational optima depend on the +random start and can plateau in distinct local minima. Leading-digit agreement (all energies ≈ −26.8) is the smoke-scale bar. The 1e-8 gate applies to the full, converged run only. @@ -114,28 +117,31 @@ the two result files' metadata before comparing, and document the hardware. documented here rather than hidden. This mainly matters at very small χ and few sweeps; see the smoke observation below. -## Smoke verification (recorded 2026-07-04, this machine) +## Smoke verification (recorded 2026-07-05, workstation) -`--smoke` (N=20, χ∈{8,16}, 4 sweeps), `JULIA_NUM_THREADS=1`, OpenBLAS, 16 BLAS threads, -Julia 1.12.6, ITensorMPS 0.4.1 / ITensors 0.9.30. Final-sweep energies vs the MPSKit smoke -reference: +`--smoke` (N=20, χ∈{8,16}, 4 sweeps / 4 TDVP measure steps), `JULIA_NUM_THREADS=1`, +BLAS threads = 1, OpenBLAS, Julia 1.12.6, ITensorMPS 0.4.1 / ITensors 0.9.30, both sides +Float64 + two-site DMRG. Final-sweep energies vs the MPSKit smoke reference: | suite | χ | ITensorMPS | MPSKit reference | |Δ| | |-------|---|-----------|------------------|-----| -| 1 (trivial) | 8 | −26.8048 | −26.8188 | 0.0140 ⚑ | -| 1 (trivial) | 16 | −26.8280 | −26.8334 | 0.0053 | -| 2 (U(1)) | 8 | −26.8188 | — | — | -| 2 (U(1)) | 16 | −26.8313 | — | — | - -All energies agree to leading digits (≈ −26.8), as expected at smoke scale. ⚑ The suite-1 -χ=8 point is 0.014 off — just over the ~1e-2 flag threshold — because that fixed-χ=8, -no-noise, two-site run **plateaued in a local minimum** by sweep 2 (trajectory: -−26.7983 → −26.8048 → −26.8048 → −26.8048) while MPSKit's single-site run kept improving to -−26.8188. This is exactly the fixed-χ local-minimum behavior ITensor's FAQ recommends -`noise` for (row 6), amplified by only 4 sweeps. It is a smoke-scale convergence artifact, -**not** a Hamiltonian/sector mismatch — note the U(1) χ=8 run reaches −26.8188 (matching the -MPSKit reference) from its Néel-seeded start. The full run at large χ is where the -1e-8 sanity gate must hold; verify it there before publishing any timing. +| 1 (trivial) | 8 | −26.8048 | −26.8188 | 0.0139 ⚑ | +| 1 (trivial) | 16 | −26.8280 | −26.8373 | 0.0093 | +| 2 (U(1)) | 8 | −26.8188 | −26.8188 | 6.0e−5 | +| 2 (U(1)) | 16 | −26.8313 | −26.8373 | 0.0060 | + +All energies agree to leading digits (≈ −26.8), the smoke-scale bar. ⚑ The suite-1 χ=8 +point: the ITensor run **plateaued in a fixed-χ local minimum** by sweep 2 while MPSKit's +run kept improving to −26.8188 from its different random start. This is exactly the +fixed-χ local-minimum behavior ITensor's FAQ recommends `noise` for (row 6), amplified by +only 4 sweeps — a smoke-scale convergence artifact, **not** a Hamiltonian/sector mismatch +(the U(1) χ=8 pair agrees to 6e−5 on the same Hamiltonian). The full run at large χ is +where the 1e-8 sanity gate must hold; verify it there before publishing any timing. + +**Suite 5 (TDVP) smoke gate**: the mid-chain ⟨Sz⟩ after the full grow+measure trajectory +agrees between the libraries to ~1e−13 at both χ (χ=8: −0.722010909330507 vs +−0.722010909330496; χ=16: −0.660446916466855 vs −0.660446916466839) — the two sides +demonstrably integrate the same quench with the same protocol. ## Do not publish numbers from an unverified run diff --git a/benchmark/comparisons/itensor/run.jl b/benchmark/comparisons/itensor/run.jl index 74be0f548..920a7fe5c 100644 --- a/benchmark/comparisons/itensor/run.jl +++ b/benchmark/comparisons/itensor/run.jl @@ -3,10 +3,14 @@ # flags, same N / χ schedules, results written into the shared `benchmark/results/`. # # Usage (from the repo root): -# julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl [--smoke] [--suite=1|2|1,2] +# julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl \ +# [--smoke] [--suite=1,2,5,7] [--blas-threads=n] # -# --smoke N = 20, χ ∈ [8, 16], 4 sweeps — a few-minute pipeline check, NOT a result. -# --suite=... comma-separated list of suites to run (default: both 1 and 2). +# --smoke N = 20, tiny χ, few steps — a few-minute pipeline check, NOT a result. +# --suite=... comma-separated list of suites to run (default: 1 and 2). +# 1 = finite DMRG (trivial), 2 = finite DMRG (U(1)), 5 = TDVP throughput, +# 7 = thread scaling (one grid point per invocation; see benchmark/slurm/). +# --blas-threads=n BLAS thread count, applied before any suite runs (default: 1). # # IMPORTANT (methodology guardrail §4.3): for the timings to be comparable, this MUST be # run with the SAME JULIA_NUM_THREADS and the SAME BLAS thread count as the MPSKit side @@ -15,6 +19,8 @@ const HERE = @__DIR__ +using LinearAlgebra + function parse_cli(args) smoke = "--smoke" in args suite_idx = findfirst(a -> startswith(a, "--suite="), args) @@ -24,20 +30,29 @@ function parse_cli(args) value = split(args[suite_idx], "=", limit = 2)[2] parse.(Int, split(value, ",")) end - return (; smoke, suites) + blas_idx = findfirst(a -> startswith(a, "--blas-threads="), args) + blas_threads = isnothing(blas_idx) ? 1 : + parse(Int, split(args[blas_idx], "=", limit = 2)[2]) + return (; smoke, suites, blas_threads) end opts = parse_cli(ARGS) resultsdir = normpath(joinpath(HERE, "..", "..", "results")) mkpath(resultsdir) +BLAS.set_num_threads(opts.blas_threads) -println("ITensorMPS benchmark harness — mode: ", opts.smoke ? "smoke" : "full", ", suites: ", opts.suites) +println( + "ITensorMPS benchmark harness — mode: ", opts.smoke ? "smoke" : "full", + ", suites: ", opts.suites, + ", julia threads: ", Threads.nthreads(), + ", blas threads: ", BLAS.get_num_threads(), +) if 1 in opts.suites include(joinpath(HERE, "suite1_dmrg_trivial.jl")) N = opts.smoke ? 20 : 100 chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] - nsweeps = opts.smoke ? 4 : 30 + nsweeps = opts.smoke ? 4 : 10 println("\n--- Suite 1: finite DMRG, no symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path1 = ITensorSuite1DMRGTrivial.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) println("Suite 1 results written to: ", path1) @@ -47,10 +62,31 @@ if 2 in opts.suites include(joinpath(HERE, "suite2_dmrg_u1.jl")) N = opts.smoke ? 20 : 100 chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] - nsweeps = opts.smoke ? 4 : 30 + nsweeps = opts.smoke ? 4 : 10 println("\n--- Suite 2: finite DMRG, U(1) symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path2 = ITensorSuite2DMRGU1.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) println("Suite 2 results written to: ", path2) end +if 5 in opts.suites + include(joinpath(HERE, "suite5_tdvp.jl")) + N = opts.smoke ? 20 : 100 + chis = opts.smoke ? [8, 16] : [64, 128, 256] + dt = 0.05 + nsteps_measure = opts.smoke ? 4 : 40 + println("\n--- Suite 5: TDVP throughput (N=$N, χ ∈ $chis, dt=$dt, $nsteps_measure measure steps) ---") + path5 = ITensorSuite5TDVP.run(; N = N, chis = chis, dt = dt, nsteps_measure = nsteps_measure, resultsdir = resultsdir) + println("Suite 5 results written to: ", path5) +end + +if 7 in opts.suites + include(joinpath(HERE, "suite7_threads.jl")) + N = opts.smoke ? 20 : 100 + chi = opts.smoke ? 16 : 512 + nsweeps = opts.smoke ? 2 : 5 + println("\n--- Suite 7: thread scaling (N=$N, χ=$chi, $nsweeps sweeps, julia=$(Threads.nthreads()), blas=$(opts.blas_threads)) ---") + path7 = ITensorSuite7Threads.run(; N = N, chi = chi, nsweeps = nsweeps, blas_threads = opts.blas_threads, resultsdir = resultsdir) + println("Suite 7 results written to: ", path7) +end + println("\nDone. These files share the schema of the MPSKit-side results in benchmark/results/.") diff --git a/benchmark/comparisons/itensor/suite5_tdvp.jl b/benchmark/comparisons/itensor/suite5_tdvp.jl new file mode 100644 index 000000000..04cb8a596 --- /dev/null +++ b/benchmark/comparisons/itensor/suite5_tdvp.jl @@ -0,0 +1,183 @@ +# Suite 5 (docs/IMPROVEMENT_PLAN.md §4.2, item 5), ITensorMPS.jl side: TDVP throughput. +# Mirrors `benchmark/suites/suite5_tdvp.jl` — global quench of the spin-1 Heisenberg chain +# (no symmetry) from the Néel product state, two-phase protocol: +# phase 1 (grow): two-site TDVP (`nsite = 2`) for `nsteps_grow` steps, growing the +# bond dimension from 1 to the χ cap (`maxdim = χ`, `cutoff = 0.0` — +# rank-only truncation, matching MPSKit's `truncrank(χ)`). +# phase 2 (measure): single-site TDVP (`nsite = 1`) at fixed bond dimensions for +# `nsteps_measure` steps — the throughput measurement. +# +# API verified against the installed ITensorMPS 0.4.1 source: +# * `tdvp(operator, t, init; time_step, nsite, maxdim, cutoff, (sweep_observer!), ...)` +# computes exp(t * operator) * init (src/solvers/tdvp.jl), so REAL-time evolution by +# exp(-i T H) is `t = -im * T`, `time_step = -im * dt` — the exact counterpart of +# MPSKit's `timestep` convention dt′ = -im * dt. +# * one tdvp "sweep" = one time step; `update_observer!(sweep_observer!; state, sweep, +# current_time, ...)` fires once per step and forwards to `measure!` for an +# `AbstractObserver` (src/solvers/alternating_update.jl, src/update_observer.jl). +# * `MPS(ComplexF64, sites, states)` product-state constructor (src/mps.jl:411): the +# state is built complex explicitly, same as the MPSKit side. +# * `expect(psi, "Sz"; sites = mid:mid)` (src/mps.jl:997). +# +# Integrator parity: both libraries implement the standard symmetric second-order +# one-timestep TDVP sweep (ITensorMPS `reverse_step = true` + `order = 2` defaults; +# MPSKit's timestep! does the analogous forward/backward half-step sweeps). +# +module ITensorSuite5TDVP + +using ITensors +using ITensorMPS +using Dates + +include(joinpath(@__DIR__, "common.jl")) +using .ITensorBenchCommon + +""" + heisenberg_opsum(N; J) -> OpSum + +Spin-1 Heisenberg `OpSum`, H = J·Σ_j S_j·S_{j+1} (identical to suites 1-2). +""" +function heisenberg_opsum(N::Int; J::Real = 1.0) + os = OpSum() + for j in 1:(N - 1) + os += J, "Sz", j, "Sz", j + 1 + os += J / 2, "S+", j, "S-", j + 1 + os += J / 2, "S-", j, "S+", j + 1 + end + return os +end + +# Per-step recorder: `update_observer!` forwards its kwargs (state, sweep, current_time, +# ...) to `measure!` for AbstractObserver subtypes. Physical time is reconstructed as +# t_offset + step * dt from the observer's own counter (current_time is complex here), +# matching the MPSKit-side recorder exactly. +mutable struct TdvpTrajectoryObserver <: AbstractObserver + t0::UInt64 + dt::Float64 + t_offset::Float64 + mid::Int + step::Int + records::Vector{Dict{String, Any}} +end +function TdvpTrajectoryObserver(t0, dt, t_offset, mid) + return TdvpTrajectoryObserver(t0, dt, t_offset, mid, 0, Dict{String, Any}[]) +end + +function ITensorMPS.measure!(o::TdvpTrajectoryObserver; kwargs...) + haskey(kwargs, :state) || return nothing # fires only for the per-step (sweep) hook + state = kwargs[:state] + o.step += 1 + push!( + o.records, Dict{String, Any}( + "t" => o.t_offset + o.step * o.dt, + "sz_mid" => real(first(expect(state, "Sz"; sites = o.mid:o.mid))), + "chi" => maxlinkdim(state), + "walltime" => (time_ns() - o.t0) / 1.0e9, + ) + ) + return nothing +end + +""" + tdvp_phase(H, psi, χ, nsite, nsteps, dt, t0, t_offset, mid) -> (psi, records) + +Run one phase (`nsite` ∈ {1, 2}) of `nsteps` steps of size `dt` and return the evolved +state plus the per-step records. The ⟨Sz⟩ measurement runs inside the timed region, as on +the MPSKit side (identical extra work on both sides). +""" +function tdvp_phase(H, psi, χ::Int, nsite::Int, nsteps::Int, dt::Real, t0, t_offset, mid) + nsteps == 0 && return psi, Dict{String, Any}[] + obs = TdvpTrajectoryObserver(t0, Float64(dt), Float64(t_offset), mid) + psi = tdvp( + H, -im * (nsteps * dt), psi; + time_step = -im * dt, nsite = nsite, + maxdim = χ, cutoff = 0.0, outputlevel = 0, + # Krylov-exponentiation tolerance parity: MPSKit's TDVP integrator defaults to + # tol = 1e-10; ITensorMPS's exponentiate_updater would inherit KrylovKit's + # tighter default (1e-12), i.e. strictly MORE Krylov work per site update. + # Matching 1e-10 keeps the local-solver work comparable and favors ITensor + # relative to its own default. (Both sides use KrylovKit.exponentiate with + # krylovdim = 30 underneath.) + updater_kwargs = (; tol = 1.0e-10), + (sweep_observer!) = obs, + ) + return psi, obs.records +end + +""" + run(; N, chis, dt, nsteps_measure, J = 1.0, resultsdir = results_dir()) + +Run the suite-5 protocol with ITensorMPS and write a timestamped JSON result file (shared +schema, `"library" => "ITensorMPS"`) to `resultsdir`. Returns the path. +""" +function run(; + N::Int, chis::AbstractVector{<:Int}, dt::Real, nsteps_measure::Int, + J::Real = 1.0, + resultsdir::AbstractString = results_dir(), + ) + sites = siteinds("S=1", N) + H = MPO(heisenberg_opsum(N; J = J), sites) + neel = [isodd(n) ? "Up" : "Dn" for n in 1:N] + mid = N ÷ 2 + d = 3 # spin-1 + + # warmup: JIT-compile both phases once at a tiny size (methodology guardrail §4.3) + let Nw = 6, sw = siteinds("S=1", Nw) + Hw = MPO(heisenberg_opsum(Nw; J = J), sw) + pw = MPS(ComplexF64, sw, [isodd(n) ? "Up" : "Dn" for n in 1:Nw]) + t0w = time_ns() + pw, _ = tdvp_phase(Hw, pw, 4, 2, 2, dt, t0w, 0.0, 3) + tdvp_phase(Hw, pw, 4, 1, 2, dt, t0w, 2dt, 3) + end + + trials = Vector{Dict{String, Any}}() + for χ in chis + nsteps_grow = ceil(Int, log(d, χ)) + 2 # same growth budget as the MPSKit side + psi = MPS(ComplexF64, sites, neel) # deterministic product state: no RNG + t0 = time_ns() + psi, rec_grow = tdvp_phase(H, psi, χ, 2, nsteps_grow, dt, t0, 0.0, mid) + psi, rec_meas = tdvp_phase(H, psi, χ, 1, nsteps_measure, dt, t0, nsteps_grow * dt, mid) + total = (time_ns() - t0) / 1.0e9 + + meas_wall = last(rec_meas)["walltime"] - last(rec_grow)["walltime"] + throughput = meas_wall / (nsteps_measure * dt) + chi_actual = maxlinkdim(psi) + + @info "suite 5 (ITensorMPS): χ = $χ done" chi_actual nsteps_grow final_sz = last(rec_meas)["sz_mid"] seconds_per_unit_time = throughput total_time = total + + push!( + trials, Dict{String, Any}( + "chi_target" => χ, + "chi_actual" => chi_actual, + "nsteps_grow" => nsteps_grow, + "nsteps_measure" => nsteps_measure, + "seconds_per_unit_time" => throughput, + "trajectory_grow" => rec_grow, + "trajectory_measure" => rec_meas, + ) + ) + end + + data = collect_metadata() + data["suite"] = "5-tdvp" + data["description"] = "TDVP throughput, global quench from Néel state, spin-1 Heisenberg chain, no symmetry (ITensorMPS)" + data["model"] = "heisenberg_XXX" + data["symmetry"] = "Trivial" + data["algorithm"] = "tdvp nsite=2 grow, then nsite=1 measure (maxdim=chi, cutoff=0)" + data["eltype"] = "ComplexF64" + data["N"] = N + data["J"] = J + data["spin"] = 1 + data["dt"] = dt + data["chi_schedule"] = collect(chis) + data["trials"] = trials + + timestamp = Dates.format(Dates.now(), "yyyymmdd-HHMMSS") + path = joinpath(resultsdir, "itensor_suite5_tdvp_N$(N)_$(timestamp).json") + write_results(path, data) + return path +end + +end # module diff --git a/benchmark/comparisons/itensor/suite7_threads.jl b/benchmark/comparisons/itensor/suite7_threads.jl new file mode 100644 index 000000000..82b001701 --- /dev/null +++ b/benchmark/comparisons/itensor/suite7_threads.jl @@ -0,0 +1,104 @@ +# Suite 7 (docs/IMPROVEMENT_PLAN.md §4.2, item 7), ITensorMPS.jl side: thread scaling. +# Mirrors `benchmark/suites/suite7_threads.jl` — the suite-1 workload (finite two-site +# DMRG, spin-1 Heisenberg, no symmetry) at ONE (N, χ) point, launched once per +# (julia-threads x blas-threads) grid point from the outside: +# +# JULIA_NUM_THREADS= julia --project=benchmark/comparisons/itensor \ +# benchmark/comparisons/itensor/run.jl --suite=7 --blas-threads= +# +# NOTE (methodology): the two libraries parallelize differently. On this dense +# (no-symmetry) workload ITensor parallelizes mostly through BLAS threads; MPSKit also +# exploits Julia threads. Scanning the same grid on both sides and reporting where each +# library peaks is the finding itself, not a parity violation. ITensor's block-sparse +# threading (`ITensors.enable_threaded_blocksparse()`) is irrelevant on this dense +# workload and is left at its default. +module ITensorSuite7Threads + +using ITensors +using ITensorMPS +using Random +using Dates +using LinearAlgebra + +include(joinpath(@__DIR__, "common.jl")) +using .ITensorBenchCommon + +""" + heisenberg_opsum(N; J) -> OpSum + +Spin-1 Heisenberg `OpSum`, H = J·Σ_j S_j·S_{j+1} (identical to suites 1-2). +""" +function heisenberg_opsum(N::Int; J::Real = 1.0) + os = OpSum() + for j in 1:(N - 1) + os += J, "Sz", j, "Sz", j + 1 + os += J / 2, "S+", j, "S-", j + 1 + os += J / 2, "S-", j, "S+", j + 1 + end + return os +end + +""" + run(; N, chi, nsweeps, blas_threads, seed = 1234, J = 1.0, resultsdir = results_dir()) + +Run the suite-1 workload once at bond dimension `chi` with BLAS limited to `blas_threads` +threads, and write a timestamped JSON result file whose filename carries the +(julia-threads, blas-threads) pair. Returns the path. +""" +function run(; + N::Int, chi::Int, nsweeps::Int, blas_threads::Int, + seed::Int = 1234, J::Real = 1.0, + resultsdir::AbstractString = results_dir(), + ) + BLAS.set_num_threads(blas_threads) + + sites = siteinds("S=1", N) + H = MPO(heisenberg_opsum(N; J = J), sites) + + # warmup: JIT-compile the pipeline once at a tiny size (methodology guardrail §4.3) + let Nw = 6, sw = siteinds("S=1", Nw) + Hw = MPO(heisenberg_opsum(Nw; J = J), sw) + dmrg_trajectory(Hw, random_mps(sw; linkdims = 4), 4; nsweeps = 2) + end + + Random.seed!(seed) + ψ₀ = random_mps(sites; linkdims = chi) + elapsed = @elapsed result = dmrg_trajectory(H, ψ₀, chi; nsweeps = nsweeps) + + nj = Threads.nthreads() + nb = BLAS.get_num_threads() + @info "suite 7 (ITensorMPS): (julia = $nj, blas = $nb) done" chi final_energy = last(result.energies) total_time = elapsed + + data = collect_metadata() + data["suite"] = "7-threads" + data["description"] = "thread scaling of the suite-1 workload (finite two-site DMRG, spin-1 Heisenberg, no symmetry) (ITensorMPS)" + data["model"] = "heisenberg_XXX" + data["symmetry"] = "Trivial" + data["algorithm"] = "dmrg (two-site), maxdim = mindim = chi, cutoff = 0, noise = 0" + data["eltype"] = "Float64" + data["N"] = N + data["J"] = J + data["spin"] = 1 + data["nsweeps"] = nsweeps + data["seed"] = seed + data["chi_schedule"] = [chi] + data["best_energy"] = minimum(result.energies) + data["trials"] = [ + Dict{String, Any}( + "chi_target" => chi, + "chi_actual" => maxlinkdim(result.psi), + "energies" => result.energies, + "walltimes" => result.walltimes, + "gctimes" => result.gctimes, + "allocd_bytes" => result.allocd, + "final_galerkin_error" => nothing, + ), + ] + + timestamp = Dates.format(Dates.now(), "yyyymmdd-HHMMSS") + path = joinpath(resultsdir, "itensor_suite7_threads_N$(N)_chi$(chi)_j$(nj)_b$(nb)_$(timestamp).json") + write_results(path, data) + return path +end + +end # module diff --git a/benchmark/run.jl b/benchmark/run.jl index fe0737351..aa9ee9078 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -1,18 +1,24 @@ # Entry point for the MPSKit competitive-benchmark harness (docs/IMPROVEMENT_PLAN.md §4). # # Usage: -# julia --project=benchmark benchmark/run.jl [--smoke] [--suite=1|2|1,2] +# julia --project=benchmark benchmark/run.jl [--smoke] [--suite=1,2,5,7] [--blas-threads=n] # # --smoke use a small chi schedule and system size so the whole harness runs in a few # minutes; verifies the pipeline works, NOT a real benchmark result. -# --suite=... comma-separated list of suites to run (default: both 1 and 2). +# --suite=... comma-separated list of suites to run (default: 1 and 2). +# 1 = finite DMRG (trivial), 2 = finite DMRG (U(1)), 5 = TDVP throughput, +# 7 = thread scaling (one grid point per invocation; see benchmark/slurm/). +# --blas-threads=n BLAS thread count, applied before any suite runs (default: 1, so a +# plain invocation is single-threaded BLAS by construction). Suite 7 exists +# to scan this together with JULIA_NUM_THREADS. # -# This script only runs the MPSKit side of suites 1-2. See `benchmark/comparisons/` for -# the (currently empty) competitor scripts, and `benchmark/plot_results.jl` to render the -# results this script produces. +# See `benchmark/comparisons/` for the competitor scripts and `benchmark/plot_results.jl` +# to render the results this script produces. const HERE = @__DIR__ +using LinearAlgebra + function parse_cli(args) smoke = "--smoke" in args suite_idx = findfirst(a -> startswith(a, "--suite="), args) @@ -22,20 +28,29 @@ function parse_cli(args) value = split(args[suite_idx], "=", limit = 2)[2] parse.(Int, split(value, ",")) end - return (; smoke, suites) + blas_idx = findfirst(a -> startswith(a, "--blas-threads="), args) + blas_threads = isnothing(blas_idx) ? 1 : + parse(Int, split(args[blas_idx], "=", limit = 2)[2]) + return (; smoke, suites, blas_threads) end opts = parse_cli(ARGS) resultsdir = joinpath(HERE, "results") mkpath(resultsdir) +BLAS.set_num_threads(opts.blas_threads) -println("MPSKit benchmark harness — mode: ", opts.smoke ? "smoke" : "full", ", suites: ", opts.suites) +println( + "MPSKit benchmark harness — mode: ", opts.smoke ? "smoke" : "full", + ", suites: ", opts.suites, + ", julia threads: ", Threads.nthreads(), + ", blas threads: ", BLAS.get_num_threads(), +) if 1 in opts.suites include(joinpath(HERE, "suites", "suite1_dmrg_trivial.jl")) N = opts.smoke ? 20 : 100 chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] - nsweeps = opts.smoke ? 4 : 30 + nsweeps = opts.smoke ? 4 : 10 println("\n--- Suite 1: finite DMRG, no symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path1 = Suite1DMRGTrivial.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) println("Suite 1 results written to: ", path1) @@ -45,10 +60,31 @@ if 2 in opts.suites include(joinpath(HERE, "suites", "suite2_dmrg_u1.jl")) N = opts.smoke ? 20 : 100 chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] - nsweeps = opts.smoke ? 4 : 30 + nsweeps = opts.smoke ? 4 : 10 println("\n--- Suite 2: finite DMRG, U(1) symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path2 = Suite2DMRGU1.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) println("Suite 2 results written to: ", path2) end +if 5 in opts.suites + include(joinpath(HERE, "suites", "suite5_tdvp.jl")) + N = opts.smoke ? 20 : 100 + chis = opts.smoke ? [8, 16] : [64, 128, 256] + dt = 0.05 + nsteps_measure = opts.smoke ? 4 : 40 + println("\n--- Suite 5: TDVP throughput (N=$N, χ ∈ $chis, dt=$dt, $nsteps_measure measure steps) ---") + path5 = Suite5TDVP.run(; N = N, chis = chis, dt = dt, nsteps_measure = nsteps_measure, resultsdir = resultsdir) + println("Suite 5 results written to: ", path5) +end + +if 7 in opts.suites + include(joinpath(HERE, "suites", "suite7_threads.jl")) + N = opts.smoke ? 20 : 100 + chi = opts.smoke ? 16 : 512 + nsweeps = opts.smoke ? 2 : 5 + println("\n--- Suite 7: thread scaling (N=$N, χ=$chi, $nsweeps sweeps, julia=$(Threads.nthreads()), blas=$(opts.blas_threads)) ---") + path7 = Suite7Threads.run(; N = N, chi = chi, nsweeps = nsweeps, blas_threads = opts.blas_threads, resultsdir = resultsdir) + println("Suite 7 results written to: ", path7) +end + println("\nDone. Render plots with: julia --project=benchmark benchmark/plot_results.jl") diff --git a/benchmark/suites/suite5_tdvp.jl b/benchmark/suites/suite5_tdvp.jl new file mode 100644 index 000000000..b4d168de8 --- /dev/null +++ b/benchmark/suites/suite5_tdvp.jl @@ -0,0 +1,175 @@ +# Suite 5 (docs/IMPROVEMENT_PLAN.md §4.2, item 5): TDVP throughput. Global quench on the +# spin-1 Heisenberg chain, no symmetry: start from the Néel product state and evolve in +# real time at fixed time step, recording per-step (observable, bond dimension, wall time) +# trajectories at a schedule of bond-dimension caps χ. +# +# Two-phase protocol, mirrored exactly on the ITensorMPS side +# (`benchmark/comparisons/itensor/suite5_tdvp.jl`): +# phase 1 (grow): two-site TDVP (`TDVP2`, `trscheme = truncrank(χ)`) for `nsteps_grow` +# steps — a product state has χ = 1 and single-site TDVP cannot grow +# the bond dimension, so two-site updates grow it to the χ cap. +# phase 2 (measure): single-site TDVP (`TDVP`, fixed bond dimensions) for +# `nsteps_measure` steps — the steady-state workload whose wall time +# per unit physical time is the throughput headline. +# +# +# Sign/element-type conventions (verified against src/algorithms/timestep/): +# * `timestep`/`time_evolve` with real dt > 0 evolve by exp(-i dt H) +# (src/algorithms/timestep/integrators.jl: dt′ = -im * dt). +# * Real-time evolution needs a complex state; the in-place `time_evolve!` does NOT +# promote automatically (only the copying `timestep` wrapper does), so the Néel state +# is built as ComplexF64 directly. The ITensor side is complex for the same reason — +# symmetric, unlike the Float64 DMRG suites. +# * `alg.finalize(t, ψ, H, envs) -> (ψ, envs)` fires once per t_span step inside +# `time_evolve!` — that is the per-step recording hook. +module Suite5TDVP + +using MPSKit +using MPSKitModels +using TensorKit +using LinearAlgebra +using Dates + +include(joinpath(@__DIR__, "common.jl")) +using .BenchCommon + +""" + neel_mps(pspace, N) -> FiniteMPS + +Néel-like product state |↑↓↑↓...⟩ (maximal/minimal Sz on alternating sites) as a χ = 1 +`FiniteMPS` with ComplexF64 scalars. The basis indices for "up"/"down" are read off the +diagonal of the `S_z` operator rather than hard-coded, so this cannot silently disagree +with the MPSKitModels basis convention. +""" +function neel_mps(N::Int; spin::Real = 1) + Sz = S_z(ComplexF64, Trivial; spin = spin) + szdiag = real.(diag(reshape(convert(Array, Sz), dim(domain(Sz)), dim(domain(Sz))))) + up, dn = argmax(szdiag), argmin(szdiag) + d = length(szdiag) + pspace = ℂ^d + tensors = map(1:N) do n + data = zeros(ComplexF64, 1, d, 1) + data[1, isodd(n) ? up : dn, 1] = 1 + return TensorMap(data, ℂ^1 ⊗ pspace ← ℂ^1) + end + return FiniteMPS(tensors) +end + +""" + tdvp_trajectory!(ψ, H, alg_factory, t0_ns, nsteps, dt, t_offset, mid, Sz) -> Vector + +Evolve `ψ` in place for `nsteps` steps of size `dt` starting at physical time `t_offset`, +using the algorithm built by `alg_factory(finalize)`. Records one entry per step: +`(t, ⟨Sz⟩ at site mid, max bond dim, cumulative wall seconds since t0_ns)`. The +measurement cost of ⟨Sz⟩ is inside the timed region on both libraries' sides — identical +extra work, negligible (O(χ²) per step vs O(N χ³) per sweep). +""" +function tdvp_trajectory!(ψ, H, alg_factory, t0_ns, nsteps, dt, t_offset, mid, Sz) + records = Vector{Dict{String, Any}}() + nsteps == 0 && return ψ, records + step = 0 + finalize = function (t, ψ′, H′, envs) + step += 1 + push!( + records, Dict{String, Any}( + "t" => t_offset + step * dt, + "sz_mid" => real(expectation_value(ψ′, mid => Sz)), + "chi" => maximum(dim(left_virtualspace(ψ′, n)) for n in 2:length(ψ′)), + "walltime" => (time_ns() - t0_ns) / 1.0e9, + ) + ) + return ψ′, envs + end + t_span = t_offset .+ (0:nsteps) .* dt + # MPSKit.time_evolve! is unexported but public-shaped (the exported `time_evolve` is + # its copying wrapper); the in-place form avoids a per-step state copy + environment + # rebuild that would pollute the throughput measurement. + ψ, = MPSKit.time_evolve!(ψ, H, t_span, alg_factory(finalize)) + return ψ, records +end + +""" + run(; N, chis, dt, nsteps_measure, J = 1.0, spin = 1, + resultsdir = BenchCommon.results_dir()) + +Run the suite-5 protocol and write a timestamped JSON result file to `resultsdir`. +Returns the path to the written file. The growth phase runs `ceil(log_d(χ)) + 2` steps +(enough two-site updates for the bond dimension to reach the cap from χ = 1). +""" +function run(; + N::Int, chis::AbstractVector{<:Int}, dt::Real, nsteps_measure::Int, + J::Real = 1.0, spin::Real = 1, + resultsdir::AbstractString = BenchCommon.results_dir() + ) + H = heisenberg_XXX(ComplexF64, Trivial, FiniteChain(N); J = J, spin = spin) + Sz = S_z(ComplexF64, Trivial; spin = spin) + mid = N ÷ 2 + d = Int(2 * spin + 1) + + alg2(χ) = fin -> TDVP2(; trscheme = truncrank(χ), finalize = fin) + alg1 = fin -> TDVP(; finalize = fin) + + # warmup: JIT-compile both phases once at a tiny size (methodology guardrail §4.3) + let Hw = heisenberg_XXX(ComplexF64, Trivial, FiniteChain(6); J = J, spin = spin) + ψw = neel_mps(6; spin = spin) + t0w = time_ns() + ψw, _ = tdvp_trajectory!(ψw, Hw, alg2(4), t0w, 2, dt, 0.0, 3, Sz) + tdvp_trajectory!(ψw, Hw, alg1, t0w, 2, dt, 2dt, 3, Sz) + end + + trials = Vector{Dict{String, Any}}() + for χ in chis + nsteps_grow = ceil(Int, log(d, χ)) + 2 + ψ = neel_mps(N; spin = spin) # deterministic product state: no RNG in this suite + t0 = time_ns() + ψ, rec_grow = tdvp_trajectory!(ψ, H, alg2(χ), t0, nsteps_grow, dt, 0.0, mid, Sz) + ψ, rec_meas = tdvp_trajectory!(ψ, H, alg1, t0, nsteps_measure, dt, nsteps_grow * dt, mid, Sz) + total = (time_ns() - t0) / 1.0e9 + + # throughput: wall seconds per unit physical time, measured over phase 2 only + meas_wall = last(rec_meas)["walltime"] - last(rec_grow)["walltime"] + throughput = meas_wall / (nsteps_measure * dt) + chi_actual = maximum(dim(left_virtualspace(ψ, n)) for n in 2:N) + + @info "suite 5: χ = $χ done" chi_actual nsteps_grow final_sz = last(rec_meas)["sz_mid"] seconds_per_unit_time = throughput total_time = total + + push!( + trials, Dict{String, Any}( + "chi_target" => χ, + "chi_actual" => chi_actual, + "nsteps_grow" => nsteps_grow, + "nsteps_measure" => nsteps_measure, + "seconds_per_unit_time" => throughput, + "trajectory_grow" => rec_grow, + "trajectory_measure" => rec_meas, + ) + ) + end + + data = collect_metadata() + data["suite"] = "5-tdvp" + data["description"] = "TDVP throughput, global quench from Néel state, spin-1 Heisenberg chain, no symmetry" + data["model"] = "heisenberg_XXX" + data["symmetry"] = "Trivial" + data["algorithm"] = "TDVP2(trscheme = truncrank(chi)) grow, then TDVP measure" + data["eltype"] = "ComplexF64" + data["N"] = N + data["J"] = J + data["spin"] = spin + data["dt"] = dt + data["chi_schedule"] = collect(chis) + data["trials"] = trials + + timestamp = Dates.format(Dates.now(), "yyyymmdd-HHMMSS") + path = joinpath(resultsdir, "suite5_tdvp_N$(N)_$(timestamp).json") + write_results(path, data) + return path +end + +end # module diff --git a/benchmark/suites/suite7_threads.jl b/benchmark/suites/suite7_threads.jl new file mode 100644 index 000000000..e8be7e481 --- /dev/null +++ b/benchmark/suites/suite7_threads.jl @@ -0,0 +1,96 @@ +# Suite 7 (docs/IMPROVEMENT_PLAN.md §4.2, item 7): thread scaling. Runs the suite-1 +# workload (finite two-site DMRG, spin-1 Heisenberg, no symmetry) at ONE representative +# (N, χ) point and records the wall-time trajectory, so that speedup vs thread count can +# be computed across runs. +# +# One process per grid point: `JULIA_NUM_THREADS` is fixed at launch and cannot be changed +# from inside Julia, so the (julia-threads x blas-threads) grid is driven from the OUTSIDE +# (see `benchmark/slurm/` and the README) by launching this suite once per grid point: +# +# JULIA_NUM_THREADS= julia --project=benchmark benchmark/run.jl --suite=7 --blas-threads= +# +# BLAS threads CAN be set at runtime (`LinearAlgebra.BLAS.set_num_threads`), which is how +# `--blas-threads` is applied before the warmup. The result JSON records both counts in +# its metadata (`nthreads_julia`, `nthreads_blas`) — the plotting side groups on those. +# +# NOTE (methodology): MPSKit and ITensorMPS parallelize differently — MPSKit uses Julia +# threads (multithreaded environments / per-sector block operations via TensorKit and +# OhMyThreads) while ITensor's dense path parallelizes mostly through BLAS threads. The +# grid therefore scans BOTH axes on BOTH libraries; observing that the libraries peak at +# different corners of the grid is a finding, not a parity violation. +module Suite7Threads + +using MPSKit +using MPSKitModels +using TensorKit +using Random +using Dates +using LinearAlgebra + +include(joinpath(@__DIR__, "common.jl")) +using .BenchCommon + +""" + run(; N, chi, nsweeps, blas_threads, seed = 1234, J = 1.0, spin = 1, + resultsdir = BenchCommon.results_dir()) + +Run the suite-1 workload once at bond dimension `chi` with `LinearAlgebra.BLAS` limited to +`blas_threads` threads, and write a timestamped JSON result file to `resultsdir` whose +filename carries the (julia-threads, blas-threads) pair. Returns the path. +""" +function run(; + N::Int, chi::Int, nsweeps::Int, blas_threads::Int, + seed::Int = 1234, J::Real = 1.0, spin::Real = 1, + resultsdir::AbstractString = BenchCommon.results_dir() + ) + BLAS.set_num_threads(blas_threads) + + H = heisenberg_XXX(Float64, Trivial, FiniteChain(N); J = J, spin = spin) + pspaces = physicalspace(H) + + # warmup: JIT-compile the pipeline once at a tiny size (methodology guardrail §4.3) + let Hw = heisenberg_XXX(Float64, Trivial, FiniteChain(6); J = J, spin = spin) + dmrg2_trajectory(FiniteMPS(Float64, physicalspace(Hw), ℂ^4), Hw; nsweeps = 2, χ = 4) + end + + Random.seed!(seed) + ψ₀ = FiniteMPS(Float64, pspaces, ℂ^chi) + elapsed = @elapsed result = dmrg2_trajectory(ψ₀, H; nsweeps = nsweeps, χ = chi) + + nj = Threads.nthreads() + nb = BLAS.get_num_threads() + @info "suite 7: (julia = $nj, blas = $nb) done" chi final_energy = last(result.energies) total_time = elapsed + + data = collect_metadata() + data["suite"] = "7-threads" + data["description"] = "thread scaling of the suite-1 workload (finite two-site DMRG, spin-1 Heisenberg, no symmetry)" + data["model"] = "heisenberg_XXX" + data["symmetry"] = "Trivial" + data["algorithm"] = "DMRG2(trscheme = truncrank(chi))" + data["eltype"] = "Float64" + data["N"] = N + data["J"] = J + data["spin"] = spin + data["nsweeps"] = nsweeps + data["seed"] = seed + data["chi_schedule"] = [chi] + data["best_energy"] = minimum(result.energies) + data["trials"] = [ + Dict{String, Any}( + "chi_target" => chi, + "chi_actual" => maximum(dim(left_virtualspace(result.ψ, n)) for n in 2:N), + "energies" => result.energies, + "walltimes" => result.walltimes, + "gctimes" => result.gctimes, + "allocd_bytes" => result.allocd, + "final_galerkin_error" => result.ϵ, + ), + ] + + timestamp = Dates.format(Dates.now(), "yyyymmdd-HHMMSS") + path = joinpath(resultsdir, "suite7_threads_N$(N)_chi$(chi)_j$(nj)_b$(nb)_$(timestamp).json") + write_results(path, data) + return path +end + +end # module From 1b1c7af10d12cf54ff0d57882d0f0aa52202a157 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 5 Jul 2026 08:35:51 -0400 Subject: [PATCH 3/8] bench: Rusty sbatch scripts, comparison/extraction plots, profiling diagnostic - benchmark/slurm/: one manually-submitted sbatch per (workload x library), single exclusive node, pinned-constraint placeholder, plus a README with the sanity-gate checklist. - plot_results.jl: cross-library trajectory overlays, sanity-gate + time-to-accuracy tables (E_ref = best energy across both libraries), TDVP observable-gate/throughput panels, thread-scaling speedup plot. - profile_sweep.jl: flat + tree profile of one DMRG2 sweep, with GC and allocation totals. - ANALYSIS.md: question-driven skeleton, to be filled only from real cluster results. Co-Authored-By: Claude Fable 5 --- .gitignore | 1 + benchmark/ANALYSIS.md | 57 +++++++ benchmark/plot_results.jl | 199 +++++++++++++++++++++++++ benchmark/profile_sweep.jl | 63 ++++++++ benchmark/slurm/README.md | 54 +++++++ benchmark/slurm/dmrg_itensor.sbatch | 20 +++ benchmark/slurm/dmrg_mpskit.sbatch | 22 +++ benchmark/slurm/tdvp_itensor.sbatch | 18 +++ benchmark/slurm/tdvp_mpskit.sbatch | 18 +++ benchmark/slurm/threads_itensor.sbatch | 24 +++ benchmark/slurm/threads_mpskit.sbatch | 26 ++++ 11 files changed, 502 insertions(+) create mode 100644 benchmark/ANALYSIS.md create mode 100644 benchmark/profile_sweep.jl create mode 100644 benchmark/slurm/README.md create mode 100644 benchmark/slurm/dmrg_itensor.sbatch create mode 100644 benchmark/slurm/dmrg_mpskit.sbatch create mode 100644 benchmark/slurm/tdvp_itensor.sbatch create mode 100644 benchmark/slurm/tdvp_mpskit.sbatch create mode 100644 benchmark/slurm/threads_itensor.sbatch create mode 100644 benchmark/slurm/threads_mpskit.sbatch diff --git a/.gitignore b/.gitignore index e2bb51443..c839e704a 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,7 @@ Manifest.toml **/dev/ benchmark/results/* !benchmark/results/.gitkeep +benchmark/slurm/logs/ # DocumenterVitepress / Node build artifacts **/node_modules/ diff --git a/benchmark/ANALYSIS.md b/benchmark/ANALYSIS.md new file mode 100644 index 000000000..7a75725d0 --- /dev/null +++ b/benchmark/ANALYSIS.md @@ -0,0 +1,57 @@ +# MPSKit vs ITensorMPS — investigation notes + +Status: **skeleton — no full-run results yet.** This file is filled in only from real +Rusty result JSONs in `results/` (never from smoke runs, never from memory). Every claim +must cite the result file(s) it comes from. + +## Questions this round must answer + +### Q1. Sanity gates (prerequisite for everything below) + +- [ ] Suites 1-2: do MPSKit and ITensorMPS converge to the same energy (~1e-8 relative) + at every matched χ? (`compare_dmrg_suite` table in `plot_results.jl` output.) +- [ ] Suite 5: do the ⟨Sz⟩(t) trajectories coincide at matched χ? +- [ ] Metadata parity: same node type, same `nthreads_julia` / `nthreads_blas` / + `blas_config` on both sides of every compared pair? + +### Q2. Headline: time-to-accuracy (suites 1-2) + +- [ ] At each χ, which library reaches E_ref + 1e-8·|E_ref| first, and by what factor? +- [ ] Does the ratio change with χ (i.e. is the crossover at small or large χ)? +- [ ] Does the U(1) suite change the picture relative to the trivial suite (block-sparse + overhead vs dense throughput)? +- [ ] Convergence-quality check: does either side plateau above the gate at some χ + (local minimum), and does that correlate with the fixed-χ no-noise protocol? + +### Q3. TDVP throughput (suite 5) + +- [ ] Seconds of wall time per unit physical time vs χ, per library — who wins, and does + the scaling exponent in χ differ (should be ~χ³; deviations are implementation + overhead)? +- [ ] How large is the growth phase (two-site) cost relative to the measure phase? + +### Q4. Thread scaling (suite 7) + +- [ ] Where does each library peak on the (julia-threads × BLAS-threads) grid? +- [ ] MPSKit: how much does pure Julia-threading ((8,1) vs (1,1)) buy on this dense + workload? ITensor: same question for BLAS ((1,8) vs (1,1)). +- [ ] Do mixed settings ((4,4), (8,8)) oversubscribe and regress? + +### Q5. Where does the time go? (diagnostics) + +- [ ] GC fraction and allocation volume per sweep, per library, vs χ (`gctimes` / + `allocd_bytes` fields). Is either side GC-bound at any χ? +- [ ] `profile_sweep.jl` on the interesting χ points (pick after Q2): what fraction of a + sweep is BLAS / SVD / eigsolve vs MPSKit's own environment and permutation code? + +## Findings + +_(empty — populate per question above, with result-file citations, after the Rusty runs.)_ + +## Follow-ups queued for later rounds + +- χ-ramp warm-start protocol (both libraries) instead of fresh random starts per χ. +- SU(2) capability suite (MPSKit-only, vs its own U(1) run). +- TeNPy comparisons; quasi-2D cylinder workload. +- Publishing: populate `docs/src/benchmarks.md` + README chart only after the sanity + gates pass and the maintainer has reviewed the parity tables. diff --git a/benchmark/plot_results.jl b/benchmark/plot_results.jl index 5f7202415..a5147af94 100644 --- a/benchmark/plot_results.jl +++ b/benchmark/plot_results.jl @@ -194,3 +194,202 @@ function plot_all_suite_results(resultsdir::AbstractString = joinpath(@__DIR__, end plot_all_suite_results(resultdir) + +# ============================================================================ +# Cross-library comparison plots and extractions (investigation layer). +# MPSKit result files have no filename prefix; ITensorMPS files carry `itensor_`. +# Every function below no-ops with an @info when one side is missing, so this script can +# always be run on a partially populated results/ directory. +# ============================================================================ + +""" + latest_result(resultsdir, prefix) -> Union{String, Nothing} + +Path of the most recent result file starting with `prefix` (exact prefix match, so +`"suite1"` does not match `"itensor_suite1"`), or `nothing`. +""" +function latest_result(resultsdir::AbstractString, prefix::AbstractString) + files = filter(readdir(resultsdir)) do fname + startswith(fname, prefix) && endswith(fname, ".json") + end + return isempty(files) ? nothing : joinpath(resultsdir, sort(files)[end]) +end + +""" + time_to_accuracy(trial, e_ref; rtol = 1e-8) -> Union{Float64, Nothing} + +First recorded wall time (s) at which the trial's energy satisfies +`E - e_ref <= rtol * |e_ref|`, or `nothing` if the trajectory never gets there. +Signed difference, not absolute: a variational energy below `e_ref` (better than the +reference) also counts as converged. +""" +function time_to_accuracy(trial, e_ref::Float64; rtol::Float64 = 1.0e-8) + energies = Float64.(trial["energies"]) + walltimes = Float64.(trial["walltimes"]) + idx = findfirst(e -> e - e_ref <= rtol * abs(e_ref), energies) + return isnothing(idx) ? nothing : walltimes[idx] +end + +""" + compare_dmrg_suite(resultsdir, prefix, label; rtol = 1e-8) + +For one DMRG suite (`prefix` ∈ {"suite1_dmrg_trivial", "suite2_dmrg_u1"}): overlay the +energy-error trajectories of the latest MPSKit and ITensorMPS runs (E_ref = best energy +across BOTH files — §4.3: the reference is produced by the suite itself), print the +sanity-gate energy differences and the time-to-accuracy table, and save the figure. +""" +function compare_dmrg_suite(resultsdir::AbstractString, prefix::AbstractString, label::AbstractString; rtol::Float64 = 1.0e-8) + mpskit_file = latest_result(resultsdir, prefix) + itensor_file = latest_result(resultsdir, "itensor_" * prefix) + if isnothing(mpskit_file) || isnothing(itensor_file) + @info "Skipping $label comparison: need both a MPSKit and an ITensorMPS result file in $resultsdir." + return nothing + end + sides = ("MPSKit" => JSON.parsefile(mpskit_file), "ITensorMPS" => JSON.parsefile(itensor_file)) + e_ref = minimum(Float64(result["best_energy"]) for (_, result) in sides) + + f = Figure(; size = (800, 550)) + ax = Axis( + f[1, 1]; + xlabel = "wall time (s)", ylabel = "E - E_ref", + xscale = log10, yscale = log10, title = label + ) + linestyles = Dict("MPSKit" => :solid, "ITensorMPS" => :dash) + for (library, result) in sides + for trial in sort(result["trials"]; by = t -> t["chi_target"]) + energies = Float64.(trial["energies"]) + walltimes = Float64.(trial["walltimes"]) + err = max.(energies .- e_ref, eps(Float64)) + scatterlines!( + ax, walltimes, err; + linestyle = linestyles[library], + label = "$library χ = $(trial["chi_target"])" + ) + end + end + axislegend(ax; position = :rt, nbanks = 2, labelsize = 10) + outbase = joinpath(resultsdir, "compare_" * prefix) + save(outbase * ".png", f) + save(outbase * ".svg", f) + println("Wrote comparison plot for $label -> $(outbase).{png,svg}") + + # sanity gate + time-to-accuracy table (stdout is the deliverable here) + println("\n$label — sanity gate and time to E - E_ref <= $rtol * |E_ref| (E_ref = $e_ref):") + println(" χ | final E (MPSKit) | final E (ITensorMPS) | ΔE | t_acc MPSKit (s) | t_acc ITensorMPS (s)") + mpskit_trials = Dict(t["chi_target"] => t for t in sides[1][2]["trials"]) + itensor_trials = Dict(t["chi_target"] => t for t in sides[2][2]["trials"]) + for χ in sort(collect(keys(mpskit_trials))) + haskey(itensor_trials, χ) || continue + tm, ti = mpskit_trials[χ], itensor_trials[χ] + em, ei = Float64(last(tm["energies"])), Float64(last(ti["energies"])) + tam = time_to_accuracy(tm, e_ref; rtol) + tai = time_to_accuracy(ti, e_ref; rtol) + fmt(x) = isnothing(x) ? "not reached" : string(round(x; sigdigits = 4)) + println(" $χ | $em | $ei | $(abs(em - ei)) | $(fmt(tam)) | $(fmt(tai))") + end + return nothing +end + +""" + compare_tdvp_suite(resultsdir; prefix = "suite5_tdvp") + +Suite 5: overlay the ⟨Sz⟩(t) trajectories of both libraries at matched χ (the TDVP +sanity gate — visible divergence means a protocol mismatch) and plot the throughput +(seconds of wall time per unit physical time, measure phase) vs χ. +""" +function compare_tdvp_suite(resultsdir::AbstractString; prefix::AbstractString = "suite5_tdvp") + mpskit_file = latest_result(resultsdir, prefix) + itensor_file = latest_result(resultsdir, "itensor_" * prefix) + if isnothing(mpskit_file) || isnothing(itensor_file) + @info "Skipping TDVP comparison: need both a MPSKit and an ITensorMPS suite-5 result file in $resultsdir." + return nothing + end + sides = ("MPSKit" => JSON.parsefile(mpskit_file), "ITensorMPS" => JSON.parsefile(itensor_file)) + + f = Figure(; size = (1100, 500)) + ax1 = Axis(f[1, 1]; xlabel = "t", ylabel = "⟨Sz⟩ mid-chain", title = "Suite 5 sanity gate: observable trajectories") + ax2 = Axis( + f[1, 2]; + xlabel = "χ", ylabel = "wall seconds per unit time", + xscale = log2, yscale = log10, title = "Suite 5: TDVP throughput (measure phase)" + ) + linestyles = Dict("MPSKit" => :solid, "ITensorMPS" => :dash) + for (library, result) in sides + chis = Float64[] + thoughputs = Float64[] + for trial in sort(result["trials"]; by = t -> t["chi_target"]) + traj = vcat(trial["trajectory_grow"], trial["trajectory_measure"]) + ts = [Float64(r["t"]) for r in traj] + szs = [Float64(r["sz_mid"]) for r in traj] + lines!(ax1, ts, szs; linestyle = linestyles[library], label = "$library χ = $(trial["chi_target"])") + push!(chis, Float64(trial["chi_target"])) + push!(thoughputs, Float64(trial["seconds_per_unit_time"])) + end + scatterlines!(ax2, chis, thoughputs; linestyle = linestyles[library], label = library) + end + axislegend(ax1; position = :rt, nbanks = 2, labelsize = 10) + axislegend(ax2; position = :lt) + outbase = joinpath(resultsdir, "compare_suite5_tdvp") + save(outbase * ".png", f) + save(outbase * ".svg", f) + println("Wrote TDVP comparison plot -> $(outbase).{png,svg}") + return nothing +end + +""" + compare_thread_scaling(resultsdir) + +Suite 7: collect ALL suite-7 result files (one per (library, julia-threads, blas-threads) +grid point, latest per point), and plot the speedup of the total workload wall time +relative to that library's own (1, 1) baseline. +""" +function compare_thread_scaling(resultsdir::AbstractString) + entries = Dict{Tuple{String, Int, Int}, Any}() # (library, nj, nb) => result, latest wins + for fname in sort(filter(f -> occursin("suite7_threads", f) && endswith(f, ".json"), readdir(resultsdir))) + result = JSON.parsefile(joinpath(resultsdir, fname)) + library = get(result, "library", "MPSKit") + entries[(library, Int(result["nthreads_julia"]), Int(result["nthreads_blas"]))] = result + end + if isempty(entries) + @info "Skipping thread-scaling plot: no suite-7 result files in $resultsdir." + return nothing + end + + total_wall(result) = Float64(last(result["trials"][1]["walltimes"])) + + f = Figure(; size = (900, 500)) + ax = Axis( + f[1, 1]; + xlabel = "(julia threads, blas threads)", ylabel = "speedup vs (1,1)", + title = "Suite 7: thread scaling of the suite-1 workload" + ) + grid = sort(unique([(nj, nb) for (_, nj, nb) in keys(entries)])) + ax.xticks = (1:length(grid), ["($nj,$nb)" for (nj, nb) in grid]) + for library in sort(unique([lib for (lib, _, _) in keys(entries)])) + haskey(entries, (library, 1, 1)) || begin + @info "Thread-scaling: no (1,1) baseline for $library; skipping its speedup line." + continue + end + baseline = total_wall(entries[(library, 1, 1)]) + xs = Int[] + ys = Float64[] + for (i, (nj, nb)) in enumerate(grid) + haskey(entries, (library, nj, nb)) || continue + push!(xs, i) + push!(ys, baseline / total_wall(entries[(library, nj, nb)])) + end + scatterlines!(ax, xs, ys; label = library) + end + hlines!(ax, [1.0]; color = :gray, linestyle = :dot) + axislegend(ax; position = :lt) + outbase = joinpath(resultsdir, "compare_suite7_threads") + save(outbase * ".png", f) + save(outbase * ".svg", f) + println("Wrote thread-scaling plot -> $(outbase).{png,svg}") + return nothing +end + +compare_dmrg_suite(resultdir, "suite1_dmrg_trivial", "Suite 1: finite two-site DMRG, no symmetry") +compare_dmrg_suite(resultdir, "suite2_dmrg_u1", "Suite 2: finite two-site DMRG, U(1) symmetry") +compare_tdvp_suite(resultdir) +compare_thread_scaling(resultdir) diff --git a/benchmark/profile_sweep.jl b/benchmark/profile_sweep.jl new file mode 100644 index 000000000..1a687de99 --- /dev/null +++ b/benchmark/profile_sweep.jl @@ -0,0 +1,63 @@ +# Diagnostic: profile WHERE the time goes in one MPSKit DMRG2 sweep of the suite-1 +# workload. MPSKit-side only — this is an investigation tool, not a benchmark; nothing it +# prints is a comparable number. +# +# Usage: +# julia --project=benchmark benchmark/profile_sweep.jl [--chi=256] [--N=100] [--maxdepth=12] +# +# Runs one warm-up sweep (JIT), then profiles a single additional sweep and prints +# 1. a flat hot-function list (depth-independent; the primary view), +# 2. a culled profile tree (deep by default: the interesting frames sit ~10 levels +# below the script's include scaffolding), and +# 3. the @timed totals (wall seconds, GC seconds, GiB allocated) for that sweep. +# Interpreting: the hot leaves are typically BLAS (gemm), the truncated SVD, and +# KrylovKit's eigsolve; how much sits above them in MPSKit's own environment/derivative +# code is the actionable part. + +using LinearAlgebra +using Profile +using Random + +const HERE = @__DIR__ + +function parse_flag(args, name, default) + idx = findfirst(a -> startswith(a, "--$name="), args) + return isnothing(idx) ? default : parse(Int, split(args[idx], "=", limit = 2)[2]) +end + +chi = parse_flag(ARGS, "chi", 256) +N = parse_flag(ARGS, "N", 100) +maxdepth = parse_flag(ARGS, "maxdepth", 30) +BLAS.set_num_threads(parse_flag(ARGS, "blas-threads", 1)) + +include(joinpath(HERE, "suites", "common.jl")) +using .BenchCommon +using MPSKit +using MPSKitModels +using TensorKit + +println("profile_sweep: N = $N, χ = $chi, julia threads = $(Threads.nthreads()), blas threads = $(BLAS.get_num_threads())") + +H = heisenberg_XXX(Float64, Trivial, FiniteChain(N); J = 1.0, spin = 1) +Random.seed!(1234) +ψ₀ = FiniteMPS(Float64, physicalspace(H), ℂ^chi) + +# one unprofiled sweep: JIT-compile everything and settle the state at the fixed point of +# the fixed-χ protocol, so the profiled sweep is a representative steady-state sweep +res = dmrg2_trajectory(ψ₀, H; nsweeps = 1, χ = chi) +ψ = res.ψ + +Profile.clear() +stats = @timed Profile.@profile dmrg2_trajectory(ψ, H; nsweeps = 1, χ = chi) + +nsamples = length(Profile.fetch()) +println("\n=== hot functions, flat view (mincount = 1% of samples) ===") +Profile.print(; format = :flat, sortedby = :count, mincount = max(1, nsamples ÷ 100)) + +println("\n=== profile tree (mincount = 2% of samples, maxdepth = $maxdepth, noisefloor = 2) ===") +Profile.print(; format = :tree, mincount = max(1, nsamples ÷ 50), maxdepth = maxdepth, noisefloor = 2) + +println("\n=== one-sweep totals (χ = $chi) ===") +println("wall seconds: ", round(stats.time; digits = 3)) +println("GC seconds: ", round(stats.gctime; digits = 3), " (", round(100 * stats.gctime / stats.time; digits = 1), "%)") +println("allocated: ", round(stats.bytes / 2^30; digits = 3), " GiB") diff --git a/benchmark/slurm/README.md b/benchmark/slurm/README.md new file mode 100644 index 000000000..0987dbe4d --- /dev/null +++ b/benchmark/slurm/README.md @@ -0,0 +1,54 @@ +# Rusty submission scripts for the comparative benchmarks + +One sbatch script per (workload x library), all single exclusive node. These are +submitted **manually by the maintainer** — nothing here self-submits. + +## Before submitting + +1. **Pin one CPU generation** — edit `--constraint=` (and `--partition=` if needed) in + ALL six scripts to the same value. Timings across different node types are not + comparable, and the metadata in the result JSONs records only hostname, not node type. +2. **Instantiate both environments on a login node** (compute-node registry access is not + guaranteed): + + ```bash + julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' + julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' + ``` + +3. Make sure `julia` is on `PATH` in batch jobs (juliaup or `module load julia` in your + shell init). + +## Submission order (from the repo root) + +Nothing depends on anything else; submit all six whenever. The DMRG pair is the priority: + +```bash +sbatch benchmark/slurm/dmrg_mpskit.sbatch # suites 1-2, ~hours (chi up to 1024) +sbatch benchmark/slurm/dmrg_itensor.sbatch +sbatch benchmark/slurm/tdvp_mpskit.sbatch # suite 5 +sbatch benchmark/slurm/tdvp_itensor.sbatch +sbatch benchmark/slurm/threads_mpskit.sbatch # suite 7, 7 grid points sequentially +sbatch benchmark/slurm/threads_itensor.sbatch +``` + +Results accumulate in `benchmark/results/` (one JSON per run/grid point); logs in +`benchmark/slurm/logs/`. Time limits are generous guesses — if a job times out, the χ +points already completed are still lost (one JSON per *suite*, written at the end), so +bump `--time` rather than resubmitting blindly. + +## Sanity-gate checklist before reading any timing + +Run `julia --project=benchmark benchmark/plot_results.jl` locally on the synced +`results/` directory, then check: + +1. **Energy gate (suites 1-2)**: at each matched χ, final MPSKit and ITensorMPS energies + agree to ~1e-8 (relative). If not, no timing from that pair means anything. +2. **Observable gate (suite 5)**: the ⟨Sz⟩(t) trajectories of the two libraries lie on + top of each other at matched χ (they solve the same Schrödinger equation with the + same integrator family; visible divergence = protocol mismatch, not physics). +3. **Thread parity (all)**: `nthreads_julia`, `nthreads_blas`, and `blas_config` match + between the two libraries' metadata for every compared pair (suite 7 pairs match + point-by-point). +4. **Same hardware**: all result files report hostnames of the same node type. +5. Only then look at wall times. Report losses as prominently as wins (§4.3). diff --git a/benchmark/slurm/dmrg_itensor.sbatch b/benchmark/slurm/dmrg_itensor.sbatch new file mode 100644 index 000000000..ce80326a7 --- /dev/null +++ b/benchmark/slurm/dmrg_itensor.sbatch @@ -0,0 +1,20 @@ +#!/bin/bash +# Suites 1-2 (finite DMRG, trivial + U(1)), ITensorMPS side. Submit from the repo root: +# sbatch benchmark/slurm/dmrg_itensor.sbatch +# Instantiate the env on the login node first: +# julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' +#SBATCH --job-name=itensor-dmrg +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=12:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR}" +mkdir -p benchmark/slurm/logs + +export JULIA_NUM_THREADS=1 + +julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --suite=1,2 --blas-threads=1 diff --git a/benchmark/slurm/dmrg_mpskit.sbatch b/benchmark/slurm/dmrg_mpskit.sbatch new file mode 100644 index 000000000..d30f1f0cd --- /dev/null +++ b/benchmark/slurm/dmrg_mpskit.sbatch @@ -0,0 +1,22 @@ +#!/bin/bash +# Suites 1-2 (finite DMRG, trivial + U(1)), MPSKit side. Submit from the repo root: +# sbatch benchmark/slurm/dmrg_mpskit.sbatch +# Instantiate the env on the login node first (compute nodes may lack registry access): +# julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' +#SBATCH --job-name=mpskit-dmrg +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=12:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR}" +mkdir -p benchmark/slurm/logs + +# Single-threaded headline comparison (methodology guardrail §4.3): identical thread +# settings on both sides. Suite 7 is the multithreaded story. +export JULIA_NUM_THREADS=1 + +julia --project=benchmark benchmark/run.jl --suite=1,2 --blas-threads=1 diff --git a/benchmark/slurm/tdvp_itensor.sbatch b/benchmark/slurm/tdvp_itensor.sbatch new file mode 100644 index 000000000..abf699dc7 --- /dev/null +++ b/benchmark/slurm/tdvp_itensor.sbatch @@ -0,0 +1,18 @@ +#!/bin/bash +# Suite 5 (TDVP throughput), ITensorMPS side. Submit from the repo root: +# sbatch benchmark/slurm/tdvp_itensor.sbatch +#SBATCH --job-name=itensor-tdvp +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=06:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR}" +mkdir -p benchmark/slurm/logs + +export JULIA_NUM_THREADS=1 + +julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --suite=5 --blas-threads=1 diff --git a/benchmark/slurm/tdvp_mpskit.sbatch b/benchmark/slurm/tdvp_mpskit.sbatch new file mode 100644 index 000000000..4acce0119 --- /dev/null +++ b/benchmark/slurm/tdvp_mpskit.sbatch @@ -0,0 +1,18 @@ +#!/bin/bash +# Suite 5 (TDVP throughput), MPSKit side. Submit from the repo root: +# sbatch benchmark/slurm/tdvp_mpskit.sbatch +#SBATCH --job-name=mpskit-tdvp +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=06:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR}" +mkdir -p benchmark/slurm/logs + +export JULIA_NUM_THREADS=1 + +julia --project=benchmark benchmark/run.jl --suite=5 --blas-threads=1 diff --git a/benchmark/slurm/threads_itensor.sbatch b/benchmark/slurm/threads_itensor.sbatch new file mode 100644 index 000000000..e9c6d5ed4 --- /dev/null +++ b/benchmark/slurm/threads_itensor.sbatch @@ -0,0 +1,24 @@ +#!/bin/bash +# Suite 7 (thread scaling), ITensorMPS side. Same structure and grid as +# threads_mpskit.sbatch. Submit from the repo root: +# sbatch benchmark/slurm/threads_itensor.sbatch +#SBATCH --job-name=itensor-threads +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=08:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR}" +mkdir -p benchmark/slurm/logs + +# Keep the grid identical to threads_mpskit.sbatch. +GRID=("1 1" "4 1" "1 4" "4 4" "8 1" "1 8" "8 8") + +for point in "${GRID[@]}"; do + read -r nj nb <<< "${point}" + echo "=== suite 7 (ITensorMPS): julia=${nj} blas=${nb} ===" + JULIA_NUM_THREADS="${nj}" julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --suite=7 --blas-threads="${nb}" +done diff --git a/benchmark/slurm/threads_mpskit.sbatch b/benchmark/slurm/threads_mpskit.sbatch new file mode 100644 index 000000000..7115d25a6 --- /dev/null +++ b/benchmark/slurm/threads_mpskit.sbatch @@ -0,0 +1,26 @@ +#!/bin/bash +# Suite 7 (thread scaling), MPSKit side. One exclusive node; one julia process per +# (JULIA_NUM_THREADS x BLAS-threads) grid point, run sequentially so every point sees an +# otherwise-idle node. Submit from the repo root: +# sbatch benchmark/slurm/threads_mpskit.sbatch +#SBATCH --job-name=mpskit-threads +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=08:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR}" +mkdir -p benchmark/slurm/logs + +# (julia threads, blas threads) grid — JULIA_NUM_THREADS is launch-time-only, hence the +# one-process-per-point loop. Keep the grid identical to threads_itensor.sbatch. +GRID=("1 1" "4 1" "1 4" "4 4" "8 1" "1 8" "8 8") + +for point in "${GRID[@]}"; do + read -r nj nb <<< "${point}" + echo "=== suite 7 (MPSKit): julia=${nj} blas=${nb} ===" + JULIA_NUM_THREADS="${nj}" julia --project=benchmark benchmark/run.jl --suite=7 --blas-threads="${nb}" +done From a9eaa6621bc5e976f61937f46841852a2390a440 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 5 Jul 2026 09:02:15 -0400 Subject: [PATCH 4/8] bench: pilot-informed schedules (chi <= 512, suite-1 6 sweeps, walltime limits) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The chi = 256, N = 100 pilot crossed the 1e-8 relative gate on sweep 1 and moved < 5e-9 absolute over 11 further sweeps, so suite 1 drops to 6 sweeps (suite 2 keeps 10: seed-growth needs ~6 sweeps to reach chi = 512 first). chi = 1024 leaves the default schedule — extrapolated single-threaded cost is days per point; extend once the chi <= 512 numbers justify it. Suite 7 moves to chi = 256 x 3 sweeps so the 7-point thread grid fits a day. sbatch limits sized accordingly (48 h DMRG, 24 h TDVP/threads); pilot also logged an odd/even sweep-cost alternation, queued in ANALYSIS.md Q5. Co-Authored-By: Claude Fable 5 --- benchmark/ANALYSIS.md | 5 +++++ benchmark/README.md | 8 +++++--- benchmark/comparisons/itensor/run.jl | 12 ++++++------ benchmark/run.jl | 12 ++++++------ benchmark/slurm/README.md | 18 ++++++++++++------ benchmark/slurm/dmrg_itensor.sbatch | 2 +- benchmark/slurm/dmrg_mpskit.sbatch | 2 +- benchmark/slurm/tdvp_itensor.sbatch | 2 +- benchmark/slurm/tdvp_mpskit.sbatch | 2 +- benchmark/slurm/threads_itensor.sbatch | 2 +- benchmark/slurm/threads_mpskit.sbatch | 2 +- 11 files changed, 40 insertions(+), 27 deletions(-) diff --git a/benchmark/ANALYSIS.md b/benchmark/ANALYSIS.md index 7a75725d0..ed1d36c34 100644 --- a/benchmark/ANALYSIS.md +++ b/benchmark/ANALYSIS.md @@ -43,6 +43,11 @@ must cite the result file(s) it comes from. `allocd_bytes` fields). Is either side GC-bound at any χ? - [ ] `profile_sweep.jl` on the interesting χ points (pick after Q2): what fraction of a sweep is BLAS / SVD / eigsolve vs MPSKit's own environment and permutation code? +- [ ] Odd/even sweep-cost alternation: the χ = 256 MPSKit pilot (2026-07-05, workstation, + 8 threads) showed alternating expensive/cheap sweeps (≈450 s vs ≈50 s). Convergence + of the local Krylov solves after the state has converged is the likely cause + (cheap sweeps do almost no eigsolve work); confirm with `gctimes`/profile on the + cluster runs and check whether ITensor shows the same pattern. ## Findings diff --git a/benchmark/README.md b/benchmark/README.md index c468257de..ed62aab43 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -18,10 +18,12 @@ julia --project=benchmark benchmark/profile_sweep.jl --chi=256 # diagnostic: wh Suites (defined in `docs/IMPROVEMENT_PLAN.md` §4.2; numbering follows the plan): -- **Suite 1** — finite two-site DMRG time-to-accuracy: spin-1 Heisenberg chain, N = 100, no symmetry, χ ∈ {64, 128, 256, 512, 1024}; per-sweep (energy, walltime, GC, allocations) trajectories. +- **Suite 1** — finite two-site DMRG time-to-accuracy: spin-1 Heisenberg chain, N = 100, no symmetry, χ ∈ {64, 128, 256, 512}; per-sweep (energy, walltime, GC, allocations) trajectories. - **Suite 2** — same with U(1) (Sz conservation). - **Suite 5** — TDVP throughput: global quench from the Néel state, same chain; two-site growth phase then single-site measure phase, per-step (⟨Sz⟩, χ, walltime) trajectories, wall seconds per unit physical time at χ ∈ {64, 128, 256}. -- **Suite 7** — thread scaling: the suite-1 workload at χ = 512, one process per (JULIA_NUM_THREADS × BLAS-threads) grid point (see `slurm/threads_*.sbatch` for the grid driver). +- **Suite 7** — thread scaling: the suite-1 workload at χ = 256, one process per (JULIA_NUM_THREADS × BLAS-threads) grid point (see `slurm/threads_*.sbatch` for the grid driver). + +χ = 1024 is deliberately absent from the default schedule: extrapolating the χ = 256 pilot (below), a single-threaded χ = 1024 point runs for days. Extend the schedule after the χ ≤ 512 numbers justify the node-time. Results land in `results/` as one JSON per run, carrying full metadata (Julia and package versions, thread counts, BLAS vendor/threads, hostname, timestamp) plus the trajectory. `slurm/` holds the Rusty submission scripts (manual submission; see its README). `ANALYSIS.md` is the investigation log — skeleton until real cluster results exist. @@ -30,7 +32,7 @@ Results land in `results/` as one JSON per run, carrying full metadata (Julia an - **Investigation first**: this round is about understanding where MPSKit is faster or slower and why; populating `docs/src/benchmarks.md` is a later round. - **Two-site DMRG on both sides**: ITensorMPS exposes only two-site `dmrg`, so both suites use MPSKit `DMRG2` with `trscheme = truncrank(χ)` against it. No single-site comparison exists. - **Float64 on both sides**: the Hamiltonian is real-symmetric; matching real arithmetic removes a 2-4x BLAS confound (MPSKit's constructor default is ComplexF64, passed explicitly as Float64 here). -- **Sweep budget 10** (was a 30-sweep placeholder); to be confirmed by a χ = 256 pilot so every χ point plateaus. Per-sweep trajectories are recorded, so time-to-accuracy extraction does not need long post-convergence tails. +- **Sweep budget: 6 (suite 1) / 10 (suite 2)**, pilot-backed. The χ = 256, N = 100 pilot (2026-07-05, 12 sweeps, MPSKit side) crossed the 1e-8 relative gate **on sweep 1** from a random full-χ start and changed by < 5e-9 absolute over the remaining 11 sweeps, so 6 sweeps carries a large margin for the trivial suite. Suite 2 keeps 10 because its bond dimension grows from a small seed (reaching χ = 512 takes ~6 two-site sweeps at d = 3) before convergence proper starts. Per-sweep trajectories are recorded, so time-to-accuracy extraction does not need long post-convergence tails. - **Fresh random start per χ** (seeded identically): symmetric across libraries and per-χ independent; χ-ramp warm-start protocols are future work. - **Full runs happen on the Rusty cluster** via the scripts in `slurm/`, submitted manually by the maintainer; local runs are smoke/plumbing only. diff --git a/benchmark/comparisons/itensor/run.jl b/benchmark/comparisons/itensor/run.jl index 920a7fe5c..2f6a235ce 100644 --- a/benchmark/comparisons/itensor/run.jl +++ b/benchmark/comparisons/itensor/run.jl @@ -51,8 +51,8 @@ println( if 1 in opts.suites include(joinpath(HERE, "suite1_dmrg_trivial.jl")) N = opts.smoke ? 20 : 100 - chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] - nsweeps = opts.smoke ? 4 : 10 + chis = opts.smoke ? [8, 16] : [64, 128, 256, 512] + nsweeps = opts.smoke ? 4 : 6 println("\n--- Suite 1: finite DMRG, no symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path1 = ITensorSuite1DMRGTrivial.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) println("Suite 1 results written to: ", path1) @@ -61,7 +61,7 @@ end if 2 in opts.suites include(joinpath(HERE, "suite2_dmrg_u1.jl")) N = opts.smoke ? 20 : 100 - chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] + chis = opts.smoke ? [8, 16] : [64, 128, 256, 512] nsweeps = opts.smoke ? 4 : 10 println("\n--- Suite 2: finite DMRG, U(1) symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path2 = ITensorSuite2DMRGU1.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) @@ -73,7 +73,7 @@ if 5 in opts.suites N = opts.smoke ? 20 : 100 chis = opts.smoke ? [8, 16] : [64, 128, 256] dt = 0.05 - nsteps_measure = opts.smoke ? 4 : 40 + nsteps_measure = opts.smoke ? 4 : 20 println("\n--- Suite 5: TDVP throughput (N=$N, χ ∈ $chis, dt=$dt, $nsteps_measure measure steps) ---") path5 = ITensorSuite5TDVP.run(; N = N, chis = chis, dt = dt, nsteps_measure = nsteps_measure, resultsdir = resultsdir) println("Suite 5 results written to: ", path5) @@ -82,8 +82,8 @@ end if 7 in opts.suites include(joinpath(HERE, "suite7_threads.jl")) N = opts.smoke ? 20 : 100 - chi = opts.smoke ? 16 : 512 - nsweeps = opts.smoke ? 2 : 5 + chi = opts.smoke ? 16 : 256 + nsweeps = opts.smoke ? 2 : 3 println("\n--- Suite 7: thread scaling (N=$N, χ=$chi, $nsweeps sweeps, julia=$(Threads.nthreads()), blas=$(opts.blas_threads)) ---") path7 = ITensorSuite7Threads.run(; N = N, chi = chi, nsweeps = nsweeps, blas_threads = opts.blas_threads, resultsdir = resultsdir) println("Suite 7 results written to: ", path7) diff --git a/benchmark/run.jl b/benchmark/run.jl index aa9ee9078..c16aa50c8 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -49,8 +49,8 @@ println( if 1 in opts.suites include(joinpath(HERE, "suites", "suite1_dmrg_trivial.jl")) N = opts.smoke ? 20 : 100 - chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] - nsweeps = opts.smoke ? 4 : 10 + chis = opts.smoke ? [8, 16] : [64, 128, 256, 512] + nsweeps = opts.smoke ? 4 : 6 println("\n--- Suite 1: finite DMRG, no symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path1 = Suite1DMRGTrivial.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) println("Suite 1 results written to: ", path1) @@ -59,7 +59,7 @@ end if 2 in opts.suites include(joinpath(HERE, "suites", "suite2_dmrg_u1.jl")) N = opts.smoke ? 20 : 100 - chis = opts.smoke ? [8, 16] : [64, 128, 256, 512, 1024] + chis = opts.smoke ? [8, 16] : [64, 128, 256, 512] nsweeps = opts.smoke ? 4 : 10 println("\n--- Suite 2: finite DMRG, U(1) symmetry (N=$N, χ ∈ $chis, $nsweeps sweeps) ---") path2 = Suite2DMRGU1.run(; N = N, chis = chis, nsweeps = nsweeps, resultsdir = resultsdir) @@ -71,7 +71,7 @@ if 5 in opts.suites N = opts.smoke ? 20 : 100 chis = opts.smoke ? [8, 16] : [64, 128, 256] dt = 0.05 - nsteps_measure = opts.smoke ? 4 : 40 + nsteps_measure = opts.smoke ? 4 : 20 println("\n--- Suite 5: TDVP throughput (N=$N, χ ∈ $chis, dt=$dt, $nsteps_measure measure steps) ---") path5 = Suite5TDVP.run(; N = N, chis = chis, dt = dt, nsteps_measure = nsteps_measure, resultsdir = resultsdir) println("Suite 5 results written to: ", path5) @@ -80,8 +80,8 @@ end if 7 in opts.suites include(joinpath(HERE, "suites", "suite7_threads.jl")) N = opts.smoke ? 20 : 100 - chi = opts.smoke ? 16 : 512 - nsweeps = opts.smoke ? 2 : 5 + chi = opts.smoke ? 16 : 256 + nsweeps = opts.smoke ? 2 : 3 println("\n--- Suite 7: thread scaling (N=$N, χ=$chi, $nsweeps sweeps, julia=$(Threads.nthreads()), blas=$(opts.blas_threads)) ---") path7 = Suite7Threads.run(; N = N, chi = chi, nsweeps = nsweeps, blas_threads = opts.blas_threads, resultsdir = resultsdir) println("Suite 7 results written to: ", path7) diff --git a/benchmark/slurm/README.md b/benchmark/slurm/README.md index 0987dbe4d..de2c507a7 100644 --- a/benchmark/slurm/README.md +++ b/benchmark/slurm/README.md @@ -24,18 +24,24 @@ submitted **manually by the maintainer** — nothing here self-submits. Nothing depends on anything else; submit all six whenever. The DMRG pair is the priority: ```bash -sbatch benchmark/slurm/dmrg_mpskit.sbatch # suites 1-2, ~hours (chi up to 1024) +sbatch benchmark/slurm/dmrg_mpskit.sbatch # suites 1-2, chi up to 512 (48 h limit) sbatch benchmark/slurm/dmrg_itensor.sbatch -sbatch benchmark/slurm/tdvp_mpskit.sbatch # suite 5 +sbatch benchmark/slurm/tdvp_mpskit.sbatch # suite 5 (24 h limit) sbatch benchmark/slurm/tdvp_itensor.sbatch -sbatch benchmark/slurm/threads_mpskit.sbatch # suite 7, 7 grid points sequentially +sbatch benchmark/slurm/threads_mpskit.sbatch # suite 7, 7 grid points sequentially (24 h limit) sbatch benchmark/slurm/threads_itensor.sbatch ``` +Walltime calibration: the 2026-07-05 MPSKit pilot (χ = 256, N = 100, two-site, 8 +workstation threads) averaged ≈ 4.5 min/sweep; single-threaded χ = 512 extrapolates to a +few hours per sweep, which sizes the 48 h DMRG limit. χ = 1024 is deliberately out of the +default schedule (days of single-threaded node-time; extend later if the 512 numbers +justify it). + Results accumulate in `benchmark/results/` (one JSON per run/grid point); logs in -`benchmark/slurm/logs/`. Time limits are generous guesses — if a job times out, the χ -points already completed are still lost (one JSON per *suite*, written at the end), so -bump `--time` rather than resubmitting blindly. +`benchmark/slurm/logs/`. If a job times out, the χ points already completed are still +lost (one JSON per *suite*, written at the end), so bump `--time` rather than +resubmitting blindly. ## Sanity-gate checklist before reading any timing diff --git a/benchmark/slurm/dmrg_itensor.sbatch b/benchmark/slurm/dmrg_itensor.sbatch index ce80326a7..5b17dc74a 100644 --- a/benchmark/slurm/dmrg_itensor.sbatch +++ b/benchmark/slurm/dmrg_itensor.sbatch @@ -6,7 +6,7 @@ #SBATCH --job-name=itensor-dmrg #SBATCH --nodes=1 #SBATCH --exclusive -#SBATCH --time=12:00:00 +#SBATCH --time=48:00:00 #SBATCH --partition=ccq # EDIT if needed #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs #SBATCH --output=benchmark/slurm/logs/%x-%j.out diff --git a/benchmark/slurm/dmrg_mpskit.sbatch b/benchmark/slurm/dmrg_mpskit.sbatch index d30f1f0cd..7920ad13d 100644 --- a/benchmark/slurm/dmrg_mpskit.sbatch +++ b/benchmark/slurm/dmrg_mpskit.sbatch @@ -6,7 +6,7 @@ #SBATCH --job-name=mpskit-dmrg #SBATCH --nodes=1 #SBATCH --exclusive -#SBATCH --time=12:00:00 +#SBATCH --time=48:00:00 #SBATCH --partition=ccq # EDIT if needed #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs #SBATCH --output=benchmark/slurm/logs/%x-%j.out diff --git a/benchmark/slurm/tdvp_itensor.sbatch b/benchmark/slurm/tdvp_itensor.sbatch index abf699dc7..6ffcc93c0 100644 --- a/benchmark/slurm/tdvp_itensor.sbatch +++ b/benchmark/slurm/tdvp_itensor.sbatch @@ -4,7 +4,7 @@ #SBATCH --job-name=itensor-tdvp #SBATCH --nodes=1 #SBATCH --exclusive -#SBATCH --time=06:00:00 +#SBATCH --time=24:00:00 #SBATCH --partition=ccq # EDIT if needed #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs #SBATCH --output=benchmark/slurm/logs/%x-%j.out diff --git a/benchmark/slurm/tdvp_mpskit.sbatch b/benchmark/slurm/tdvp_mpskit.sbatch index 4acce0119..1db5f5fe3 100644 --- a/benchmark/slurm/tdvp_mpskit.sbatch +++ b/benchmark/slurm/tdvp_mpskit.sbatch @@ -4,7 +4,7 @@ #SBATCH --job-name=mpskit-tdvp #SBATCH --nodes=1 #SBATCH --exclusive -#SBATCH --time=06:00:00 +#SBATCH --time=24:00:00 #SBATCH --partition=ccq # EDIT if needed #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs #SBATCH --output=benchmark/slurm/logs/%x-%j.out diff --git a/benchmark/slurm/threads_itensor.sbatch b/benchmark/slurm/threads_itensor.sbatch index e9c6d5ed4..1496db986 100644 --- a/benchmark/slurm/threads_itensor.sbatch +++ b/benchmark/slurm/threads_itensor.sbatch @@ -5,7 +5,7 @@ #SBATCH --job-name=itensor-threads #SBATCH --nodes=1 #SBATCH --exclusive -#SBATCH --time=08:00:00 +#SBATCH --time=24:00:00 #SBATCH --partition=ccq # EDIT if needed #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs #SBATCH --output=benchmark/slurm/logs/%x-%j.out diff --git a/benchmark/slurm/threads_mpskit.sbatch b/benchmark/slurm/threads_mpskit.sbatch index 7115d25a6..e0590e16a 100644 --- a/benchmark/slurm/threads_mpskit.sbatch +++ b/benchmark/slurm/threads_mpskit.sbatch @@ -6,7 +6,7 @@ #SBATCH --job-name=mpskit-threads #SBATCH --nodes=1 #SBATCH --exclusive -#SBATCH --time=08:00:00 +#SBATCH --time=24:00:00 #SBATCH --partition=ccq # EDIT if needed #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs #SBATCH --output=benchmark/slurm/logs/%x-%j.out From e9713028d3e27a5c0ebf20282a508757a9849971 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 5 Jul 2026 09:34:47 -0400 Subject: [PATCH 5/8] bench: single-node run_all.sbatch with ThreadPinning core discipline Replace the six per-workload sbatch scripts with one driver so every compared number comes from the same machine by construction: - phases 1-2 (DMRG, TDVP) run the two libraries concurrently, each single-threaded and taskset-pinned to one core on opposite sockets (2 of 96 cores active; no shared core, L3, or memory controller); - phase 3 (thread grid) runs one point at a time on the idle node; - inside every process, pin_cores() uses ThreadPinning.jl to pin Julia threads 1:1 within the affinity mask (pinthreads(:affinitymask)) and gives OpenBLAS workers disjoint mask cores when the pool is real. SMOKE=1 bash benchmark/slurm/run_all.sbatch rehearses the whole driver locally as plain bash; rehearsal passed on a 2-socket workstation with no pinning warnings, and disjoint 1:1 placement was asserted directly. Co-Authored-By: Claude Fable 5 --- benchmark/Project.toml | 1 + benchmark/README.md | 6 +- benchmark/comparisons/itensor/Project.toml | 9 +- benchmark/comparisons/itensor/README.md | 5 +- benchmark/comparisons/itensor/run.jl | 31 +++++++ benchmark/profile_sweep.jl | 6 ++ benchmark/run.jl | 31 +++++++ benchmark/slurm/README.md | 83 ++++++++++------- benchmark/slurm/dmrg_itensor.sbatch | 20 ---- benchmark/slurm/dmrg_mpskit.sbatch | 22 ----- benchmark/slurm/run_all.sbatch | 102 +++++++++++++++++++++ benchmark/slurm/tdvp_itensor.sbatch | 18 ---- benchmark/slurm/tdvp_mpskit.sbatch | 18 ---- benchmark/slurm/threads_itensor.sbatch | 24 ----- benchmark/slurm/threads_mpskit.sbatch | 26 ------ 15 files changed, 227 insertions(+), 175 deletions(-) delete mode 100644 benchmark/slurm/dmrg_itensor.sbatch delete mode 100644 benchmark/slurm/dmrg_mpskit.sbatch create mode 100644 benchmark/slurm/run_all.sbatch delete mode 100644 benchmark/slurm/tdvp_itensor.sbatch delete mode 100644 benchmark/slurm/tdvp_mpskit.sbatch delete mode 100644 benchmark/slurm/threads_itensor.sbatch delete mode 100644 benchmark/slurm/threads_mpskit.sbatch diff --git a/benchmark/Project.toml b/benchmark/Project.toml index d319fee6f..59ba9e163 100644 --- a/benchmark/Project.toml +++ b/benchmark/Project.toml @@ -15,6 +15,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitManifolds = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684" +ThreadPinning = "811555cd-349b-4f26-b7bc-1f208b848042" [sources] MPSKit = {path = ".."} diff --git a/benchmark/README.md b/benchmark/README.md index ed62aab43..ee3d85ab8 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -21,11 +21,11 @@ Suites (defined in `docs/IMPROVEMENT_PLAN.md` §4.2; numbering follows the plan) - **Suite 1** — finite two-site DMRG time-to-accuracy: spin-1 Heisenberg chain, N = 100, no symmetry, χ ∈ {64, 128, 256, 512}; per-sweep (energy, walltime, GC, allocations) trajectories. - **Suite 2** — same with U(1) (Sz conservation). - **Suite 5** — TDVP throughput: global quench from the Néel state, same chain; two-site growth phase then single-site measure phase, per-step (⟨Sz⟩, χ, walltime) trajectories, wall seconds per unit physical time at χ ∈ {64, 128, 256}. -- **Suite 7** — thread scaling: the suite-1 workload at χ = 256, one process per (JULIA_NUM_THREADS × BLAS-threads) grid point (see `slurm/threads_*.sbatch` for the grid driver). +- **Suite 7** — thread scaling: the suite-1 workload at χ = 256, one process per (JULIA_NUM_THREADS × BLAS-threads) grid point (see `slurm/run_all.sbatch` for the grid driver). χ = 1024 is deliberately absent from the default schedule: extrapolating the χ = 256 pilot (below), a single-threaded χ = 1024 point runs for days. Extend the schedule after the χ ≤ 512 numbers justify the node-time. -Results land in `results/` as one JSON per run, carrying full metadata (Julia and package versions, thread counts, BLAS vendor/threads, hostname, timestamp) plus the trajectory. `slurm/` holds the Rusty submission scripts (manual submission; see its README). `ANALYSIS.md` is the investigation log — skeleton until real cluster results exist. +Results land in `results/` as one JSON per run, carrying full metadata (Julia and package versions, thread counts, BLAS vendor/threads, hostname, timestamp) plus the trajectory. `slurm/run_all.sbatch` runs everything as one single-node Rusty job with pinned cores (manual submission; see `slurm/README.md`). `ANALYSIS.md` is the investigation log — skeleton until real cluster results exist. ## Methodology decisions (2026-07-05, maintainer) @@ -34,7 +34,7 @@ Results land in `results/` as one JSON per run, carrying full metadata (Julia an - **Float64 on both sides**: the Hamiltonian is real-symmetric; matching real arithmetic removes a 2-4x BLAS confound (MPSKit's constructor default is ComplexF64, passed explicitly as Float64 here). - **Sweep budget: 6 (suite 1) / 10 (suite 2)**, pilot-backed. The χ = 256, N = 100 pilot (2026-07-05, 12 sweeps, MPSKit side) crossed the 1e-8 relative gate **on sweep 1** from a random full-χ start and changed by < 5e-9 absolute over the remaining 11 sweeps, so 6 sweeps carries a large margin for the trivial suite. Suite 2 keeps 10 because its bond dimension grows from a small seed (reaching χ = 512 takes ~6 two-site sweeps at d = 3) before convergence proper starts. Per-sweep trajectories are recorded, so time-to-accuracy extraction does not need long post-convergence tails. - **Fresh random start per χ** (seeded identically): symmetric across libraries and per-χ independent; χ-ramp warm-start protocols are future work. -- **Full runs happen on the Rusty cluster** via the scripts in `slurm/`, submitted manually by the maintainer; local runs are smoke/plumbing only. +- **Full runs happen on the Rusty cluster** via `slurm/run_all.sbatch` — a single job on a single node so all pairs share the machine — submitted manually by the maintainer; local runs are smoke/plumbing only. ## Methodology guardrails (non-negotiable before publishing any number) diff --git a/benchmark/comparisons/itensor/Project.toml b/benchmark/comparisons/itensor/Project.toml index 77984754a..03928decf 100644 --- a/benchmark/comparisons/itensor/Project.toml +++ b/benchmark/comparisons/itensor/Project.toml @@ -1,16 +1,9 @@ -# Standalone environment for the ITensorMPS.jl side of the MPSKit competitive benchmark. -# -# This directory is deliberately NOT a member of the MPSKit root workspace: the root -# Project.toml lists `[workspace] projects = ["test", "docs", "examples", "benchmark"]` -# explicitly, and `benchmark/comparisons/itensor` is not among them, so this env resolves -# independently and never drags ITensors into the MPSKit manifest (or vice versa). [deps] ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +ThreadPinning = "811555cd-349b-4f26-b7bc-1f208b848042" -# Pinned to the versions used for the recorded runs (also captured per-run in each result -# file's metadata). Loosen deliberately, re-instantiate, and re-verify the sanity gate. [compat] ITensorMPS = "0.4" ITensors = "0.9" diff --git a/benchmark/comparisons/itensor/README.md b/benchmark/comparisons/itensor/README.md index 4b5e46cc2..37e8b3a96 100644 --- a/benchmark/comparisons/itensor/README.md +++ b/benchmark/comparisons/itensor/README.md @@ -7,7 +7,7 @@ finite DMRG **time-to-accuracy** for the spin-1 Heisenberg chain, no symmetry (s with U(1) Sz conservation (suite 2); **TDVP throughput** for a global quench (suite 5, see the parity notes in `suite5_tdvp.jl`, including the matched Krylov-exponentiation tolerance); and **thread scaling** of the suite-1 workload (suite 7, one process per -thread-grid point — see `benchmark/slurm/threads_*.sbatch`). +thread-grid point — see `benchmark/slurm/run_all.sbatch`). Everything here uses the **official, documented ITensorMPS idiom** so that a reader from the ITensor community cannot call it a strawman. Every nontrivial API call carries a comment @@ -39,7 +39,8 @@ From the repository root: JULIA_NUM_THREADS=1 julia --project=benchmark/comparisons/itensor \ benchmark/comparisons/itensor/run.jl --smoke -# full run (N=100, χ∈[64,128,256,512,1024], 10 sweeps) +# full run (N=100, χ∈[64,128,256,512]; suite 1: 6 sweeps, suite 2: 10) — full runs +# belong in the single-node cluster job, benchmark/slurm/run_all.sbatch JULIA_NUM_THREADS=1 julia --project=benchmark/comparisons/itensor \ benchmark/comparisons/itensor/run.jl diff --git a/benchmark/comparisons/itensor/run.jl b/benchmark/comparisons/itensor/run.jl index 2f6a235ce..4e925a6dd 100644 --- a/benchmark/comparisons/itensor/run.jl +++ b/benchmark/comparisons/itensor/run.jl @@ -20,6 +20,7 @@ const HERE = @__DIR__ using LinearAlgebra +using ThreadPinning function parse_cli(args) smoke = "--smoke" in args @@ -36,10 +37,40 @@ function parse_cli(args) return (; smoke, suites, blas_threads) end +""" + pin_cores(blas_threads) + +Pin Julia threads 1:1 to distinct cores inside the affinity mask this process was +launched with (`taskset` in `benchmark/slurm/run_all.sbatch`, or Slurm's cpu mask), via +ThreadPinning's `:affinitymask` policy. When the OpenBLAS pool has real workers +(`blas_threads > 1`), pin those to mask cores disjoint from the Julia threads. +No-oversubscription is the launcher's job: every concurrent process gets a disjoint mask. +Pinning failure is a warning, not an error, so smoke runs on unsupported platforms still +work — but cluster timings must never run unpinned, so check the log. +""" +function pin_cores(blas_threads::Int) + try + pinthreads(:affinitymask) + if blas_threads > 1 + mask_ids = findall(==(1), getaffinity(; cutoff = nothing)) .- 1 + free_ids = setdiff(mask_ids, getcpuids()) + if length(free_ids) < blas_threads + @warn "affinity mask too small to give OpenBLAS workers their own cores; they will share with Julia threads" mask_ids blas_threads + free_ids = vcat(free_ids, mask_ids) + end + openblas_pinthreads(free_ids[1:blas_threads]) + end + catch err + @warn "thread pinning failed; timings from this run are not publishable" err + end + return nothing +end + opts = parse_cli(ARGS) resultsdir = normpath(joinpath(HERE, "..", "..", "results")) mkpath(resultsdir) BLAS.set_num_threads(opts.blas_threads) +pin_cores(opts.blas_threads) println( "ITensorMPS benchmark harness — mode: ", opts.smoke ? "smoke" : "full", diff --git a/benchmark/profile_sweep.jl b/benchmark/profile_sweep.jl index 1a687de99..4d9c238b8 100644 --- a/benchmark/profile_sweep.jl +++ b/benchmark/profile_sweep.jl @@ -17,6 +17,7 @@ using LinearAlgebra using Profile using Random +using ThreadPinning const HERE = @__DIR__ @@ -29,6 +30,11 @@ chi = parse_flag(ARGS, "chi", 256) N = parse_flag(ARGS, "N", 100) maxdepth = parse_flag(ARGS, "maxdepth", 30) BLAS.set_num_threads(parse_flag(ARGS, "blas-threads", 1)) +try + pinthreads(:affinitymask) # same placement policy as benchmark/run.jl +catch err + @warn "thread pinning failed" err +end include(joinpath(HERE, "suites", "common.jl")) using .BenchCommon diff --git a/benchmark/run.jl b/benchmark/run.jl index c16aa50c8..8ba91ef7a 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -18,6 +18,7 @@ const HERE = @__DIR__ using LinearAlgebra +using ThreadPinning function parse_cli(args) smoke = "--smoke" in args @@ -34,10 +35,40 @@ function parse_cli(args) return (; smoke, suites, blas_threads) end +""" + pin_cores(blas_threads) + +Pin Julia threads 1:1 to distinct cores inside the affinity mask this process was +launched with (`taskset` in `benchmark/slurm/run_all.sbatch`, or Slurm's cpu mask), via +ThreadPinning's `:affinitymask` policy. When the OpenBLAS pool has real workers +(`blas_threads > 1`), pin those to mask cores disjoint from the Julia threads. +No-oversubscription is the launcher's job: every concurrent process gets a disjoint mask. +Pinning failure is a warning, not an error, so smoke runs on unsupported platforms still +work — but cluster timings must never run unpinned, so check the log. +""" +function pin_cores(blas_threads::Int) + try + pinthreads(:affinitymask) + if blas_threads > 1 + mask_ids = findall(==(1), getaffinity(; cutoff = nothing)) .- 1 + free_ids = setdiff(mask_ids, getcpuids()) + if length(free_ids) < blas_threads + @warn "affinity mask too small to give OpenBLAS workers their own cores; they will share with Julia threads" mask_ids blas_threads + free_ids = vcat(free_ids, mask_ids) + end + openblas_pinthreads(free_ids[1:blas_threads]) + end + catch err + @warn "thread pinning failed; timings from this run are not publishable" err + end + return nothing +end + opts = parse_cli(ARGS) resultsdir = joinpath(HERE, "results") mkpath(resultsdir) BLAS.set_num_threads(opts.blas_threads) +pin_cores(opts.blas_threads) println( "MPSKit benchmark harness — mode: ", opts.smoke ? "smoke" : "full", diff --git a/benchmark/slurm/README.md b/benchmark/slurm/README.md index de2c507a7..64359ed22 100644 --- a/benchmark/slurm/README.md +++ b/benchmark/slurm/README.md @@ -1,60 +1,75 @@ -# Rusty submission scripts for the comparative benchmarks +# Rusty submission script for the comparative benchmarks -One sbatch script per (workload x library), all single exclusive node. These are -submitted **manually by the maintainer** — nothing here self-submits. +**One script, one job, one node**: `run_all.sbatch` runs all four suites for both +libraries inside a single exclusive-node allocation, so every compared pair of numbers +comes from the same machine by construction (no cross-job node-type mismatch to check +away afterwards). It is submitted **manually by the maintainer** — nothing here +self-submits. + +## Core discipline (96-core nodes) + +- **Phases 1–2 (DMRG suites 1–2, then TDVP suite 5)**: the MPSKit and ITensorMPS sides + run *concurrently*, each single-threaded (`JULIA_NUM_THREADS=1`, BLAS = 1) and + `taskset`-pinned to one core on *opposite sockets* — 2 of 96 cores active, no shared + core, L3, or memory controller, and the compared pair sees identical machine state. +- **Phase 3 (thread-scaling suite 7)**: one (JULIA_NUM_THREADS × BLAS-threads) grid point + at a time with the node otherwise idle; the largest point uses 16 threads. +- **Inside every Julia process**, ThreadPinning.jl pins Julia threads 1:1 to distinct + cores of the process affinity mask (`pinthreads(:affinitymask)`), and pins OpenBLAS + workers to further disjoint mask cores when the pool is real (`pin_cores` in + `benchmark/run.jl` / the ITensor `run.jl`). Nothing is oversubscribed at any point. + If a log ever shows the "thread pinning failed" warning, the timings from that run are + not publishable. ## Before submitting -1. **Pin one CPU generation** — edit `--constraint=` (and `--partition=` if needed) in - ALL six scripts to the same value. Timings across different node types are not - comparable, and the metadata in the result JSONs records only hostname, not node type. -2. **Instantiate both environments on a login node** (compute-node registry access is not - guaranteed): +1. Edit `--partition` / `--constraint` in `run_all.sbatch` (any single CPU generation). +2. Check the socket→core numbering assumption (`lscpu -e` on the target node): the script + assumes socket 0 = cores 0..47, socket 1 = cores 48..95 and pins the concurrent pair + to cores 0 and 48. Edit `CORE_MPSKIT` / `CORE_ITENSOR` if the numbering interleaves. +3. Instantiate both environments on a login node: ```bash julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' ``` -3. Make sure `julia` is on `PATH` in batch jobs (juliaup or `module load julia` in your - shell init). +4. Optional local rehearsal (same phases at smoke scale, plain bash, no Slurm): -## Submission order (from the repo root) + ```bash + SMOKE=1 bash benchmark/slurm/run_all.sbatch + ``` -Nothing depends on anything else; submit all six whenever. The DMRG pair is the priority: +Then, from the repo root: ```bash -sbatch benchmark/slurm/dmrg_mpskit.sbatch # suites 1-2, chi up to 512 (48 h limit) -sbatch benchmark/slurm/dmrg_itensor.sbatch -sbatch benchmark/slurm/tdvp_mpskit.sbatch # suite 5 (24 h limit) -sbatch benchmark/slurm/tdvp_itensor.sbatch -sbatch benchmark/slurm/threads_mpskit.sbatch # suite 7, 7 grid points sequentially (24 h limit) -sbatch benchmark/slurm/threads_itensor.sbatch +sbatch benchmark/slurm/run_all.sbatch ``` -Walltime calibration: the 2026-07-05 MPSKit pilot (χ = 256, N = 100, two-site, 8 -workstation threads) averaged ≈ 4.5 min/sweep; single-threaded χ = 512 extrapolates to a -few hours per sweep, which sizes the 48 h DMRG limit. χ = 1024 is deliberately out of the -default schedule (days of single-threaded node-time; extend later if the 512 numbers -justify it). +## Walltime calibration + +The 2026-07-05 MPSKit pilot (χ = 256, N = 100, two-site, 8 workstation threads) averaged +≈ 4.5 min/sweep; single-threaded χ = 512 extrapolates to a few hours per sweep. Phases +run back-to-back: DMRG (dominant, tens of hours) → TDVP → thread grid; the 96 h limit is +sized for that with margin. χ = 1024 is deliberately out of the default schedule (days of +single-threaded node-time; extend later if the 512 numbers justify it). If the job times +out, note that each *suite* writes its JSON only at the end — completed χ points inside +an unfinished suite are lost, so bump `--time` rather than resubmitting blindly. -Results accumulate in `benchmark/results/` (one JSON per run/grid point); logs in -`benchmark/slurm/logs/`. If a job times out, the χ points already completed are still -lost (one JSON per *suite*, written at the end), so bump `--time` rather than -resubmitting blindly. +Per-phase, per-library logs land in `benchmark/slurm/logs/` (`--.log`); +results accumulate in `benchmark/results/` (one JSON per suite run / grid point). ## Sanity-gate checklist before reading any timing -Run `julia --project=benchmark benchmark/plot_results.jl` locally on the synced -`results/` directory, then check: +Run `julia --project=benchmark benchmark/plot_results.jl` on the synced `results/` +directory, then check: 1. **Energy gate (suites 1-2)**: at each matched χ, final MPSKit and ITensorMPS energies agree to ~1e-8 (relative). If not, no timing from that pair means anything. 2. **Observable gate (suite 5)**: the ⟨Sz⟩(t) trajectories of the two libraries lie on - top of each other at matched χ (they solve the same Schrödinger equation with the - same integrator family; visible divergence = protocol mismatch, not physics). + top of each other at matched χ (smoke-scale agreement was ~1e-13; visible divergence + = protocol mismatch, not physics). 3. **Thread parity (all)**: `nthreads_julia`, `nthreads_blas`, and `blas_config` match between the two libraries' metadata for every compared pair (suite 7 pairs match - point-by-point). -4. **Same hardware**: all result files report hostnames of the same node type. -5. Only then look at wall times. Report losses as prominently as wins (§4.3). + point-by-point), and no log contains the "thread pinning failed" warning. +4. Only then look at wall times. Report losses as prominently as wins (§4.3). diff --git a/benchmark/slurm/dmrg_itensor.sbatch b/benchmark/slurm/dmrg_itensor.sbatch deleted file mode 100644 index 5b17dc74a..000000000 --- a/benchmark/slurm/dmrg_itensor.sbatch +++ /dev/null @@ -1,20 +0,0 @@ -#!/bin/bash -# Suites 1-2 (finite DMRG, trivial + U(1)), ITensorMPS side. Submit from the repo root: -# sbatch benchmark/slurm/dmrg_itensor.sbatch -# Instantiate the env on the login node first: -# julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' -#SBATCH --job-name=itensor-dmrg -#SBATCH --nodes=1 -#SBATCH --exclusive -#SBATCH --time=48:00:00 -#SBATCH --partition=ccq # EDIT if needed -#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs -#SBATCH --output=benchmark/slurm/logs/%x-%j.out - -set -euo pipefail -cd "${SLURM_SUBMIT_DIR}" -mkdir -p benchmark/slurm/logs - -export JULIA_NUM_THREADS=1 - -julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --suite=1,2 --blas-threads=1 diff --git a/benchmark/slurm/dmrg_mpskit.sbatch b/benchmark/slurm/dmrg_mpskit.sbatch deleted file mode 100644 index 7920ad13d..000000000 --- a/benchmark/slurm/dmrg_mpskit.sbatch +++ /dev/null @@ -1,22 +0,0 @@ -#!/bin/bash -# Suites 1-2 (finite DMRG, trivial + U(1)), MPSKit side. Submit from the repo root: -# sbatch benchmark/slurm/dmrg_mpskit.sbatch -# Instantiate the env on the login node first (compute nodes may lack registry access): -# julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' -#SBATCH --job-name=mpskit-dmrg -#SBATCH --nodes=1 -#SBATCH --exclusive -#SBATCH --time=48:00:00 -#SBATCH --partition=ccq # EDIT if needed -#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs -#SBATCH --output=benchmark/slurm/logs/%x-%j.out - -set -euo pipefail -cd "${SLURM_SUBMIT_DIR}" -mkdir -p benchmark/slurm/logs - -# Single-threaded headline comparison (methodology guardrail §4.3): identical thread -# settings on both sides. Suite 7 is the multithreaded story. -export JULIA_NUM_THREADS=1 - -julia --project=benchmark benchmark/run.jl --suite=1,2 --blas-threads=1 diff --git a/benchmark/slurm/run_all.sbatch b/benchmark/slurm/run_all.sbatch new file mode 100644 index 000000000..a48f12ee9 --- /dev/null +++ b/benchmark/slurm/run_all.sbatch @@ -0,0 +1,102 @@ +#!/bin/bash +# The ENTIRE comparative benchmark (suites 1, 2, 5, 7 x both libraries) as ONE job on ONE +# node, so every published pair of numbers comes from the same machine by construction. +# Submit from the repo root: +# +# sbatch benchmark/slurm/run_all.sbatch +# +# Local end-to-end check (runs the same phases at smoke scale, as plain bash): +# +# SMOKE=1 bash benchmark/slurm/run_all.sbatch +# +# Instantiate both environments on a login node first (compute nodes may lack registry +# access): +# julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' +# julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' +# +# Core discipline (96-core node, nothing else running: --exclusive): +# * Phase 1 (DMRG) and phase 2 (TDVP) run the MPSKit and ITensorMPS sides +# CONCURRENTLY, one single-threaded process per SOCKET (disjoint `taskset` masks, so +# 2 of 96 cores are active and the two sides never share a core, an L3, or a memory +# controller). The compared pair therefore also sees identical machine state. +# * Phase 3 (thread-scaling grid) runs ONE grid point at a time with the node to +# itself; the largest grid point uses 8 + 8 = 16 threads, far below 96. +# * Inside every Julia process, ThreadPinning.jl pins Julia threads 1:1 to distinct +# cores of the process's affinity mask (`pinthreads(:affinitymask)`) and, when the +# OpenBLAS pool has real workers, pins those to further disjoint cores (see +# `pin_cores` in benchmark/run.jl). Nothing is ever oversubscribed. +#SBATCH --job-name=mpskit-bench-all +#SBATCH --nodes=1 +#SBATCH --exclusive +#SBATCH --time=96:00:00 +#SBATCH --partition=ccq # EDIT if needed +#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation (any single one) +#SBATCH --output=benchmark/slurm/logs/%x-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR:-$PWD}" +LOGDIR="benchmark/slurm/logs" +mkdir -p "${LOGDIR}" +TAG="${SLURM_JOB_ID:-local}" + +SMOKE_FLAG="" +if [[ "${SMOKE:-0}" == "1" ]]; then + SMOKE_FLAG="--smoke" + echo ">>> SMOKE mode: pipeline check only, timings meaningless" +fi + +MPSKIT=(julia --project=benchmark benchmark/run.jl) +ITENSOR=(julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl) + +# One core per socket for the concurrent single-threaded pairs. Assumes the common +# contiguous numbering (socket 0 = cores 0..N/2-1, socket 1 = cores N/2..N-1); verify +# with `lscpu -e` on the target node and EDIT if the numbering is interleaved. +NCPU=$(nproc) +CORE_MPSKIT=0 +CORE_ITENSOR=$(( NCPU / 2 )) +echo ">>> $(hostname): ${NCPU} CPU threads; MPSKit on core ${CORE_MPSKIT}, ITensorMPS on core ${CORE_ITENSOR}" + +# Run one MPSKit/ITensorMPS pair concurrently, each single-threaded and pinned to its own +# socket. Fails the job if either side fails. +run_pair() { + local suites="$1" + local phase="$2" + echo ">>> phase ${phase}: suites ${suites}, both libraries concurrently ($(date))" + JULIA_NUM_THREADS=1 taskset -c "${CORE_MPSKIT}" \ + "${MPSKIT[@]}" ${SMOKE_FLAG} --suite="${suites}" --blas-threads=1 \ + > "${LOGDIR}/${TAG}-${phase}-mpskit.log" 2>&1 & + local pid_mpskit=$! + JULIA_NUM_THREADS=1 taskset -c "${CORE_ITENSOR}" \ + "${ITENSOR[@]}" ${SMOKE_FLAG} --suite="${suites}" --blas-threads=1 \ + > "${LOGDIR}/${TAG}-${phase}-itensor.log" 2>&1 & + local pid_itensor=$! + wait "${pid_mpskit}" + wait "${pid_itensor}" + echo ">>> phase ${phase} done ($(date))" +} + +run_pair "1,2" "dmrg" +run_pair "5" "tdvp" + +# Phase 3: thread-scaling grid, one point at a time, node otherwise idle. The affinity +# mask is the whole node; ThreadPinning places the point's threads compactly within it. +GRID=("1 1" "4 1" "1 4" "4 4" "8 1" "1 8" "8 8") +if [[ -n "${SMOKE_FLAG}" ]]; then + GRID=("1 1" "2 2") +fi +for side in mpskit itensor; do + for point in "${GRID[@]}"; do + read -r nj nb <<< "${point}" + echo ">>> phase threads (${side}): julia=${nj} blas=${nb} ($(date))" + if [[ "${side}" == mpskit ]]; then + JULIA_NUM_THREADS="${nj}" "${MPSKIT[@]}" ${SMOKE_FLAG} --suite=7 --blas-threads="${nb}" \ + > "${LOGDIR}/${TAG}-threads-${side}-j${nj}b${nb}.log" 2>&1 + else + JULIA_NUM_THREADS="${nj}" "${ITENSOR[@]}" ${SMOKE_FLAG} --suite=7 --blas-threads="${nb}" \ + > "${LOGDIR}/${TAG}-threads-${side}-j${nj}b${nb}.log" 2>&1 + fi + done +done + +echo ">>> all phases done ($(date)). Results in benchmark/results/, logs in ${LOGDIR}/." +echo ">>> next: julia --project=benchmark benchmark/plot_results.jl" diff --git a/benchmark/slurm/tdvp_itensor.sbatch b/benchmark/slurm/tdvp_itensor.sbatch deleted file mode 100644 index 6ffcc93c0..000000000 --- a/benchmark/slurm/tdvp_itensor.sbatch +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/bash -# Suite 5 (TDVP throughput), ITensorMPS side. Submit from the repo root: -# sbatch benchmark/slurm/tdvp_itensor.sbatch -#SBATCH --job-name=itensor-tdvp -#SBATCH --nodes=1 -#SBATCH --exclusive -#SBATCH --time=24:00:00 -#SBATCH --partition=ccq # EDIT if needed -#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs -#SBATCH --output=benchmark/slurm/logs/%x-%j.out - -set -euo pipefail -cd "${SLURM_SUBMIT_DIR}" -mkdir -p benchmark/slurm/logs - -export JULIA_NUM_THREADS=1 - -julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --suite=5 --blas-threads=1 diff --git a/benchmark/slurm/tdvp_mpskit.sbatch b/benchmark/slurm/tdvp_mpskit.sbatch deleted file mode 100644 index 1db5f5fe3..000000000 --- a/benchmark/slurm/tdvp_mpskit.sbatch +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/bash -# Suite 5 (TDVP throughput), MPSKit side. Submit from the repo root: -# sbatch benchmark/slurm/tdvp_mpskit.sbatch -#SBATCH --job-name=mpskit-tdvp -#SBATCH --nodes=1 -#SBATCH --exclusive -#SBATCH --time=24:00:00 -#SBATCH --partition=ccq # EDIT if needed -#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs -#SBATCH --output=benchmark/slurm/logs/%x-%j.out - -set -euo pipefail -cd "${SLURM_SUBMIT_DIR}" -mkdir -p benchmark/slurm/logs - -export JULIA_NUM_THREADS=1 - -julia --project=benchmark benchmark/run.jl --suite=5 --blas-threads=1 diff --git a/benchmark/slurm/threads_itensor.sbatch b/benchmark/slurm/threads_itensor.sbatch deleted file mode 100644 index 1496db986..000000000 --- a/benchmark/slurm/threads_itensor.sbatch +++ /dev/null @@ -1,24 +0,0 @@ -#!/bin/bash -# Suite 7 (thread scaling), ITensorMPS side. Same structure and grid as -# threads_mpskit.sbatch. Submit from the repo root: -# sbatch benchmark/slurm/threads_itensor.sbatch -#SBATCH --job-name=itensor-threads -#SBATCH --nodes=1 -#SBATCH --exclusive -#SBATCH --time=24:00:00 -#SBATCH --partition=ccq # EDIT if needed -#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs -#SBATCH --output=benchmark/slurm/logs/%x-%j.out - -set -euo pipefail -cd "${SLURM_SUBMIT_DIR}" -mkdir -p benchmark/slurm/logs - -# Keep the grid identical to threads_mpskit.sbatch. -GRID=("1 1" "4 1" "1 4" "4 4" "8 1" "1 8" "8 8") - -for point in "${GRID[@]}"; do - read -r nj nb <<< "${point}" - echo "=== suite 7 (ITensorMPS): julia=${nj} blas=${nb} ===" - JULIA_NUM_THREADS="${nj}" julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl --suite=7 --blas-threads="${nb}" -done diff --git a/benchmark/slurm/threads_mpskit.sbatch b/benchmark/slurm/threads_mpskit.sbatch deleted file mode 100644 index e0590e16a..000000000 --- a/benchmark/slurm/threads_mpskit.sbatch +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash -# Suite 7 (thread scaling), MPSKit side. One exclusive node; one julia process per -# (JULIA_NUM_THREADS x BLAS-threads) grid point, run sequentially so every point sees an -# otherwise-idle node. Submit from the repo root: -# sbatch benchmark/slurm/threads_mpskit.sbatch -#SBATCH --job-name=mpskit-threads -#SBATCH --nodes=1 -#SBATCH --exclusive -#SBATCH --time=24:00:00 -#SBATCH --partition=ccq # EDIT if needed -#SBATCH --constraint=icelake # EDIT: pin ONE CPU generation, same for ALL jobs -#SBATCH --output=benchmark/slurm/logs/%x-%j.out - -set -euo pipefail -cd "${SLURM_SUBMIT_DIR}" -mkdir -p benchmark/slurm/logs - -# (julia threads, blas threads) grid — JULIA_NUM_THREADS is launch-time-only, hence the -# one-process-per-point loop. Keep the grid identical to threads_itensor.sbatch. -GRID=("1 1" "4 1" "1 4" "4 4" "8 1" "1 8" "8 8") - -for point in "${GRID[@]}"; do - read -r nj nb <<< "${point}" - echo "=== suite 7 (MPSKit): julia=${nj} blas=${nb} ===" - JULIA_NUM_THREADS="${nj}" julia --project=benchmark benchmark/run.jl --suite=7 --blas-threads="${nb}" -done From 5338eb7ec102a89caba01e9c842e544f1880a17e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 6 Jul 2026 10:20:19 -0400 Subject: [PATCH 6/8] bench: fix cluster julia resolution in run_all.sbatch MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two failures from the first Rusty submissions: - job 6579930: batch shells have no interactive PATH, so 'taskset ... julia' died with 'No such file or directory'; - job 6579944: 'module load julia' is 1.11.2, which silently ignores the repo's Pkg [workspace] section, so --project=benchmark cannot see the root manifest and errors with a misleading 'Package ThreadPinning ... does not seem to be installed'. The script now resolves julia itself: $JULIA override, else the newest 1.12+ binary in the shared-home juliaup depot (~/.julia/juliaup, GPFS, identical on all nodes and the same install that resolved the manifests), else PATH — then prints the choice and fail-fasts on julia < 1.12. Also adds an ERR trap pointing at the per-process logs. Verified by rehearsing the full driver under env -i with no julia on PATH (the cluster batch-shell condition): all phases pass. Keeps the user's partition=gen edit; drops the module load. Co-Authored-By: Claude Fable 5 --- benchmark/slurm/README.md | 21 +++++++++++++++++---- benchmark/slurm/run_all.sbatch | 27 +++++++++++++++++++++++---- 2 files changed, 40 insertions(+), 8 deletions(-) diff --git a/benchmark/slurm/README.md b/benchmark/slurm/README.md index 64359ed22..91a7d6f79 100644 --- a/benchmark/slurm/README.md +++ b/benchmark/slurm/README.md @@ -25,16 +25,29 @@ self-submits. 1. Edit `--partition` / `--constraint` in `run_all.sbatch` (any single CPU generation). 2. Check the socket→core numbering assumption (`lscpu -e` on the target node): the script - assumes socket 0 = cores 0..47, socket 1 = cores 48..95 and pins the concurrent pair - to cores 0 and 48. Edit `CORE_MPSKIT` / `CORE_ITENSOR` if the numbering interleaves. -3. Instantiate both environments on a login node: + assumes contiguous numbering (socket 0 = first half, socket 1 = second half) and pins + the concurrent pair to cores 0 and `nproc/2`. Edit `CORE_MPSKIT` / `CORE_ITENSOR` if + the numbering interleaves. +3. **Julia ≥ 1.12 is required** (`benchmark/` is a Pkg `[workspace]` member; 1.11 and + older silently ignore `[workspace]` and then claim the packages are "not installed"). + Do **not** rely on `module load julia` (currently 1.11.2) or on the batch shell's + PATH (it lacks your interactive PATH entirely). The script resolves julia by itself: + `$JULIA` env var if set, else the newest 1.12+ binary in the shared-home juliaup + depot (`~/.julia/juliaup/julia-1.1x*/bin/julia`, visible from all nodes), else PATH — + and fail-fasts with a version check. To pin explicitly: + + ```bash + JULIA=~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/bin/julia sbatch --export=ALL benchmark/slurm/run_all.sbatch + ``` + +4. Instantiate both environments on a login node **with that same julia**: ```bash julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' ``` -4. Optional local rehearsal (same phases at smoke scale, plain bash, no Slurm): +5. Optional local rehearsal (same phases at smoke scale, plain bash, no Slurm): ```bash SMOKE=1 bash benchmark/slurm/run_all.sbatch diff --git a/benchmark/slurm/run_all.sbatch b/benchmark/slurm/run_all.sbatch index a48f12ee9..581debf38 100644 --- a/benchmark/slurm/run_all.sbatch +++ b/benchmark/slurm/run_all.sbatch @@ -14,7 +14,8 @@ # julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' # julia --project=benchmark/comparisons/itensor -e 'using Pkg; Pkg.instantiate()' # -# Core discipline (96-core node, nothing else running: --exclusive): +# Core discipline (node sizes vary, e.g. 64-96 cores — the script adapts via nproc; +# nothing else running: --exclusive): # * Phase 1 (DMRG) and phase 2 (TDVP) run the MPSKit and ITensorMPS sides # CONCURRENTLY, one single-threaded process per SOCKET (disjoint `taskset` masks, so # 2 of 96 cores are active and the two sides never share a core, an L3, or a memory @@ -29,7 +30,7 @@ #SBATCH --nodes=1 #SBATCH --exclusive #SBATCH --time=96:00:00 -#SBATCH --partition=ccq # EDIT if needed +#SBATCH --partition=gen #SBATCH --constraint=icelake # EDIT: pin ONE CPU generation (any single one) #SBATCH --output=benchmark/slurm/logs/%x-%j.out @@ -38,6 +39,24 @@ cd "${SLURM_SUBMIT_DIR:-$PWD}" LOGDIR="benchmark/slurm/logs" mkdir -p "${LOGDIR}" TAG="${SLURM_JOB_ID:-local}" +trap 'echo ">>> FAILED — check ${LOGDIR}/${TAG}-*.log for the per-process error"' ERR + +# --- Julia resolution ----------------------------------------------------------------- +# Julia >= 1.12 is REQUIRED: benchmark/ is a Pkg [workspace] member of the repo root, and +# older julia ignores [workspace] entirely — the symptom is a bogus +# "Package X ... does not seem to be installed: Run Pkg.instantiate()" (seen with +# `module load julia` = 1.11.2, job 6579944). Do NOT `module load julia` here; batch +# shells also lack the interactive PATH (job 6579930: "taskset: failed to execute +# julia"). Instead: honor an explicit $JULIA, else pick the newest 1.12+ binary in the +# shared-home juliaup depot (~/.julia/juliaup — GPFS, so identical on workstations and +# cluster nodes, and the same install that resolved the manifests), else fall back to +# PATH and let the version gate below catch it. +if [[ -z "${JULIA:-}" ]]; then + JULIA=$(ls -d "${HOME}"/.julia/juliaup/julia-1.1[2-9]*/bin/julia 2>/dev/null | sort -V | tail -1 || true) + [[ -n "${JULIA}" && -x "${JULIA}" ]] || JULIA=julia +fi +echo ">>> julia: $("${JULIA}" --version) ($(command -v "${JULIA}" || echo "${JULIA}"))" +"${JULIA}" -e 'VERSION >= v"1.12" || (println(stderr, "ERROR: julia >= 1.12 required (benchmark/ is a Pkg workspace member; older julia cannot resolve it)"); exit(1))' SMOKE_FLAG="" if [[ "${SMOKE:-0}" == "1" ]]; then @@ -45,8 +64,8 @@ if [[ "${SMOKE:-0}" == "1" ]]; then echo ">>> SMOKE mode: pipeline check only, timings meaningless" fi -MPSKIT=(julia --project=benchmark benchmark/run.jl) -ITENSOR=(julia --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl) +MPSKIT=("${JULIA}" --project=benchmark benchmark/run.jl) +ITENSOR=("${JULIA}" --project=benchmark/comparisons/itensor benchmark/comparisons/itensor/run.jl) # One core per socket for the concurrent single-threaded pairs. Assumes the common # contiguous numbering (socket 0 = cores 0..N/2-1, socket 1 = cores N/2..N-1); verify From 4df5fa0071e52bc9b20465972d6b4816e8443fa8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 6 Jul 2026 21:04:21 -0400 Subject: [PATCH 7/8] bench: latest_result picks newest file by mtime, not filename 'suite1_..._N100_.json' sorts lexically before 'suite1_..._N20_.json', so the first cluster analysis silently compared stale N=20 smoke files. Co-Authored-By: Claude Fable 5 --- benchmark/plot_results.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/benchmark/plot_results.jl b/benchmark/plot_results.jl index a5147af94..157f2aae0 100644 --- a/benchmark/plot_results.jl +++ b/benchmark/plot_results.jl @@ -206,13 +206,15 @@ plot_all_suite_results(resultdir) latest_result(resultsdir, prefix) -> Union{String, Nothing} Path of the most recent result file starting with `prefix` (exact prefix match, so -`"suite1"` does not match `"itensor_suite1"`), or `nothing`. +`"suite1"` does not match `"itensor_suite1"`), or `nothing`. "Most recent" is by file +mtime, NOT by name: filenames embed N before the timestamp, and lexically "N100..." sorts +before "N20...", which once made this function silently compare stale smoke files. """ function latest_result(resultsdir::AbstractString, prefix::AbstractString) - files = filter(readdir(resultsdir)) do fname - startswith(fname, prefix) && endswith(fname, ".json") + files = filter(readdir(resultsdir; join = true)) do path + startswith(basename(path), prefix) && endswith(path, ".json") end - return isempty(files) ? nothing : joinpath(resultsdir, sort(files)[end]) + return isempty(files) ? nothing : argmax(mtime, files) end """ From 6ce45ddc37790f0f1f1d090f15550ad43db36df4 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 7 Jul 2026 13:32:03 -0400 Subject: [PATCH 8/8] bench: round-1 ANALYSIS from Rusty job 6581850 + pin_cores readback fix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit All sanity gates pass (energies <= 4.6e-9 relative at matched chi; TDVP Sz trajectories agree to ~1.5e-8). Headline findings: - TDVP (matched inner Krylov tol): MPSKit 2-4x faster at all chi. - Two-site DMRG: ITensor 1.5-1.6x faster dense, ~12x faster per U(1) sweep at chi=512 — attributed (REVIEW-flagged) to inner-eigensolver policy: ITensor caps local solves at krylovdim=3/maxiter=1 while MPSKit's dynamic tolerances drive krylovdim-30 solves toward 1e-14 under the run-past-convergence protocol; per-sweep costs and allocation volumes match that story. - U(1) grow-from-seed reaches the 1e-8 gate in ~2 s at any chi cap (workflow win, labelled as such). - Thread grid: MPSKit scales to 4.8x at (8,8) vs ITensor's 2.2x, from an 8x slower single-core baseline on this workload. pin_cores fix: read the affinity mask BEFORE pinthreads narrows the calling thread (JULIA_NUM_THREADS=1 previously saw a 1-core mask and bailed out of BLAS pinning — grid points (1,4)/(1,8) ran with floating BLAS workers), and re-pin Julia threads after openblas_pinthreads, whose pool includes the calling thread at nthreads=1. Verified locally for both regimes. Co-Authored-By: Claude Fable 5 --- benchmark/ANALYSIS.md | 198 +++++++++++++++++++-------- benchmark/comparisons/itensor/run.jl | 27 +++- benchmark/run.jl | 27 +++- 3 files changed, 183 insertions(+), 69 deletions(-) diff --git a/benchmark/ANALYSIS.md b/benchmark/ANALYSIS.md index ed1d36c34..b132420b7 100644 --- a/benchmark/ANALYSIS.md +++ b/benchmark/ANALYSIS.md @@ -1,62 +1,146 @@ # MPSKit vs ITensorMPS — investigation notes -Status: **skeleton — no full-run results yet.** This file is filled in only from real -Rusty result JSONs in `results/` (never from smoke runs, never from memory). Every claim -must cite the result file(s) it comes from. - -## Questions this round must answer - -### Q1. Sanity gates (prerequisite for everything below) - -- [ ] Suites 1-2: do MPSKit and ITensorMPS converge to the same energy (~1e-8 relative) - at every matched χ? (`compare_dmrg_suite` table in `plot_results.jl` output.) -- [ ] Suite 5: do the ⟨Sz⟩(t) trajectories coincide at matched χ? -- [ ] Metadata parity: same node type, same `nthreads_julia` / `nthreads_blas` / - `blas_config` on both sides of every compared pair? - -### Q2. Headline: time-to-accuracy (suites 1-2) - -- [ ] At each χ, which library reaches E_ref + 1e-8·|E_ref| first, and by what factor? -- [ ] Does the ratio change with χ (i.e. is the crossover at small or large χ)? -- [ ] Does the U(1) suite change the picture relative to the trivial suite (block-sparse - overhead vs dense throughput)? -- [ ] Convergence-quality check: does either side plateau above the gate at some χ - (local minimum), and does that correlate with the fixed-χ no-noise protocol? - -### Q3. TDVP throughput (suite 5) - -- [ ] Seconds of wall time per unit physical time vs χ, per library — who wins, and does - the scaling exponent in χ differ (should be ~χ³; deviations are implementation - overhead)? -- [ ] How large is the growth phase (two-site) cost relative to the measure phase? - -### Q4. Thread scaling (suite 7) - -- [ ] Where does each library peak on the (julia-threads × BLAS-threads) grid? -- [ ] MPSKit: how much does pure Julia-threading ((8,1) vs (1,1)) buy on this dense - workload? ITensor: same question for BLAS ((1,8) vs (1,1)). -- [ ] Do mixed settings ((4,4), (8,8)) oversubscribe and regress? - -### Q5. Where does the time go? (diagnostics) - -- [ ] GC fraction and allocation volume per sweep, per library, vs χ (`gctimes` / - `allocd_bytes` fields). Is either side GC-bound at any χ? -- [ ] `profile_sweep.jl` on the interesting χ points (pick after Q2): what fraction of a - sweep is BLAS / SVD / eigsolve vs MPSKit's own environment and permutation code? -- [ ] Odd/even sweep-cost alternation: the χ = 256 MPSKit pilot (2026-07-05, workstation, - 8 threads) showed alternating expensive/cheap sweeps (≈450 s vs ≈50 s). Convergence - of the local Krylov solves after the state has converged is the likely cause - (cheap sweeps do almost no eigsolve work); confirm with `gctimes`/profile on the - cluster runs and check whether ITensor shows the same pattern. - -## Findings - -_(empty — populate per question above, with result-file citations, after the Rusty runs.)_ +Status: **round-1 findings recorded** (2026-07-07), from Rusty job 6581850: one exclusive +node (`worker6192`, 64 CPU threads, 2 sockets), all four suites, both libraries, julia +1.12.6, OpenBLAS (LBT ILP64), N = 100 spin-1 Heisenberg throughout. Result files cited +inline; raw JSONs in `results/`, per-process logs in `slurm/logs/6581850-*`. + +## Q1. Sanity gates — ALL PASS + +- **Energy gate (suites 1–2)**: at every matched χ ∈ {64, 128, 256, 512}, final energies + agree to ≤ 6.4e-7 absolute on |E| ≈ 138.94, i.e. relative ≤ 4.6e-9 < 1e-8. (Tables + printed by `plot_results.jl`; files `suite{1,2}_*_N100_20260706-*.json` vs + `itensor_suite{1,2}_*_N100_20260706-*.json`.) +- **Observable gate (suite 5)**: mid-chain ⟨Sz⟩(t) agrees to ~1.5e-8 at all χ after + t = 1.3–1.4 of evolution (`suite5_tdvp_N100_20260707-015719.json` vs + `itensor_suite5_tdvp_N100_20260707-093904.json`); trajectories are visually identical + (`results/compare_suite5_tdvp.png`). +- **Thread/BLAS parity**: same node, same julia, `nthreads_julia`/`nthreads_blas` match + per compared pair; suites 1/2/5 ran the pair concurrently pinned to opposite sockets. + **Caveat**: the four grid points (1,4)/(1,8) × both libraries ran with *unpinned* BLAS + workers (floating on the otherwise-idle node) due to a pin_cores readback bug, fixed + after the run (see Q4). Symmetric across libraries; placement noise only. + +## Q2. Time-to-accuracy, finite two-site DMRG (suites 1–2) + +Time to E − E_ref ≤ 1e-8·|E_ref| (E_ref = best energy across both libraries), seconds: + +| χ | MPSKit (trivial) | ITensor (trivial) | MPSKit (U(1)) | ITensor (U(1)) | +|---|---|---|---|---| +| 64 | 17.2 | 16.9 | 1.5 | 11.2 | +| 128 | 122 | 100 | 1.7 | 23.7 | +| 256 | 1017 | 653 | 2.3 | 80.4 | +| 512 | 5660 | 3685 | 2.3 | 301 | + +- **Trivial suite (matched random full-χ starts — the apples-to-apples comparison): + ITensor is ~1.5–1.6× faster at χ ≥ 128**; parity at χ = 64. +- **U(1) suite: MPSKit's numbers are protocol-dominated, not throughput.** Its + grow-from-seed start crosses the gate during the cheap early growth sweeps (this gapped + model already meets 1e-8-relative at effective χ ≈ 64), giving 2-second + "time-to-accuracy" at every χ cap. Legitimate as an idiomatic-workflow comparison — + each library runs its natural protocol, and practitioners do grow — but the per-sweep + numbers below are the implementation comparison. +- Steady-state per-sweep cost at full χ = 512 (last sweeps): MPSKit U(1) ≈ 570–590 s vs + ITensor ≈ 49 s — **ITensor ~12× faster per U(1) two-site sweep**; trivial: MPSKit + erratic 308–7178 s (!) vs ITensor flat ≈ 608 s. + +### Why: inner-eigensolver work policy, not kernels (root-cause hypothesis) + +Verified from both sources: + +- ITensorMPS `dmrg` caps the local solve at `eigsolve_krylovdim = 3`, + `eigsolve_maxiter = 1` (src/dmrg.jl:168–170) — a fixed, tiny Krylov budget per site + update; convergence comes from many cheap sweeps. +- MPSKit's `Defaults.alg_eigsolve` uses `krylovdim = 30` with **dynamic tolerances** + (`DynamicTol`, `tol_min = 1e-14`, `tol_factor = 1e-3`; src/utility/defaults.jl:26–61): + the inner tolerance tracks the outer error × 1e-3, floored at 1e-14. Under this + suite's run-past-convergence protocol (`tol = 0`), the inner solves are driven toward + machine precision with a 30-vector Krylov space — enormous and *variable* per-sweep + work. This matches the observed erratic sweep costs (308 → 7178 → 4480 s at χ = 512) + and allocation swings (120 → 2195 GiB/sweep) where ITensor is flat (608 s, 456 GiB). + +Corroboration: in TDVP (Q3), where the inner Krylov tolerance is **matched** (1e-10 both +sides), MPSKit is uniformly 2–4× *faster* — the tensor-contraction stack is not the +bottleneck; the DMRG inner-solve policy is. + + +## Q3. TDVP throughput (suite 5) — MPSKit wins 2–4× + +Wall-seconds per unit physical time, single-site measure phase (matched Krylov tol +1e-10, ComplexF64 both sides): + +| χ | MPSKit | ITensorMPS | ratio | +|---|---|---|---| +| 64 | 217 | 890 | 4.1× | +| 128 | 1237 | 3972 | 3.2× | +| 256 | 8490 | 18924 | 2.2× | + +The two-site growth phase shows the same ordering (χ = 256: 6758 s vs 17803 s, 2.6×). +The ratio shrinks with χ as O(χ³) BLAS (identical on both sides) increasingly dominates. + +## Q4. Thread scaling (suite 7, χ = 256, 3 sweeps from random full-χ start) + +Total workload walltime (s) and speedup vs the library's own (1,1): + +| (julia, blas) | MPSKit | speedup | ITensorMPS | speedup | +|---|---|---|---|---| +| (1,1) | 2670 | 1.00 | 335 | 1.00 | +| (4,1) | 2144 | 1.25 | 296 | 1.13 | +| (8,1) | 1837 | 1.45 | 278 | 1.21 | +| (1,4) | 1230 | 2.17 | 231 | 1.45 | +| (1,8) | 863 | 3.10 | 213 | 1.57 | +| (4,4) | 804 | 3.32 | 170 | 1.97 | +| (8,8) | 554 | 4.82 | 151 | 2.21 | + +- MPSKit scales notably better (4.8× at (8,8) vs 2.2×) but from an 8× worse + single-thread baseline on this workload (which Q2's inner-solver policy explains — the + extra Krylov work is BLAS-heavy and parallelizes). Even at (8,8) MPSKit (554 s) does + not catch ITensor's single-core 335 s here. +- For both libraries BLAS threads beat Julia threads on this dense two-site workload; + mixed (8,8) is best and does not oversubscribe. +- (1,4)/(1,8) points: BLAS workers were unpinned (floating on the idle node) due to the + `pin_cores` bug — it read the affinity mask *after* `pinthreads` had narrowed the + calling thread to one core, so BLAS pinning bailed out (only manifest at + `JULIA_NUM_THREADS = 1`; with > 1 the toplevel task sits on an unpinned thread and read + the full mask). Fixed post-run: mask read before pinning + Julia threads re-pinned + after `openblas_pinthreads` (whose pool includes the calling thread when + `JULIA_NUM_THREADS = 1` and would otherwise relocate it). Symmetric noise; rerun of + those four points is optional (~30 min total on one node). + +## Q5. Where the time goes (diagnostics) + +- GC fractions at χ = 512 (steady state): MPSKit 7% (trivial) / 18% (U(1)); ITensor 12% + (trivial) / 31% (U(1)). ITensor pays a *higher* GC share yet is much faster per sweep — + more evidence the gap is solver work, not memory management. +- Allocation volume tracks the Krylov work: MPSKit's expensive sweeps allocate ~2.2 + TiB/sweep (trivial χ = 512) vs ITensor's flat 456 GiB. +- `profile_sweep.jl` (workstation, χ = 32 shakedown) shows the sweep dominated by + KrylovKit `eigsolve`/Lanczos over `tensorcontract!` — consistent with the policy story; + a cluster-node profile at χ = 256–512 is queued as follow-up. + +## Round-1 conclusions + +1. Both libraries solve the same problems to the same accuracy (all gates pass); the + harness, parity choices, and pinning discipline held up on the cluster. +2. **MPSKit's tensor stack is fast** — 2–4× ahead in TDVP at matched inner tolerances. +3. **MPSKit's DMRG defaults spend far more inner-eigensolver work per sweep** than + ITensor's (krylovdim 30 + dynamic tol → machine precision, vs krylovdim 3 × 1 sweep), + which costs it the two-site DMRG comparison (1.5–1.6× dense, ~12× per-sweep U(1)) + under this fixed-sweep protocol. +4. The U(1) grow-from-seed workflow is a massive practical win (gate in ~2 s at any χ + cap) and worth advertising as such, clearly labelled as a workflow comparison. ## Follow-ups queued for later rounds -- χ-ramp warm-start protocol (both libraries) instead of fresh random starts per χ. -- SU(2) capability suite (MPSKit-only, vs its own U(1) run). -- TeNPy comparisons; quasi-2D cylinder workload. -- Publishing: populate `docs/src/benchmarks.md` + README chart only after the sanity - gates pass and the maintainer has reviewed the parity tables. +- **Decisive experiment**: rerun suite 1 (and the U(1) per-sweep measurement) with + MPSKit `DMRG2(alg_eigsolve = ...)` matched to ITensor's inner budget (krylovdim = 3, + maxiter = 1, fixed tol) — if the gap closes, the round-1 story is confirmed and a + defaults discussion for MPSKit's ground-state solvers is warranted. +- Rerun grid points (1,4)/(1,8) with the fixed pin_cores (optional, ~30 min). +- Cluster-node `profile_sweep.jl` at χ = 256/512 to quantify eigsolve vs contraction vs + SVD shares directly. +- χ-ramp warm-start protocol; SU(2) capability suite; TeNPy; quasi-2D cylinder. +- Publishing: populate `docs/src/benchmarks.md` + README chart only after the losses/wins + above are reviewed and the follow-up experiments settle the DMRG story. diff --git a/benchmark/comparisons/itensor/run.jl b/benchmark/comparisons/itensor/run.jl index 4e925a6dd..a747c9f8c 100644 --- a/benchmark/comparisons/itensor/run.jl +++ b/benchmark/comparisons/itensor/run.jl @@ -50,15 +50,30 @@ work — but cluster timings must never run unpinned, so check the log. """ function pin_cores(blas_threads::Int) try + # Read the launch mask BEFORE pinning: pinthreads narrows the calling thread's + # affinity, and with JULIA_NUM_THREADS=1 a later getaffinity() would see only + # that single core (bug found in cluster job 6581850, grid points (1,4)/(1,8): + # the readback returned [0] and BLAS pinning crashed out). + mask_ids = findall(==(1), getaffinity(; cutoff = nothing)) .- 1 pinthreads(:affinitymask) if blas_threads > 1 - mask_ids = findall(==(1), getaffinity(; cutoff = nothing)) .- 1 - free_ids = setdiff(mask_ids, getcpuids()) - if length(free_ids) < blas_threads - @warn "affinity mask too small to give OpenBLAS workers their own cores; they will share with Julia threads" mask_ids blas_threads - free_ids = vcat(free_ids, mask_ids) + julia_ids = getcpuids() + free_ids = setdiff(mask_ids, julia_ids) + if length(free_ids) >= blas_threads + openblas_pinthreads(free_ids[1:blas_threads]) + else + # Not enough mask cores to give BLAS workers their own: leave them + # unpinned. They still inherit the process affinity mask, so they float + # within it (no oversubscription of anything outside the mask) — noted + # as placement noise, not a disqualifying failure. + @warn "mask has too few free cores to pin OpenBLAS workers; they float within the mask" mask_ids blas_threads end - openblas_pinthreads(free_ids[1:blas_threads]) + # When Julia runs single-threaded, OpenBLAS counts the CALLING thread as one + # of its pool slots, so openblas_pinthreads just moved the main Julia thread + # onto a BLAS core. Re-pin the Julia threads to the placement chosen above + # (explicit-cpuid variant; a second :affinitymask call would re-read the + # already-narrowed mask). + pinthreads(julia_ids) end catch err @warn "thread pinning failed; timings from this run are not publishable" err diff --git a/benchmark/run.jl b/benchmark/run.jl index 8ba91ef7a..889156f91 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -48,15 +48,30 @@ work — but cluster timings must never run unpinned, so check the log. """ function pin_cores(blas_threads::Int) try + # Read the launch mask BEFORE pinning: pinthreads narrows the calling thread's + # affinity, and with JULIA_NUM_THREADS=1 a later getaffinity() would see only + # that single core (bug found in cluster job 6581850, grid points (1,4)/(1,8): + # the readback returned [0] and BLAS pinning crashed out). + mask_ids = findall(==(1), getaffinity(; cutoff = nothing)) .- 1 pinthreads(:affinitymask) if blas_threads > 1 - mask_ids = findall(==(1), getaffinity(; cutoff = nothing)) .- 1 - free_ids = setdiff(mask_ids, getcpuids()) - if length(free_ids) < blas_threads - @warn "affinity mask too small to give OpenBLAS workers their own cores; they will share with Julia threads" mask_ids blas_threads - free_ids = vcat(free_ids, mask_ids) + julia_ids = getcpuids() + free_ids = setdiff(mask_ids, julia_ids) + if length(free_ids) >= blas_threads + openblas_pinthreads(free_ids[1:blas_threads]) + else + # Not enough mask cores to give BLAS workers their own: leave them + # unpinned. They still inherit the process affinity mask, so they float + # within it (no oversubscription of anything outside the mask) — noted + # as placement noise, not a disqualifying failure. + @warn "mask has too few free cores to pin OpenBLAS workers; they float within the mask" mask_ids blas_threads end - openblas_pinthreads(free_ids[1:blas_threads]) + # When Julia runs single-threaded, OpenBLAS counts the CALLING thread as one + # of its pool slots, so openblas_pinthreads just moved the main Julia thread + # onto a BLAS core. Re-pin the Julia threads to the placement chosen above + # (explicit-cpuid variant; a second :affinitymask call would re-read the + # already-narrowed mask). + pinthreads(julia_ids) end catch err @warn "thread pinning failed; timings from this run are not publishable" err