Reuse exponential! workspaces across phiv!/expv! calls by ChrisRackauckas-Claude · Pull Request #255 · SciML/ExponentialUtilities.jl

ChrisRackauckas-Claude · 2026-07-17T09:30:46Z

Summary

Stacked on #254 (includes its commit; review only the last commit here until #254 merges).

phiv_dense! (called by every Krylov phiv!) and the non-Hermitian branch of expv! called exponential! without a preallocated workspace, so every call re-allocated the ExpMethodHigham2005Base scratch matrices plus a LinearSolve Padé workspace for the (m+p)×(m+p) Hessenberg exponential. After #254, this was the dominant remaining per-step allocation in the Krylov-based exponential integrators.

PhivCache and ExpvCache now lazily hold alloc_mem workspaces — one per distinct extended-matrix size and exponential method requested, because a single integrator step evaluates phi functions of several orders (ETDRK4 alternates k=1/k=3; phiv_timestep! inside the EPIRK methods uses different column counts per stage), so a single-slot cache would be invalidated on almost every call. The store is bounded (reset past 16 entries) in case adaptive Krylov wanders over many subspace sizes. phiv! gains an expmethod keyword forwarded to phiv_dense! together with the workspace; matrix types without an alloc_mem preallocation (e.g. GPU) keep the previous two-argument exponential! path.

Measurements (local, Julia 1.12.6)

Per call (n=400, m=30, k=3, preallocated caches): phiv! 90,128 → 1,152 bytes; expv! → 752 bytes.

Per step on the Kuramoto–Sivashinsky pseudospectral problem (N=128, non-allocating nonlinearity), combined with #254 (baseline = released v1.33.0):

method	released	after
ETDRK2(krylov=true, m=20)	114 KB	3.0 KB
ETDRK4(krylov=true, m=20)	405 KB	12 KB
EPIRK4s3A	636 KB	11 KB
EPIRK5P2	981 KB	16 KB
Exp4	1180 KB	23 KB
Exprb32	250 KB	3.7 KB
Exprb43	506 KB	8.4 KB

Caching-mode (non-Krylov) ETDRK methods were already 0 bytes/step and are unchanged.

Tests

New testset: alternating phi orders against one cache reuse exactly one workspace per size and match the uncached results; expv! workspace identity is stable across calls.
Ran Pkg.test() locally: Core/basictests.jl 687 pass, 1 broken (pre-existing @test_skip).

Note: please ignore this PR until reviewed by @ChrisRackauckas.

🤖 Generated with Claude Code

https://claude.ai/code/session_01EDVGmmovzD5Aos3fFoPwMY

ChrisRackauckas-Claude · 2026-07-22T21:24:30Z

Both addressed in the latest commit.

