performance boosts

Author: ChrisRackauckas

src/integrators/integrator_utils.jl

@@ -32,7 +32,7 @@ end
 
 @inline function modify_dt_for_tstops!(integrator)
   tstops = integrator.opts.tstops
-  if !isempty(tstops)
+  @fastmath if !isempty(tstops)
     if integrator.opts.adaptive
       if integrator.tdir > 0
         integrator.dt = min(abs(integrator.dt),abs(top(tstops)-integrator.t)) # step! to the end
@@ -43,7 +43,7 @@ end
       integrator.dt = integrator.tdir*abs(top(tstops)-integrator.t)
     elseif integrator.dtchangeable && !integrator.force_stepfail
       # always try to step! with dtcache, but lower if a tstops
-      integrator.dt = integrator.tdir*min(abs(integrator.dtcache),abs(top(tstops)-integrator.t)) # step! to the end
+      integrator.dt = @fastmath integrator.tdir*min(abs(integrator.dtcache),abs(top(tstops)-integrator.t)) # step! to the end
     end
   end
 end
@@ -119,7 +119,8 @@ end
       end
     end
   end
-  if force_save || (integrator.opts.save_everystep && integrator.iter%integrator.opts.timeseries_steps==0)
+  if force_save || (integrator.opts.save_everystep &&
+     integrator.iter%integrator.opts.timeseries_steps==0)
     integrator.saveiter += 1
     if integrator.opts.save_idxs == nothing
       copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
@@ -155,7 +156,7 @@ end
       integrator.tprev = integrator.t
       if typeof(integrator.t)<:AbstractFloat && !isempty(integrator.opts.tstops)
         tstop = top(integrator.opts.tstops)
-        abs(ttmp - tstop) < 10eps(typeof(integrator.EEst)) ? (integrator.t = tstop) : (integrator.t = ttmp)
+        @fastmath abs(ttmp - tstop) < 10eps(typeof(integrator.EEst)) ? (integrator.t = tstop) : (integrator.t = ttmp)
       else
         integrator.t = ttmp
       end
@@ -166,7 +167,9 @@ end
     integrator.tprev = integrator.t
     if typeof(integrator.t)<:AbstractFloat && !isempty(integrator.opts.tstops)
       tstop = top(integrator.opts.tstops)
-      abs(ttmp - tstop) < 10eps(integrator.t) ? (integrator.t = tstop) : (integrator.t = ttmp)
+      # For some reason 10eps(integrator.t) is slow here
+      # TODO: Allow higher precision but profile
+      @fastmath abs(ttmp - tstop) < 1e-15 ? (integrator.t = tstop) : (integrator.t = ttmp)
     else
       integrator.t = ttmp
     end
@@ -260,7 +263,7 @@ end
   modify_dt_for_tstops!(integrator)
   accept_step!(integrator.W,integrator.dt)
   integrator.dt = integrator.W.dt
-  integrator.sqdt = sqrt(abs(integrator.dt)) # It can change dt, like in RSwM1
+  integrator.sqdt = @fastmath sqrt(abs(integrator.dt)) # It can change dt, like in RSwM1
 end
 
 @inline function handle_tstop!(integrator)

src/perform_step/low_order.jl

@@ -1,9 +1,9 @@
 @muladd function perform_step!(integrator,cache::EMConstantCache,f=integrator.f)
   @unpack t,dt,uprev,u,W = integrator
-  if is_diagonal_noise(integrator.sol.prob)
-    noise = integrator.g(t,uprev).*W.dW
-  else
+  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
     noise = integrator.g(t,uprev)*W.dW
+  else
+    noise = integrator.g(t,uprev).*W.dW
   end
   u = @muladd uprev + dt*integrator.f(t,uprev) + noise
   integrator.u = u

src/perform_step/sra.jl

@@ -191,7 +191,7 @@ end
 
   if integrator.opts.adaptive
     E₁ = dt*(k1 + k2 + k3)
-    integrator.EEst = integrator.opts.internalnorm(@muladd(integrator.opts.delta*E₁+E₂)./@muladd(integrator.opts.abstol + max.(integrator.opts.internalnorm.(uprev),integrator.opts.internalnorm.(u))*integrator.opts.reltol))
+    integrator.EEst = integrator.opts.internalnorm(@muladd(integrator.opts.delta*E₁+E₂)/@muladd(integrator.opts.abstol + max(integrator.opts.internalnorm(uprev),integrator.opts.internalnorm(u))*integrator.opts.reltol))
   end
   integrator.u = u
 end

src/solve.jl

