Skip to content

Bio: isolate loser methods to show large cost without bloating main WP #1623

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ChrisRackauckas merged 1 commit into from
Jul 24, 2026
Merged

Bio: isolate loser methods to show large cost without bloating main WP #1623

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 74 additions & 61 deletions benchmarks/Bio/BCR.jmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,8 +150,7 @@ reltols = 1.0 ./ 10.0 .^ (5:8);

## Failures

Before proceeding to the results, we note the notable omissions. CVODE with KLU diverges in the solution, and
thus it is omitted from the results:
CVODE with KLU diverges on this problem and is omitted from all suites:

```julia
try
Expand All @@ -161,78 +160,37 @@ catch e
end
```

## Work-Precision Diagrams (CVODE and lsoda solvers)
## Work-Precision Diagrams (competitive CVODE / lsoda)

#### Declare solvers.
Main suite: methods that remain competitive on this large sparse stiff system.
Dense Lapack and bare GMRES (no preconditioner) are timed only in the loser
section below.

```julia
setups = [
Dict(:alg=>lsoda(), :prob_choice => 1),
Dict(:alg=>CVODE_BDF(), :prob_choice => 1),
Dict(:alg=>CVODE_BDF(linear_solver = :LapackDense), :prob_choice => 1),
Dict(:alg=>CVODE_BDF(linear_solver = :GMRES), :prob_choice => 1),
Dict(
:alg=>CVODE_BDF(linear_solver = :GMRES, prec = precilu, psetup = psetupilu, prec_side = 1),
:prob_choice => 2)
];
```

#### Plot Work-Precision Diagram.

```julia
wp = WorkPrecisionSet(
[oprob, oprob_sparse, sparsejacprob], abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = [test_sol, test_sol, test_sol], maxiters = Int(1e6), numruns = 1)

names = ["lsoda" "CVODE_BDF" "CVODE_BDF (LapackDense)" "CVODE_BDF (GMRES)" "CVODE_BDF (GMRES, iLU)" "CVODE_BDF (KLU, sparse jac)"]
plot(wp; label = names)
```

## Work-Precision Diagrams (various Julia solvers)

#### Declare solvers (using default linear solver).

```julia
setups = [
Dict(:alg=>TRBDF2(autodiff = AutoFiniteDiff())),
Dict(:alg=>QNDF(autodiff = AutoFiniteDiff())),
Dict(:alg=>FBDF(autodiff = AutoFiniteDiff())),
Dict(:alg=>KenCarp4(autodiff = AutoFiniteDiff()))
];
```

#### Plot Work-Precision Diagram (using default linear solver).

```julia
wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)

names = ["TRBDF2" "QNDF" "FBDF" "KenCarp4"]
names = ["lsoda" "CVODE_BDF" "CVODE_BDF (GMRES, iLU)"]
plot(wp; label = names)
```

#### Declare solvers (using GMRES linear solver).

```julia
setups = [
Dict(:alg=>TRBDF2(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff())),
Dict(:alg=>QNDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff())),
Dict(:alg=>FBDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff())),
Dict(:alg=>KenCarp4(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
];
```

#### Plot Work-Precision Diagram (using GMRES linear solver).
## Work-Precision Diagrams (competitive Julia solvers)

```julia
wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)

names = ["TRBDF2 (GMRES)" "QNDF (GMRES)" "FBDF (GMRES)" "KenCarp4 (GMRES)"]
plot(wp; label = names)
```
Only preconditioned GMRES and sparse KLU variants; default dense factorizations
and unpreconditioned GMRES are in the loser section.

#### Declare solvers (using GMRES linear solver, with pre-conditioner).
#### GMRES + incomplete LU

```julia
setups = [
Expand All @@ -245,8 +203,6 @@ setups = [
];
```

#### Plot Work-Precision Diagram (using GMRES linear solver, with pre-conditioner).

```julia
wp = WorkPrecisionSet(sparsejacprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)
Expand All @@ -255,9 +211,7 @@ names = ["TRBDF2 (GMRES, iLU)" "QNDF (GMRES, iLU)" "FBDF (GMRES, iLU)" "KenCarp4
plot(wp; label = names)
```

#### Declare solvers (using sparse jacobian)

We designate the solvers we wish to use.
#### Sparse Jacobian + KLU

```julia
setups = [
Expand All @@ -268,10 +222,6 @@ setups = [
];
```

#### Plot Work-Precision Diagram (using sparse jacobian)

Finally, we generate a work-precision diagram for the selection of solvers.

