Skip to content

Use Krylov for matrix-free split ExpRK operators#3588

Open
AshtonSBradley wants to merge 1 commit into
SciML:masterfrom
AshtonSBradley:fix-matrix-free-split-exprk
Open

Use Krylov for matrix-free split ExpRK operators#3588
AshtonSBradley wants to merge 1 commit into
SciML:masterfrom
AshtonSBradley:fix-matrix-free-split-exprk

Conversation

@AshtonSBradley
Copy link
Copy Markdown

@AshtonSBradley AshtonSBradley commented May 1, 2026

Summary

Teaches fixed-step exponential RK methods to use the existing Krylov path when a SplitODEProblem linear operator is matrix-free.

This builds on the has_concretization support merged in SciMLOperators PR: SciML/SciMLOperators.jl#361

Fixes #2078.

Motivation

Operator-splitting problems using SciMLOperators matrix-free compositions can fail during cache construction because the non-Krylov path eagerly calls convert(AbstractMatrix, A).

Changes

  • During prepare_alg, detect split linear operators with !has_concretization(A).
  • Automatically reconstruct the same algorithm with krylov = true for LawsonEuler, NorsettEuler, ETDRK2, ETDRK3, ETDRK4, and HochOst4.
  • Preserve m, iop, autodiff, and concrete_jac.
  • Add a regression test for a matrix-free FunctionOperator \ DiagonalOperator * FunctionOperator split.

Type stability / allocations

The new check happens once during algorithm preparation. It does not add per-step work. Matrix-backed operators keep the existing dense cached path. Matrix-free operators move onto the already-existing Krylov path, avoiding failed matrix materialization.

Tests

On current master, with SciMLOperators v1.18.0:

julia --project=lib/OrdinaryDiffEqExponentialRK -e 'using OrdinaryDiffEqExponentialRK, LinearAlgebra, Test; using SciMLBase: successful_retcode; u0 = ComplexF64[1.0 + 0.5im, -0.5 + 0.25im, 0.75 - 0.125im, -0.25 - 0.5im]; λ = ComplexF64[-1.0, -2.0, -3.0, -4.0]; F = FunctionOperator((v, u, p, t) -> v, u0; T = ComplexF64, islinear = true, op_inverse = (v, u, p, t) -> v, opnorm = (_ -> 1.0)); L = cache_operator(F \ DiagonalOperator(λ) * F, u0); @test !has_concretization(L); prob = SplitODEProblem(L, (u, p, t) -> zero(u), u0, (0.0, 0.1)); sol = solve(prob, ETDRK4(), dt = 0.01, save_everystep = false); @test successful_retcode(sol); @test sol.alg.krylov; @test sol.u[end] ≈ exp.(0.1 .* λ) .* u0 rtol = 1.0e-6'

The OrdinaryDiffEq matrix-free split tests for this change pass on the rebased branch.

Pkg.test() for OrdinaryDiffEqExponentialRK was also rerun on Julia 1.12.5 after the rebase. The new Matrix-free SciMLOperator split test passes; the run still reports the existing ExpRK Autodiff error, which is independent of this change and left out of scope for this PR.

@AshtonSBradley
Copy link
Copy Markdown
Author

Codex was used to help prepare and validate this PR.

@AshtonSBradley AshtonSBradley marked this pull request as ready for review May 1, 2026 20:56
@AshtonSBradley AshtonSBradley marked this pull request as draft May 5, 2026 10:09
@AshtonSBradley AshtonSBradley force-pushed the fix-matrix-free-split-exprk branch from 9c337c0 to 7599142 Compare May 5, 2026 18:33
@AshtonSBradley AshtonSBradley marked this pull request as ready for review May 5, 2026 18:49
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

SplitODEProblem and SciMLOperators error

1 participant