Merge pull request #2498 from SciML/fastpower

Change Core, Extrapolation, and FIRK to using FastPower.jl

Author: ChrisRackauckas

lib/OrdinaryDiffEqCore/Project.toml

@@ -14,6 +14,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
+FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -53,6 +54,7 @@ DocStringExtensions = "0.9"
 EnumX = "1"
 FastBroadcast = "0.2, 0.3"
 FastClosures = "0.3"
+FastPower = "1"
 FillArrays = "1.9"
 FunctionWrappersWrappers = "0.1"
 InteractiveUtils = "1.9"

lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl

@@ -19,6 +19,8 @@ using PrecompileTools
 
 import FillArrays: Trues, Falses
 
+import FastPower
+
 # Interfaces
 import DiffEqBase: solve!, step!, initialize!, isadaptive
 

lib/OrdinaryDiffEqCore/src/integrators/controllers.jl

@@ -68,7 +68,7 @@ end
         q = inv(qmax)
     else
         expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
-        qtmp = DiffEqBase.fastpow(EEst, expo) / gamma
+        qtmp = FastPower.fastpower(EEst, expo) / gamma
         @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
         # TODO: Shouldn't this be in `step_accept_controller!` as for the PI controller?
         integrator.qold = DiffEqBase.value(integrator.dt) / q
@@ -141,8 +141,8 @@ end
     if iszero(EEst)
         q = inv(qmax)
     else
-        q11 = DiffEqBase.fastpow(EEst, float(beta1))
-        q = q11 / DiffEqBase.fastpow(qold, float(beta2))
+        q11 = FastPower.fastpower(EEst, float(beta1))
+        q = q11 / FastPower.fastpower(qold, float(beta2))
         integrator.q11 = q11
         @fastmath q = max(inv(qmax), min(inv(qmin), q / gamma))
     end

lib/OrdinaryDiffEqExtrapolation/Project.toml

@@ -6,6 +6,7 @@ version = "1.1.0"
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
@@ -18,6 +19,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
 FastBroadcast = "0.3.5"
+FastPower = "1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"
 OrdinaryDiffEqCore = "1.1"

lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl

@@ -20,6 +20,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or
                            differentiation_rk_docstring
 using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools, LinearSolve
 import OrdinaryDiffEqCore
+import FastPower
 import OrdinaryDiffEqDifferentiation: TimeDerivativeWrapper, UDerivativeWrapper, calc_J,
                                       WOperator, TimeGradientWrapper, UJacobianWrapper,
                                       build_grad_config,

lib/OrdinaryDiffEqExtrapolation/src/controllers.jl

