using OrdinaryDiffEqLowStorageRK, DiffEqDevTools, Test, Random
import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear

Random.seed!(100)

testTol = 0.25

f = (u, p, t) -> cos(t)
prob_ode_sin = ODEProblem(ODEFunction(f; analytic = (u0, p, t) -> sin(t)), 0.0, (0.0, 1.0))

f = (du, u, p, t) -> du[1] = cos(t)
prob_ode_sin_inplace = ODEProblem(
    ODEFunction(f; analytic = (u0, p, t) -> [sin(t)]), [0.0],
    (0.0, 1.0)
)

f = (u, p, t) -> sin(u)
prob_ode_nonlinear = ODEProblem(
    ODEFunction(
        f;
        analytic = (u0, p, t) -> 2 * acot(
            exp(-t) *
                cot(0.5)
        )
    ), 1.0,
    (0.0, 0.5)
)

f = (du, u, p, t) -> du[1] = sin(u[1])
prob_ode_nonlinear_inplace = ODEProblem(
    ODEFunction(
        f;
        analytic = (u0, p, t) -> [
            2 * acot(exp(-t) * cot(0.5)),
        ]
    ),
    [1.0], (0.0, 0.5)
)

test_problems_only_time = [prob_ode_sin, prob_ode_sin_inplace]
test_problems_linear = [prob_ode_linear, prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear]
test_problems_nonlinear = [prob_ode_nonlinear, prob_ode_nonlinear_inplace]

# Test the memory usage, cf. #640
# Note: Basically, the size of the integrator should be the size of the cache
# plus the size of the initial condition (stored is integ.sol.prob.u0) if the
# keyword argument `alias_u0` is not set to `true` (default).
# Note: The memory requirements of the 2N methods can be reduced if an assignment
# of the form `tmp = A2end[i]*tmp + dt*f(u, p, t+c2end[i]*dt)` can be carried out
# without saving `f(u, p, t+c2end[i]*dt)` as `k`.
u0_large = rand(10^6)
prob_ode_large = ODEProblem((du, u, p, t) -> du .= u, u0_large, (0.0, 1.0))

println("SHLDDRK_2N")
dts_SHLDDRK_2N = (1 / 2) .^ (0:3)
alg = SHLDDRK_2N()
for prob in test_problems_only_time
    sim = test_convergence(dts_SHLDDRK_2N, prob, alg)
    @test sim.𝒪est[:final] ≈ 4 atol = 0.46
end
for prob in test_problems_linear
    sim = test_convergence(dts_SHLDDRK_2N, prob, alg)
    @test sim.𝒪est[:final] ≈ 4 atol = 0.46
end
for prob in test_problems_nonlinear
    sim = test_convergence(dts_SHLDDRK_2N, prob, alg)
    @test sim.𝒪est[:final] ≈ 4 atol = 1
    # due to unusual saturation towards high dts(0.5 and onwards) and
    # saturation towards low dts due to less precision in the provided values of weights , tolerance is kept so high
end

println("SHLDDRK52")
dts = 1 .// 2 .^ (8:-1:4)
alg = SHLDDRK52()
for prob in test_problems_only_time
    sim = test_convergence(dts, prob, alg)
    @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
end
for prob in test_problems_linear
    sim = test_convergence(dts, prob, alg)
    @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
end
for prob in test_problems_nonlinear
    sim = test_convergence(dts, prob, alg)
    @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
end

@testset "ORK256" begin
    alg = ORK256(; williamson_condition = true)
    alg2 = ORK256(; williamson_condition = false)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "RK46NL" begin
    alg = RK46NL(; williamson_condition = true)
    alg2 = RK46NL(; williamson_condition = false)
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CarpenterKennedy2N54" begin
    alg = CarpenterKennedy2N54(; williamson_condition = true)
    alg2 = CarpenterKennedy2N54(; williamson_condition = false)
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end