@@ -14,10 +14,11 @@ function init(
   recompile::Type{Val{recompile_flag}}=Val{true};
   dt = tType(0),
   timeseries_steps::Int = 1,
-  save_noise = true,
   saveat = tType[],tstops = tType[],d_discontinuities= tType[],
   save_timeseries = nothing,
   save_everystep = isempty(saveat),
+  save_noise = save_everystep && typeof(prob.f) <: Tuple ?
+               has_analytic(prob.f[1]) : has_analytic(prob.f),
   save_idxs = nothing,
   save_start = true,save_end = true,
   dense = save_everystep,
@@ -151,6 +152,16 @@ function init(
     timeseries = convert(Vector{typeof(u_initial)},timeseries_init)
   end
   ts = convert(Vector{tType},ts_init)
+
+  if !adaptive
+    dt == 0 ? steps = length(tstops) : steps = round(Int,(tspan[2]-tspan[1])/dt,RoundUp)
+    sizehint!(timeseries,steps+1)
+    sizehint!(ts,steps+1)
+  else
+    sizehint!(timeseries,50)
+    sizehint!(ts,50)
+  end
+
   #ks = convert(Vector{ksEltype},ks_init)
   alg_choice = Int[]
 
@@ -220,31 +231,33 @@ function init(
   seed == 0 ? (prob.seed == 0 ? _seed = rand(UInt64) : _seed = prob.seed) : _seed = seed
 
   if typeof(prob.noise) <: Void
+    isadaptive(alg) ? rswm = RSWM(adaptivealg=:RSwM3) : rswm = RSWM(adaptivealg=:RSwM1)
     if isinplace
-      isadaptive(alg) ? rswm = RSWM(adaptivealg=:RSwM3) : rswm = RSWM(adaptivealg=:RSwM1)
       if alg_needs_extra_process(alg)
         W = WienerProcess!(t,rand_prototype,rand_prototype,
-                           save_everystep=save_everystep,
+                           save_everystep=save_noise,
                            timeseries_steps=timeseries_steps,
                            rswm=rswm,
                            rng = Xorshifts.Xoroshiro128Plus(_seed))
       else
         W = WienerProcess!(t,rand_prototype,
-                           save_everystep=save_everystep,
+                           save_everystep=save_noise,
                            timeseries_steps=timeseries_steps,
                            rswm=rswm,
                            rng = Xorshifts.Xoroshiro128Plus(_seed))
       end
     else
       if alg_needs_extra_process(alg)
         W = WienerProcess(t,rand_prototype,rand_prototype,
-                           save_everystep=save_everystep,
+                           save_everystep=save_noise,
                            timeseries_steps=timeseries_steps,
+                           rswm=rswm,
                            rng = Xorshifts.Xoroshiro128Plus(_seed))
       else
         W = WienerProcess(t,rand_prototype,
-                           save_everystep=save_everystep,
+                           save_everystep=save_noise,
                            timeseries_steps=timeseries_steps,
+                           rswm=rswm,
                            rng = Xorshifts.Xoroshiro128Plus(_seed))
       end
     end

test/first_rand_test.jl

@@ -4,21 +4,21 @@ f1(t,u) = 0.
 g1(t,u) = 1.
 dt = 1//2^(4)
 prob1 = SDEProblem{false}(f1,g1,0.,(0.0,1.0))
-integrator = init(prob1,EM(),dt=dt)
+integrator = init(prob1,EM(),dt=dt,save_noise=true)
 
 k = integrator.W.dW
 @test integrator.W.dW != 0
 solve!(integrator)
 @test integrator.sol.W[2] == k
 
 prob1 = SDEProblem{false}(f1,g1,zeros(4),(0.0,1.0))
-sol = solve(prob1,EM(),dt=dt)
+sol = solve(prob1,EM(),dt=dt,save_noise=true)
 @test sol.W[2] != zeros(4)
 
 
 f1(t,u,du) = du.=0.
 g1(t,u,du) = du.=1.
 dt = 1//2^(4)
 prob1 = SDEProblem(f1,g1,zeros(4),(0.0,1.0))
-sol = solve(prob1,EM(),dt=dt)
+sol = solve(prob1,EM(),dt=dt,save_noise=true)
 @test sol.W[2] != zeros(4)

test/split_tests.jl

@@ -8,7 +8,7 @@ f2(t,u) = (1.01)/2 * u
 
 prob = SDEProblem{false}((f1,f2),σ,1/2,(0.0,1.0))
 
-sol = solve(prob,SplitEM(),dt=1/10)
+sol = solve(prob,SplitEM(),dt=1/10,save_noise=true)
 
 prob = SDEProblem{false}(f,σ,1/2,(0.0,1.0),noise = NoiseWrapper(sol.W))
 
@@ -19,7 +19,7 @@ sol2 = solve(prob,EM(),dt=1/10)
 u0 = rand(4)
 prob = SDEProblem{false}((f1,f2),σ,u0,(0.0,1.0))
 
-sol = solve(prob,SplitEM(),dt=1/10)
+sol = solve(prob,SplitEM(),dt=1/10,save_noise=true)
 
 prob = SDEProblem{false}(f,σ,u0,(0.0,1.0),noise = NoiseWrapper(sol.W))