Any → concrete. You were right — my "must be heterogeneous" reasoning was wrong. The workspace type is size-independent (typeof(alloc_mem(30×30)) === typeof(alloc_mem(400×400)), since LinearSolve's DefaultLinearSolver is a size-independent wrapper); it depends only on the element type and the method. So expcache is now Vector{Tuple{Int, W}} with W the concrete workspace type, recovered at construction via Base.promote_op (no allocation; Nothing for the fallback). get_expcache! now dispatches on the method type, so its return is concretely typed (W or Nothing, never a union).

Why there were still allocations. They were not the cache — they were entirely inside ExpMethodHigham2005Base itself (confirmed by comparison: the non-Base ExpMethodHigham2005 was already 0 alloc/call). Two sources: LAPACK.gebal! allocated a balancing vector per call, and the Padé coefficient arrays C = T[...] were allocated fresh each call. Fixed by using gebal_noalloc! with a preallocated scale buffer, and const coefficient tuples fed through a type-stable function barrier.

Per-call (warm cache, Julia 1.12), released → now:

	released	now
`phiv!`	~90 KB	0 B
`expv!`	~1 KB	96 B
`phiv_timestep!`	~3 KB	352 B

phiv! is now fully non-allocating. exponential!(_, ExpMethodHigham2005Base, work) itself goes 464 → 0 B/call, correct to machine precision across all norm ranges and Float64/ComplexF64/Float32. The remaining expv! 96 B (a column-view that does not elide) and phiv_timestep! 352 B (its own [t]/coeffs timestepping temporaries) are small, size-independent, and separate from both the cache and the exp method — I can chase them in a follow-up if you want them at zero too.

Local: Core 693 pass / QA 20 pass (1 pre-existing broken each).

🤖 Generated with Claude Code

phiv_dense! (called by every Krylov phiv!) and the non-Hermitian branch of expv! previously called exponential! without a preallocated workspace, so every call re-allocated the Higham2005Base scratch matrices and a LinearSolve Pade workspace for the (m+p)x(m+p) Hessenberg exponential. For the Krylov-based exponential integrators this was the dominant per-step allocation after the PhivCache resize fix. PhivCache and ExpvCache now lazily hold alloc_mem workspaces, one per distinct extended-matrix size and exponential method requested (a single integrator step evaluates phi functions of several orders, so a single-slot cache would be invalidated on every call). The store is bounded and reset if adaptive Krylov wanders over many sizes. phiv! also gains an expmethod keyword, forwarded to phiv_dense! together with the workspace; matrix types without an alloc_mem preallocation (e.g. GPU) fall back to the previous two-argument exponential! path. Measured per call (n=400, m=30, k=3, preallocated caches, Julia 1.12.6): phiv! 90128 -> 1152 bytes; expv! -> 752 bytes. Per step on the Kuramoto-Sivashinsky pseudospectral problem (N=128, non-allocating nonlinearity), together with the PhivCache resize fix: - ETDRK2(krylov=true, m=20): 114 KB -> 3.0 KB - ETDRK4(krylov=true, m=20): 405 KB -> 12 KB - EPIRK4s3A: 636 KB -> 11 KB - EPIRK5P2: 981 KB -> 16 KB - Exp4: 1180 KB -> 23 KB - Exprb32: 250 KB -> 3.7 KB - Exprb43: 506 KB -> 8.4 KB Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

C2 is a strided reshape of the flat PhivCache buffer, which BLAS gemm consumes directly; the adapt (which copies m x (k+1)) is only needed when w has a different storage type (e.g. GPU). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…n test - expcache field on ExpvCache/PhivCache was typed ::Any; make it ::Vector{Any}. The container is now concrete; the entries stay boxed because each exponential! workspace's type is size-dependent (LinearSolve's default algorithm choice varies with matrix size) and one cache holds workspaces of several sizes at once, so a single concrete element type is not possible. Drops the now-redundant runtime isa(Vector{Any}) guard in get_expcache!. - Add a runtime allocation regression test asserting phiv!/expv!/ phiv_timestep! per-call allocations are independent of the problem size n after warmup (1488/1088/3328 bytes at n=64 and n=1024 alike), with a small absolute ceiling. This is the meaningful guard for the workspace-reuse fix: AllocCheck's static check_allocs cannot certify it, since the workspace is allocated lazily on the first cache miss and reused afterward (a temporal property), and LinearSolve/BLAS/the dynamic dispatch over the heterogeneous expcache store are flagged regardless. Uses minimum-over-repeats to filter one-off JIT allocations. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…allocation-free Addresses review feedback on the expcache typing and the residual per-call allocations. Concrete cache typing: - ExpvCache/PhivCache gain a workspace-type parameter W, and expcache is now Vector{Tuple{Int, W}} (was Vector{Any}). The exponential! workspace type is fixed by the element type and the default ExpMethodHigham2005Base and is independent of matrix size (LinearSolve's DefaultLinearSolver is a size-independent wrapper), so a cache spanning several subspace sizes stores one concrete type; only the instances (one per size) differ. W is recovered at construction via Base.promote_op (no allocation; yields Nothing for element types without an alloc_mem preallocation). - get_expcache! dispatches on the method type so its return is concretely typed (W or Nothing, never a union): the default method uses the typed store, every other method falls back to two-arg exponential!. ExpMethodHigham2005Base allocation-free (the source of the residual phiv!/expv! allocations, confirmed by comparison with the already-0-alloc non-Base ExpMethodHigham2005): - Pade coefficients are const NTuples fed through a type-stable function barrier (_pade_evaluate!) instead of a per-call "C = T[...]" vector. - Balancing uses gebal_noalloc! with a preallocated scale buffer added to the Base cache, instead of the allocating LAPACK.gebal!. Per-call allocations (KS-like operator, m=30, warm cache, Julia 1.12), released -> after: phiv! ~90 KB -> 0; expv! -> 96 B; phiv_timestep! -> 352 B. All size-independent; the size-independence allocation test's ceilings are tightened accordingly. exponential!(Higham2005Base) itself goes 464 -> 0 bytes/call, correct to machine precision across norm ranges and Float64/ComplexF64/Float32. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Moves the remaining per-call temporaries into the caches (or removes them), so the cached Krylov phi/exp entry points allocate nothing after warmup. - arnoldi! firststep!: replace the `@. V[:, 1] = b/beta` broadcast, whose materialized column view did not elide, with an explicit loop. - phiv!: the `errest ? (w, err) : w` return made the return type a union, which boxed at internal callers (phiv_timestep!). Split the body into a type-stable internal `_phiv!` that always returns `(w, err)`; the public phiv! keeps its contract. phiv_timestep! now calls `_phiv!`. - expv!: copy the first column of exp(H) into a new ExpvCache.expcol buffer instead of taking a non-eliding column view of the reshaped mem buffer. - phiv_timestep!: reuse a new PhivCache.coeffs buffer for the t^l/l! scratch instead of `ones(tType, p)` per call. - scalar-time phiv_timestep!/expv_timestep!: the matrix method now accepts a AbstractVecOrMat output and a vector B via small dispatch helpers, so the scalar wrappers pass u (and b) straight through with no reshape, and take the length-1 ts from a new PhivCache.ts1 buffer. No OrdinaryDiffEq change is needed: the new buffers live in PhivCache/ExpvCache, which the integrators construct through the public constructors and thread in via the existing caches tuple. Per-call allocations (warm cache, Julia 1.12), before -> after this commit: phiv! 0 -> 0; expv! 96 -> 0; phiv_timestep! (vector ts) 240 -> 0; phiv_timestep! (scalar t) 352 -> 0; expv_timestep! ~160 -> 0. All size-independent; the allocation regression test's ceilings are tightened. End to end this drops the EPIRK integrators' ExponentialUtilities-side per-step allocations to zero (e.g. EPIRK5P2 16 KB -> 80 B/step). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…tracker Reconciles the branch with master, which had diverged: it merged the PhivCache resize fix (SciML#254, dropped here as already-applied), added API docstrings, migrated QA to SciMLTesting 2.4, and replaced the LAPACK matrix-balancing path with GenericSchur.balance! in the Higham2005 exp methods. GenericSchur.balance! allocates ~15 KB per call even for a small (Hessenberg- sized) matrix, so on master every Krylov phiv!/expv! step allocates ~15 KB from the reduced matrix exponential. This restores a non-allocating LAPACK GEBAL wrapper (`gebal_noalloc!`) that writes the balancing scale into the cache buffer, used only for the strided-BlasFloat CPU exponential; non- BlasFloat / GPU inputs keep the GenericSchur fallback. This is what makes the Krylov exp/phi hot path allocation-free again. QA: re-add the ExplicitImports ignores for the LAPACK GEBAL symbols (@blasfunc, BlasInt, chkfinite, chklapackerror, chkstride1, libblastrampoline) and Base.promote_op, in the SciMLTesting 2.4 `ei_kwargs` form. AllocCheck: wire AllocCheck into the test suite as a documented tracker. `check_allocs` cannot report zero for these entry points and that is expected -- the first call for a given subspace size lazily builds the exponential! workspace through LinearSolve.init (hundreds of static sites), and the Pade denominator solve goes through LinearSolve.solve!; inlining attributes both to the caller and they cannot be filtered from the per-call code. The runtime size-independence testset is the authoritative allocation guard; the AllocCheck testset documents the static result (broken = true). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