```julia
wp = WorkPrecisionSet(sparsejacprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e6), numruns = 1)
Expand All @@ -280,6 +230,69 @@ names = ["TRBDF2 (KLU, sparse jac)" "QNDF (KLU, sparse jac)" "FBDF (KLU, sparse
plot(wp; label = names)
```

## Loser methods (large cost in isolation)

On this ~1122-ODE sparse chemistry system, the following are not competitive:
dense Lapack factorization, default dense Julia Newton/linear solves, and GMRES
without a preconditioner. We do **not** fold them into the main work-precision
suites. Instead each is timed **once**, in isolation, at a fixed tolerance, next
to a competitive sparse reference so the wall-time gap is obvious.

```julia
const _loser_tol = 1e-6
const _loser_maxiters = Int(1e6)
_solve_kwargs = (; abstol = _loser_tol, reltol = _loser_tol, maxiters = _loser_maxiters,
save_everystep = false)

loser_labels = String[]
loser_elapsed = Float64[]

function _time_loser!(label, prob, alg)
println("--- $label ---")
t = @elapsed sol = solve(prob, alg; _solve_kwargs...)
@show sol.retcode
println("elapsed = ", t, " s")
push!(loser_labels, label)
push!(loser_elapsed, t)
return sol
end

# Competitive reference (sparse KLU)
_time_loser!("FBDF + KLU (reference)", sparsejacprob,
FBDF(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff()))

# Dense CVODE Lapack
_time_loser!("CVODE_BDF LapackDense", oprob, CVODE_BDF(linear_solver = :LapackDense))

# Bare CVODE GMRES (no preconditioner)
_time_loser!("CVODE_BDF GMRES (no prec)", oprob, CVODE_BDF(linear_solver = :GMRES))

# Default dense Julia factorizations on the non-sparse problem
_time_loser!("TRBDF2 (default dense)", oprob, TRBDF2(autodiff = AutoFiniteDiff()))
_time_loser!("QNDF (default dense)", oprob, QNDF(autodiff = AutoFiniteDiff()))
_time_loser!("FBDF (default dense)", oprob, FBDF(autodiff = AutoFiniteDiff()))
_time_loser!("KenCarp4 (default dense)", oprob, KenCarp4(autodiff = AutoFiniteDiff()))

# Unpreconditioned GMRES on the dense residual problem
_time_loser!("TRBDF2 GMRES (no prec)", oprob,
TRBDF2(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
_time_loser!("QNDF GMRES (no prec)", oprob,
QNDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
_time_loser!("FBDF GMRES (no prec)", oprob,
FBDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
_time_loser!("KenCarp4 GMRES (no prec)", oprob,
KenCarp4(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()))
```

```julia
# Relative cost vs the sparse KLU reference (first entry)
ref_t = loser_elapsed[1]
bar(loser_labels, loser_elapsed ./ ref_t; xrotation = 45, legend = false,
ylabel = "wall time / (FBDF+KLU reference)",
title = "BCR loser isolation (tol=$_loser_tol, one solve each)",
size = (900, 500), left_margin = 5Plots.mm, bottom_margin = 15Plots.mm)
```

## Summary of results

Finally, we compute a single diagram comparing the various solvers used.
Expand Down
39 changes: 32 additions & 7 deletions benchmarks/Bio/egfr_net.jmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -233,16 +233,41 @@ names = ["lsoda" "CVODE_Adams" "Tsit5" "BS5" "VCABM" "Vern6" "Vern7" "Vern8" "Ve
plot(wp; label = names)
```

#### Additional explicit solvers
## Loser methods (large cost in isolation)

One additional explicit solver, `ROCK2`, performs noticeably worse as compared to the other ones.
`ROCK2` is not competitive with the explicit set above. It is **not** folded into
the multi-tolerance / multi-run work-precision suites; a single-tolerance isolation
timing against a competitive explicit (`Vern6`) shows the wall-time gap.

```julia
setups = [Dict(:alg=>ROCK2())];
wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e9), numruns = 200)
names = ["ROCK2"]
plot(wp; label = names)
const _loser_tol = 1e-6
const _loser_maxiters = Int(1e7)
_solve_kwargs = (; abstol = _loser_tol, reltol = _loser_tol, maxiters = _loser_maxiters,
save_everystep = false)

loser_labels = String[]
loser_elapsed = Float64[]

function _time_loser!(label, alg)
println("--- $label ---")
t = @elapsed sol = solve(oprob, alg; _solve_kwargs...)
@show sol.retcode
println("elapsed = ", t, " s")
push!(loser_labels, label)
push!(loser_elapsed, t)
return sol
end