@testset "SHLDDRK64" begin
    alg = SHLDDRK64(; williamson_condition = true)
    alg2 = SHLDDRK64(; williamson_condition = true)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "DGLDDRK73_C" begin
    alg = DGLDDRK73_C(; williamson_condition = true)
    alg2 = DGLDDRK73_C(; williamson_condition = false)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "DGLDDRK84_C" begin
    alg = DGLDDRK84_C(; williamson_condition = true)
    alg2 = DGLDDRK84_C(; williamson_condition = false)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "DGLDDRK84_F" begin
    alg = DGLDDRK84_F(; williamson_condition = true)
    alg2 = DGLDDRK84_F(; williamson_condition = false)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "NDBLSRK124" begin
    alg = NDBLSRK124(; williamson_condition = true)
    alg2 = NDBLSRK124(; williamson_condition = false)
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "NDBLSRK134" begin
    alg = NDBLSRK134(; williamson_condition = true)
    alg2 = NDBLSRK134(; williamson_condition = false)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "NDBLSRK144" begin
    alg = NDBLSRK144(; williamson_condition = true)
    alg2 = NDBLSRK144(; williamson_condition = false)
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
        sim = test_convergence(dts, prob, alg2)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CFRLDDRK64" begin
    alg = CFRLDDRK64()
    dts = 1 ./ 2 .^ (7:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "TSLDDRK74" begin
    alg = TSLDDRK74()
    dts = 1 ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

# Methods from Carpenter, Kennedy, Lewis (2000)

function RemakeNew(p::ODEProblem)
    u1 = @. BigFloat(p.u0)
    tsp1 = @. BigFloat(p.tspan)
    return remake(p; u0 = u1, tspan = tsp1)
end

test_problems_only_time_BigFloat = @. RemakeNew(test_problems_only_time)
test_problems_linear_BigFloat = @. RemakeNew(test_problems_linear)
f = (u, p, t) -> sin(u)
prob_nonlinear_A = ODEProblem(
    ODEFunction(
        f;
        analytic = (u0, p, t) -> 2 * acot(
            exp(-t) *
                cot(BigFloat(0.5))
        )
    ),
    BigFloat(1.0), (BigFloat(0.0), BigFloat(0.5))
)

f = (du, u, p, t) -> du[1] = sin(u[1])
prob_nonlinear_B = ODEProblem(
    ODEFunction(
        f;
        analytic = (u0, p, t) -> [
            2 * acot(exp(-t) * cot(BigFloat(0.5))),
        ]
    ),
    [BigFloat(1.0)],
    (BigFloat(0.0), BigFloat(0.5))
)
test_problems_nonlinear_BigFloat = [prob_nonlinear_A, prob_nonlinear_B]

@testset "CKLLSRK43_2" begin
    alg = CKLLSRK43_2()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol    # This scheme has linear order of 4
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK54_3C" begin
    alg = CKLLSRK54_3C()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ 1 atol = testTol          # The CI plot is linear but the evaluated order is 1
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK95_4S" begin
    alg = CKLLSRK95_4S()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK95_4C" begin
    alg = CKLLSRK95_4C()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK95_4M" begin
    alg = CKLLSRK95_4M()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK54_3C_3R" begin
    alg = CKLLSRK54_3C_3R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK54_3M_3R" begin
    alg = CKLLSRK54_3M_3R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 0.5 atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK54_3N_3R" begin
    alg = CKLLSRK54_3N_3R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK85_4C_3R" begin
    alg = CKLLSRK85_4C_3R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK85_4M_3R" begin
    alg = CKLLSRK85_4M_3R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK85_4P_3R" begin
    alg = CKLLSRK85_4P_3R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 2 atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK54_3N_4R" begin
    alg = CKLLSRK54_3N_4R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK54_3M_4R" begin
    alg = CKLLSRK54_3M_4R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 0.5 atol = testTol                              # This scheme has linear orderof 4.5
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK65_4M_4R" begin
    alg = CKLLSRK65_4M_4R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK85_4FM_4R" begin
    alg = CKLLSRK85_4FM_4R()
    dts = BigFloat(1) ./ 2 .^ (10:-1:6)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "CKLLSRK75_4M_5R" begin
    alg = CKLLSRK75_4M_5R()
    dts = BigFloat(1) ./ 2 .^ (8:-1:4)
    for prob in test_problems_only_time_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear_BigFloat
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, adaptive = false, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 13
    integ = init(
        prob_ode_large, alg, adaptive = true, dt = 1.0e-2, save_start = false,
        save_end = false, save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 14
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 13
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

# Methods from Parsani, Ketcheson, Deconinck (2013)

@testset "ParsaniKetchesonDeconinck3S32" begin
    alg = ParsaniKetchesonDeconinck3S32()
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S82" begin
    alg = ParsaniKetchesonDeconinck3S82()
    dts = 1 ./ 2 .^ (8:-1:5)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S53" begin
    alg = ParsaniKetchesonDeconinck3S53()
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S173" begin
    alg = ParsaniKetchesonDeconinck3S173()
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1 ./ 2 .^ (6:-1:3)
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = 1
    end

    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S94" begin
    alg = ParsaniKetchesonDeconinck3S94()
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S184" begin
    alg = ParsaniKetchesonDeconinck3S184()
    dts = 1 ./ 2 .^ (6:-1:2)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1 ./ 2 .^ (7:-1:2)
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S105" begin
    alg = ParsaniKetchesonDeconinck3S105()
    dts = 1 ./ 1.95 .^ (5:-1:1)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1 ./ 2 .^ (5:-1:2)
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1.5 ./ 2 .^ (5:-1:2)
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "ParsaniKetchesonDeconinck3S205" begin
    alg = ParsaniKetchesonDeconinck3S205()
    dts = 1 ./ 1.95 .^ (5:-1:1)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1 ./ 2 .^ (5:-1:2)
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = 0.33
    end
    dts = 1.5 ./ 2 .^ (5:-1:2)
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

# Methods from Al Jahdali, Dalcin, Boukharfane, Nolasco, Keyes, Parsani (2022)

@testset "AlJahdali 3S* methods" begin
    aljahdali_algs = [
        AlJahdaliAdv3S42(), AlJahdaliAdv3S82(), AlJahdaliAdv3S53(), AlJahdaliAdv3S113(),
        AlJahdaliAdv3S64(), AlJahdaliAdv3S154(), AlJahdaliAdv3S85(), AlJahdaliAdv3S165(),
        AlJahdaliVor3S42(), AlJahdaliVor3S82(), AlJahdaliVor3S53(), AlJahdaliVor3S113(),
        AlJahdaliVor3S64(), AlJahdaliVor3S154(), AlJahdaliVor3S85(), AlJahdaliVor3S165(),
    ]
    for alg in aljahdali_algs
        order = OrdinaryDiffEqLowStorageRK.alg_order(alg)
        # order-5 methods saturate at machine precision on the fine end of the
        # standard dts range (same effect documented for ParsaniKetchesonDeconinck3S105/205
        # above), so they get the same coarser dts.
        dts = order == 5 ? 1 ./ 1.95 .^ (5:-1:1) : 1 ./ 2 .^ (7:-1:3)
        for prob in test_problems_only_time
            sim = test_convergence(dts, prob, alg)
            # atol = 1.25, not testTol: unlike ParsaniKetcheson3S*, these don't show a
            # consistent +1 pure-quadrature bump (some do, some measure close to their
            # claimed order, some land between) at this dts range, so a tight per-method
            # fit isn't reliable here. This still catches an actually wrong coefficient
            # set, which would land a full order or more off.
            @test sim.𝒪est[:final] ≈ order atol = 1.25
        end
    end
end

# Methods from Ranocha, Dalcin, Parsani, Ketcheson (2021)

@testset "RDPK3Sp35" begin
    alg = RDPK3Sp35()
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "RDPK3Sp49" begin
    alg = RDPK3Sp49()
    dts = 1 ./ 2 .^ (5:-1:2)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1 ./ 2 .^ (8:-1:2)
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "RDPK3Sp510" begin
    alg = RDPK3Sp510()
    dts = 1 ./ 2 .^ (4.5:-1:1.5)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "RDPK3SpFSAL35" begin
    alg = RDPK3SpFSAL35()
    dts = 1 ./ 2 .^ (7:-1:3)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "RDPK3SpFSAL49" begin
    alg = RDPK3SpFSAL49()
    dts = 1 ./ 2 .^ (5:-1:2)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        # higher order as pure quadrature
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) + 1 atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    dts = 1 ./ 2 .^ (8:-1:2)
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

@testset "RDPK3SpFSAL510" begin
    alg = RDPK3SpFSAL510()
    dts = 1 ./ 2 .^ (4.5:-1:1.5)
    for prob in test_problems_only_time
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_linear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    for prob in test_problems_nonlinear
        sim = test_convergence(dts, prob, alg)
        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
    end
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
    integ = init(
        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
    # test whether aliasing u0 is bad
    new_prob_ode_nonlinear_inplace = ODEProblem(
        prob_ode_nonlinear_inplace.f, [1.0],
        (0.0, 0.5)
    )
    sol_old = solve(
        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false
    )
    sol_new = solve(
        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
    )
    @test sol_old.u[end] ≈ sol_new.u[end]
end

## Convergence Testing
testTol = 0.2

@testset "Implicit Solver Convergence Tests ($(["out-of-place", "in-place"][i]))" for i in 1:2
    prob = (
        prob_ode_linear,
        prob_ode_2Dlinear,
    )[i]
    sim20 = test_convergence(dts, prob, RK46NL())
    @test sim20.𝒪est[:final] ≈ 4 atol = testTol
    @test sim20.𝒪est[:l2] ≈ 4 atol = testTol
    @test sim20.𝒪est[:l∞] ≈ 4 atol = testTol
end

@testset "VectorOfArray/StructArray compatibility" begin
    using RecursiveArrayTools, StaticArrays, StructArrays

    function rhs!(du_voa, u_voa, p, t)
        du = parent(du_voa)
        u = parent(u_voa)
        du .= u
    end

    # StructArray storage
    u = StructArray{SVector{1, Float64}}(ntuple(_ -> [1.0, 2.0], 1))
    ode = ODEProblem(rhs!, VectorOfArray(u), (0, 0.7))
    sol_SA = solve(ode, RDPK3SpFSAL35())

    # Vector{<:SVector} storage
    u = SVector{1, Float64}.([1.0, 2.0])
    ode = ODEProblem(rhs!, VectorOfArray(u), (0, 0.7))
    sol_SV = solve(ode, RDPK3SpFSAL35())

    @test sol_SA ≈ sol_SV
    @test sol_SV.stats.naccept == sol_SA.stats.naccept

    # Plain vector
    u = [1.0, 2.0]
    ode = ODEProblem(rhs!, u, (0, 0.7))
    sol = solve(ode, RDPK3SpFSAL35())
    @test sol.stats.naccept == sol_SA.stats.naccept
    @test sol.t ≈ sol_SA.t
    for i in eachindex(sol_SA.u), j in eachindex(u)
        @test sol.u[i][j] ≈ sol_SA.u[i][j][1]
    end
end