ChrisRackauckas-Claude · 2026-07-22T22:47:10Z

Rebased onto master. Two things worth flagging from the rebase:

1. GenericSchur.balance! is a ~15 KB/call allocation regression on the CPU exp path. Master (via the GPU-fallback work) replaced the LAPACK balancing in the ExpMethodHigham2005/Base matrix exponentials with GenericSchur.balance!. Measured on a 33×33 Hessenberg-sized matrix:

GenericSchur.balance! alone:        15168 bytes/call
exponential!(_, Higham2005Base):    15280 bytes/call

So every Krylov phiv!/expv!/phiv_timestep! step currently allocates ~15 KB from the reduced matrix exponential on master. This PR restores a non-allocating LAPACK gebal_noalloc! (writing the scale into the cache buffer) scoped to the strided-BlasFloat CPU method only — non-BlasFloat / GPU inputs keep the GenericSchur.balance! fallback, so GPU support is unaffected. With it, exponential!(_, ExpMethodHigham2005Base, work) is 0 bytes/call and the Krylov hot path is allocation-free. Happy to split this into its own PR if you'd prefer to review the balancing change separately.

2. AllocCheck (per request): wired in as a test dependency and a testset, but it cannot report a static zero here — check_allocs finds 346 sites for phiv! (even arnoldi! has 12), all tracing into LinearSolve.init (lazy first-call workspace build) and LinearSolve.solve!, which inlining attributes to the caller and shallow backtraces can't filter out. So the AllocCheck testset is a documented tracker (broken = true, the pattern OrdinaryDiffEqExponentialRK uses), and the runtime size-independence testset is the authoritative allocation guard.