@@ -31,10 +31,10 @@ function stepsize_controller_internal!(integrator,
     else
         # Update gamma and beta1
         controller.beta1 = typeof(controller.beta1)(1 // (2integrator.cache.n_curr + 1))
-        integrator.opts.gamma = DiffEqBase.fastpow(typeof(integrator.opts.gamma)(1 // 4),
+        integrator.opts.gamma = FastPower.fastpower(typeof(integrator.opts.gamma)(1 // 4),
             controller.beta1)
         # Compute new stepsize scaling
-        qtmp = DiffEqBase.fastpow(integrator.EEst, controller.beta1) / integrator.opts.gamma
+        qtmp = FastPower.fastpower(integrator.EEst, controller.beta1) / integrator.opts.gamma
         @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
     end
     integrator.cache.Q[integrator.cache.n_curr - alg.min_order + 1] = q
@@ -57,10 +57,10 @@ function stepsize_predictor!(integrator,
         s_new = stage_number[n_new - alg.min_order + 1]
         # Update gamma and beta1
         controller.beta1 = typeof(controller.beta1)(1 // (2integrator.cache.n_curr + 1))
-        integrator.opts.gamma = DiffEqBase.fastpow(typeof(integrator.opts.gamma)(1 // 4),
+        integrator.opts.gamma = FastPower.fastpower(typeof(integrator.opts.gamma)(1 // 4),
             controller.beta1)
         # Compute new stepsize scaling
-        qtmp = EEst * DiffEqBase.fastpow(DiffEqBase.fastpow(tol, (1.0 - s_curr / s_new)),
+        qtmp = EEst * FastPower.fastpower(FastPower.fastpower(tol, (1.0 - s_curr / s_new)),
             controller.beta1) / integrator.opts.gamma
         @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
     end
@@ -158,12 +158,12 @@ function stepsize_controller_internal!(integrator,
                 controller.beta1 = typeof(controller.beta1)(1 //
                                                             (integrator.cache.n_curr - 1))
             end
-            integrator.opts.gamma = DiffEqBase.fastpow(
+            integrator.opts.gamma = FastPower.fastpower(
                 typeof(integrator.opts.gamma)(65 //
                                               100),
                 controller.beta1)
             # Compute new stepsize scaling
-            qtmp = DiffEqBase.fastpow(integrator.EEst, controller.beta1) /
+            qtmp = FastPower.fastpower(integrator.EEst, controller.beta1) /
                    (integrator.opts.gamma)
             @fastmath q = max(inv(integrator.opts.qmax),
                 min(inv(integrator.opts.qmin), qtmp))
@@ -175,12 +175,12 @@ function stepsize_controller_internal!(integrator,
         else
             # Update gamma and beta1
             controller.beta1 = typeof(controller.beta1)(1 // (2integrator.cache.n_curr + 1))
-            integrator.opts.gamma = DiffEqBase.fastpow(
+            integrator.opts.gamma = FastPower.fastpower(
                 typeof(integrator.opts.gamma)(65 //
                                               100),
                 controller.beta1)
             # Compute new stepsize scaling
-            qtmp = DiffEqBase.fastpow(integrator.EEst, controller.beta1) /
+            qtmp = FastPower.fastpower(integrator.EEst, controller.beta1) /
                    integrator.opts.gamma
             @fastmath q = max(inv(integrator.opts.qmax),
                 min(inv(integrator.opts.qmin), qtmp))

lib/OrdinaryDiffEqFIRK/Project.toml

@@ -6,6 +6,7 @@ version = "1.1.1"
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -25,6 +26,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
 FastBroadcast = "0.3.5"
+FastPower = "1"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"

lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl

@@ -24,6 +24,7 @@ using LinearAlgebra: I, UniformScaling, mul!, lu
 import LinearSolve
 import FastBroadcast: @..
 import OrdinaryDiffEqCore
+import FastPower
 
 using OrdinaryDiffEqDifferentiation: UJacobianWrapper, build_J_W, build_jac_config,
                                      UDerivativeWrapper, calc_J!, dolinsolve, calc_J,

lib/OrdinaryDiffEqFIRK/src/controllers.jl

@@ -17,7 +17,7 @@
             fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
         end
         expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
-        qtmp = DiffEqBase.fastpow(EEst, expo) / fac
+        qtmp = FastPower.fastpower(EEst, expo) / fac
         @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
         integrator.qold = q
     end
@@ -31,7 +31,7 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
     if integrator.success_iter > 0
         expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
         qgus = (integrator.dtacc / integrator.dt) *
-               DiffEqBase.fastpow((EEst^2) / integrator.erracc, expo)
+               FastPower.fastpower((EEst^2) / integrator.erracc, expo)
         qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
         qacc = max(q, qgus)
     else

test/interface/composite_algorithm_test.jl

@@ -47,11 +47,11 @@ v = @inferred OrdinaryDiffEqCore.ode_extrapolant(
     alg_mixed_r = CompositeAlgorithm((ABM54(), Tsit5()), reverse_choice)
     alg_mixed2 = CompositeAlgorithm((Tsit5(), ABM54()), reverse_choice)
 
-    @test_throws ErrorException solve(prob_ode_linear, alg_mixed)
+    @test_throws ArgumentError solve(prob_ode_linear, alg_mixed)
     sol2 = solve(prob_ode_linear, Tsit5())
-    sol3 = solve(prob_ode_linear, alg_mixed; dt = 0.05)
-    sol4 = solve(prob_ode_linear, alg_mixed_r; dt = 0.05)
-    sol5 = solve(prob_ode_linear, alg_mixed2; dt = 0.05)
+    sol3 = solve(prob_ode_linear, alg_mixed; dt = 0.05, adaptive=false)
+    sol4 = solve(prob_ode_linear, alg_mixed_r; dt = 0.05, adaptive=false)
+    sol5 = solve(prob_ode_linear, alg_mixed2; dt = 0.05, adaptive=false)
     @test sol3.t == sol4.t && sol3.u == sol4.u
     @test sol3(0.8)≈sol2(0.8) atol=1e-4
     @test sol5(0.8)≈sol2(0.8) atol=1e-4