_time_loser!("Vern6 (reference)", Vern6())
_time_loser!("ROCK2 (loser)", ROCK2())
```

```julia
ref_t = loser_elapsed[1]
bar(loser_labels, loser_elapsed ./ ref_t; legend = false,
ylabel = "wall time / Vern6 reference",
title = "egfr_net ROCK2 isolation (tol=$_loser_tol, one solve each)",
size = (600, 400))
```

## Summary of results
Expand Down
39 changes: 32 additions & 7 deletions benchmarks/Bio/multisite2.jmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,16 +200,41 @@ names = ["lsoda" "CVODE_Adams" "Tsit5" "BS5" "VCABM" "Vern6" "Vern7" "Vern8" "Ve
plot(wp; label = names)
```

#### Additional explicit solvers
## Loser methods (large cost in isolation)

One additional explicit solver, `ROCK2`, performs noticeably worse as compared to the other ones.
`ROCK2` is not competitive with the explicit set above. It is **not** folded into
the multi-tolerance / multi-run work-precision suites; a single-tolerance isolation
timing against a competitive explicit (`Vern6`) shows the wall-time gap.

```julia
setups = [Dict(:alg=>ROCK2())];
wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e9), numruns = 200)
names = ["ROCK2"]
plot(wp; label = names)
const _loser_tol = 1e-6
const _loser_maxiters = Int(1e7)
_solve_kwargs = (; abstol = _loser_tol, reltol = _loser_tol, maxiters = _loser_maxiters,
save_everystep = false)

loser_labels = String[]
loser_elapsed = Float64[]

function _time_loser!(label, alg)
println("--- $label ---")
t = @elapsed sol = solve(oprob, alg; _solve_kwargs...)
@show sol.retcode
println("elapsed = ", t, " s")
push!(loser_labels, label)
push!(loser_elapsed, t)
return sol
end

_time_loser!("Vern6 (reference)", Vern6())
_time_loser!("ROCK2 (loser)", ROCK2())
```

```julia
ref_t = loser_elapsed[1]
bar(loser_labels, loser_elapsed ./ ref_t; legend = false,
ylabel = "wall time / Vern6 reference",
title = "multisite2 ROCK2 isolation (tol=$_loser_tol, one solve each)",
size = (600, 400))
```

## Summary of results
Expand Down
41 changes: 34 additions & 7 deletions benchmarks/Bio/multistate.jmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,16 +199,43 @@ names = ["lsoda" "Tsit5" "BS5" "VCABM" "Vern6" "Vern7" "Vern8" "Vern9" "ROCK4"]
plot(wp; label = names)
```

#### Additional explicit solvers
## Loser methods (large cost in isolation)

Two additional explicit solvers, `CVODE_Adams` and `ROCK2`, perform noticeably worse as compared to the other ones.
`CVODE_Adams` and `ROCK2` are not competitive with the explicit set above. They
are **not** folded into the multi-tolerance / multi-run work-precision suites; a
single-tolerance isolation timing against a competitive explicit (`Vern7`) shows
the wall-time gap.

```julia
setups = [Dict(:alg=>CVODE_Adams()), Dict(:alg=>ROCK2())];
wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate = :l2,
saveat = tf/10000.0, appxsol = test_sol, maxiters = Int(1e9), numruns = 200)
names = ["CVODE_Adams" "ROCK2"]
plot(wp; label = names)
const _loser_tol = 1e-6
const _loser_maxiters = Int(1e7)
_solve_kwargs = (; abstol = _loser_tol, reltol = _loser_tol, maxiters = _loser_maxiters,
save_everystep = false)

loser_labels = String[]
loser_elapsed = Float64[]

function _time_loser!(label, alg)
println("--- $label ---")
t = @elapsed sol = solve(oprob, alg; _solve_kwargs...)
@show sol.retcode
println("elapsed = ", t, " s")
push!(loser_labels, label)
push!(loser_elapsed, t)
return sol
end

_time_loser!("Vern7 (reference)", Vern7())
_time_loser!("CVODE_Adams (loser)", CVODE_Adams())
_time_loser!("ROCK2 (loser)", ROCK2())
```

```julia
ref_t = loser_elapsed[1]
bar(loser_labels, loser_elapsed ./ ref_t; legend = false, xrotation = 20,
ylabel = "wall time / Vern7 reference",
title = "multistate loser isolation (tol=$_loser_tol, one solve each)",
size = (700, 400), bottom_margin = 8Plots.mm)
```

## Summary of results
Expand Down
Loading