Local: Core 695 pass / 3 broken (1 pre-existing + 2 AllocCheck markers); QA 21 pass.

🤖 Generated with Claude Code

…n-gated ceiling On Julia 1.10 the concrete workspace-cache type parameter (a ~1200-character type dominated by the LinearSolve DefaultLinearSolver LinearCache) exceeds the compiler's inference-complexity limits, so the Krylov entry points box a small, size-independent constant (~8 KB/call) that they do not on 1.11+. The absolute byte-ceiling assertion therefore failed on the LTS and Downgrade jobs. Replace the byte ceiling as the primary guard with a structural reuse check that is independent of Julia/dependency version: after warmup, repeated calls must keep the PhivCache/ExpvCache flat buffer identity and not grow the exponential! workspace store. This directly tests the SciML#254/SciML#255 reallocation fixes -- which a byte count cannot reliably catch anyway, since a reallocated workspace is itself a size-independent ~15 KB constant. Keep size-independence on all versions; assert the tight byte ceiling only on Julia >= 1.11. Verified locally: full suite passes on Julia 1.10.11 (696 pass, 3 broken) and 1.12.6. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

On Julia 1.10 the ~1200-character concrete workspace-cache type defeats inference and each Krylov call boxes a size-independent ~8 KB constant, so the tight per-call byte bound does not hold there. Rather than gate the assertion out on 1.10, mark it broken = (VERSION < v"1.11"): it passes normally on 1.11+ and is recorded as an accepted broken test on 1.10 (visible, and flags if the allocation is ever fixed there). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…size-independence guards The tight per-call byte ceiling is not reliable across platforms: the per-call constant is 0 on modern x86 but a few KB wherever the compiler does not fully optimize the very large (~1200-char) concrete workspace-cache type -- Julia 1.10 everywhere, and also macOS on 1.11/pre. That constant is size-independent and is not a reallocation regression, so the ceiling produced false failures. Remove the byte ceiling. The workspace/buffer-identity structural check catches the per-call reallocation regressions (SciML#254/SciML#255) directly and portably, and size-independence (equal allocation at small and large n) catches O(n^2) growth on every version and platform. Verified locally on Julia 1.10.11 and 1.12.6. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

The explicit first-column loop added to firststep! for a non-allocating CPU path did scalar indexing on the Krylov basis, which is disallowed on GPU arrays (the GPU CI hit a scalar-indexing error via arnoldi!). Dispatch on ArrayInterface.fast_scalar_indexing: use the 0-allocation loop on CPU arrays and a broadcast into the column on GPU / non-scalar-indexing arrays. Verified on CPU: fast_scalar_indexing is true, arnoldi!/phiv_timestep! stay allocation-free. Core 696 pass / 3 broken, QA 21 pass. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…A group - ExpMethodHigham2005Base is defined for every T <: BlasFloat, but only Float64 was exercised. Add a testset covering Float64, Float32, ComplexF64 and ComplexF32 across every Pade norm branch and the scaling-squaring path, checked against a high-precision reference. Confirms the const Float64 coefficient tuples (converted to T in _pade_evaluate!) and the per-BLAS-type gebal_noalloc! methods are correct for all element types; no extra coefficients are needed. - Move the AllocCheck static-analysis testset from Core (basictests.jl, which runs on lts/1/pre) to the QA group (test/qa/qa.jl, lts/1 only) so AllocCheck is not exercised on the `pre` channel. The AllocCheck dependency moves from the main Project.toml test target to test/qa/Project.toml accordingly. Local: Core 736 pass / 1 broken; QA 21 pass / 2 broken. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…ot typeof(w) Commit 9bd1e4b (SciML#258) changed the three `Adapt.adapt(parameterless_type(w), x)` calls in the GPU expv!/phiv! paths to `Adapt.adapt(typeof(w), x)`. `typeof(w)` is the fully concrete GPU array type (e.g. CuArray{ComplexF64, 1, DeviceMemory}), whose constructor cannot take a host array, so adapting the CPU-side reduced exponential result throws `MethodError: no method matching CuArray{ComplexF64, 1, ...}(::Matrix{ComplexF64})`. `parameterless_type(w)` (e.g. CuArray) is the correct target for Adapt. Revert all three to parameterless_type (the pre-SciML#258 code) and re-import it from ArrayInterface; add the corresponding ExplicitImports ignore (parameterless_type is not public). This was surfaced here because the firststep! GPU fix lets the GPU test run get past the earlier scalar-indexing error to reach these adapt calls. Local: Core 736 pass / 1 broken, QA 21 pass / 2 broken. (GPU path itself needs the self-hosted runner; the fix restores the known-good pre-SciML#258 code.) Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

ChrisRackauckas-Claude mentioned this pull request Jul 17, 2026

ExponentialRK: remove per-step allocations from EPIRK/Exp4/EXPRB53s3 steppers SciML/OrdinaryDiffEq.jl#3928

Merged

ChrisRackauckas-Claude force-pushed the phiv-expcache-reuse branch from 2fe839e to edb7e04 Compare July 22, 2026 21:24

ChrisRackauckas added 6 commits July 22, 2026 18:07

ChrisRackauckas-Claude force-pushed the phiv-expcache-reuse branch from edb7e04 to a20834f Compare July 22, 2026 22:46

ChrisRackauckas added 6 commits July 22, 2026 19:25

ChrisRackauckas marked this pull request as ready for review July 23, 2026 08:10

ChrisRackauckas closed this Jul 23, 2026

ChrisRackauckas reopened this Jul 23, 2026

ChrisRackauckas merged commit c1be5b2 into SciML:master Jul 23, 2026
36 checks passed

Navigation Menu

Uh oh!

Uh oh!

Reuse exponential! workspaces across phiv!/expv! calls - #255

Reuse exponential! workspaces across phiv!/expv! calls#255
ChrisRackauckas merged 12 commits into
SciML:masterfrom
ChrisRackauckas-Claude:phiv-expcache-reuse

ChrisRackauckas-Claude commented Jul 17, 2026

Uh oh!

ChrisRackauckas-Claude commented Jul 22, 2026

Uh oh!

ChrisRackauckas-Claude commented Jul 22, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Uh oh!

Conversation

ChrisRackauckas-Claude commented Jul 17, 2026

Summary

Measurements (local, Julia 1.12.6)

Tests

Uh oh!

ChrisRackauckas-Claude commented Jul 22, 2026

Uh oh!

ChrisRackauckas-Claude commented Jul 22, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants