From 5ed41ce579e9a299ecf776c295a5bc3702abc846 Mon Sep 17 00:00:00 2001
From: Michele Mesiti <michele.mesiti@kit.edu>
Date: Wed, 4 Dec 2024 00:40:17 +0100
Subject: [PATCH 1/3] Added threading to ABM algorithms using @..

1. in *_perform_step.jl
   in functions using @..
   in unpacking of the "cache" object
   add a "thread" variable
2. in *_caches.jl
   for the cache types affected
   (determined in 1.)
   - add parametric type Thread
   - and a field thread::Thread
   - in the alg_cache function,
     add 'alg.thread' as argument to the cache type constructor
3. in *_algorithms.jl,
   for the alg type affected
   (determined in 2.)
   - add parametric type Thread
   - add field thread::Thread
   - add outer constructor taking a thread argument
     with default value False(),
     that builds the respective algorithm object
     forwarding the thread parameter to it.
---
 .../src/adams_bashforth_moulton_caches.jl     |  67 ++++---
 .../adams_bashforth_moulton_perform_step.jl   | 165 +++++++++---------
 .../src/algorithms.jl                         |  93 ++++++++--
 3 files changed, 207 insertions(+), 118 deletions(-)

diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl
index d417f9cdad..283b4b6b18 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl
@@ -4,7 +4,7 @@ get_fsalfirstlast(cache::ABMMutableCache, u) = (cache.fsalfirst, cache.k)
 function get_fsalfirstlast(cache::ABMVariableCoefficientMutableCache, u)
     (cache.fsalfirst, cache.k4)
 end
-@cache mutable struct AB3Cache{uType, rateType} <: ABMMutableCache
+@cache mutable struct AB3Cache{uType, rateType, Thread} <: ABMMutableCache
     u::uType
     uprev::uType
     fsalfirst::rateType
@@ -14,6 +14,7 @@ end
     k::rateType
     tmp::uType
     step::Int
+    thread::Thread
 end
 
 @cache mutable struct AB3ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
@@ -32,7 +33,7 @@ function alg_cache(alg::AB3, u, rate_prototype, ::Type{uEltypeNoUnits},
     ralk2 = zero(rate_prototype)
     k = zero(rate_prototype)
     tmp = zero(u)
-    AB3Cache(u, uprev, fsalfirst, k2, k3, ralk2, k, tmp, 1)
+    AB3Cache(u, uprev, fsalfirst, k2, k3, ralk2, k, tmp, 1, alg.thread)
 end
 
 function alg_cache(alg::AB3, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -44,7 +45,7 @@ function alg_cache(alg::AB3, u, rate_prototype, ::Type{uEltypeNoUnits},
     AB3ConstantCache(k2, k3, 1)
 end
 
-@cache mutable struct ABM32Cache{uType, rateType} <: ABMMutableCache
+@cache mutable struct ABM32Cache{uType, rateType, Thread} <: ABMMutableCache
     u::uType
     uprev::uType
     fsalfirst::rateType
@@ -54,6 +55,7 @@ end
     k::rateType
     tmp::uType
     step::Int
+    thread::Thread
 end
 
 @cache mutable struct ABM32ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
@@ -72,7 +74,7 @@ function alg_cache(alg::ABM32, u, rate_prototype, ::Type{uEltypeNoUnits},
     ralk2 = zero(rate_prototype)
     k = zero(rate_prototype)
     tmp = zero(u)
-    ABM32Cache(u, uprev, fsalfirst, k2, k3, ralk2, k, tmp, 1)
+    ABM32Cache(u, uprev, fsalfirst, k2, k3, ralk2, k, tmp, 1, alg.thread)
 end
 
 function alg_cache(alg::ABM32, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -84,7 +86,7 @@ function alg_cache(alg::ABM32, u, rate_prototype, ::Type{uEltypeNoUnits},
     ABM32ConstantCache(k2, k3, 1)
 end
 
-@cache mutable struct AB4Cache{uType, rateType} <: ABMMutableCache
+@cache mutable struct AB4Cache{uType, rateType, Thread} <: ABMMutableCache
     u::uType
     uprev::uType
     fsalfirst::rateType
@@ -98,6 +100,7 @@ end
     t3::rateType
     t4::rateType
     step::Int
+    thread::Thread
 end
 
 @cache mutable struct AB4ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
@@ -121,7 +124,7 @@ function alg_cache(alg::AB4, u, rate_prototype, ::Type{uEltypeNoUnits},
     t2 = zero(rate_prototype)
     t3 = zero(rate_prototype)
     t4 = zero(rate_prototype)
-    AB4Cache(u, uprev, fsalfirst, k2, k3, k4, ralk2, k, tmp, t2, t3, t4, 1)
+    AB4Cache(u, uprev, fsalfirst, k2, k3, k4, ralk2, k, tmp, t2, t3, t4, 1, alg.thread)
 end
 
 function alg_cache(alg::AB4, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -134,7 +137,7 @@ function alg_cache(alg::AB4, u, rate_prototype, ::Type{uEltypeNoUnits},
     AB4ConstantCache(k2, k3, k4, 1)
 end
 
-@cache mutable struct ABM43Cache{uType, rateType} <: ABMMutableCache
+@cache mutable struct ABM43Cache{uType, rateType, Thread} <: ABMMutableCache
     u::uType
     uprev::uType
     fsalfirst::rateType
@@ -151,6 +154,7 @@ end
     t6::rateType
     t7::rateType
     step::Int
+    thread::Thread
 end
 
 @cache mutable struct ABM43ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
@@ -177,7 +181,8 @@ function alg_cache(alg::ABM43, u, rate_prototype, ::Type{uEltypeNoUnits},
     t5 = zero(rate_prototype)
     t6 = zero(rate_prototype)
     t7 = zero(rate_prototype)
-    ABM43Cache(u, uprev, fsalfirst, k2, k3, k4, ralk2, k, tmp, t2, t3, t4, t5, t6, t7, 1)
+    ABM43Cache(u, uprev, fsalfirst, k2, k3, k4, ralk2, k,
+        tmp, t2, t3, t4, t5, t6, t7, 1, alg.thread)
 end
 
 function alg_cache(alg::ABM43, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -190,7 +195,7 @@ function alg_cache(alg::ABM43, u, rate_prototype, ::Type{uEltypeNoUnits},
     ABM43ConstantCache(k2, k3, k4, 1)
 end
 
-@cache mutable struct AB5Cache{uType, rateType} <: ABMMutableCache
+@cache mutable struct AB5Cache{uType, rateType, Thread} <: ABMMutableCache
     u::uType
     uprev::uType
     fsalfirst::rateType
@@ -204,6 +209,7 @@ end
     t3::rateType
     t4::rateType
     step::Int
+    thread::Thread
 end
 
 @cache mutable struct AB5ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
@@ -228,7 +234,7 @@ function alg_cache(alg::AB5, u, rate_prototype, ::Type{uEltypeNoUnits},
     t2 = zero(rate_prototype)
     t3 = zero(rate_prototype)
     t4 = zero(rate_prototype)
-    AB5Cache(u, uprev, fsalfirst, k2, k3, k4, k5, k, tmp, t2, t3, t4, 1)
+    AB5Cache(u, uprev, fsalfirst, k2, k3, k4, k5, k, tmp, t2, t3, t4, 1, alg.thread)
 end
 
 function alg_cache(alg::AB5, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -242,7 +248,7 @@ function alg_cache(alg::AB5, u, rate_prototype, ::Type{uEltypeNoUnits},
     AB5ConstantCache(k2, k3, k4, k5, 1)
 end
 
-@cache mutable struct ABM54Cache{uType, rateType} <: ABMMutableCache
+@cache mutable struct ABM54Cache{uType, rateType, Thread} <: ABMMutableCache
     u::uType
     uprev::uType
     fsalfirst::rateType
@@ -260,6 +266,7 @@ end
     t7::rateType
     t8::rateType
     step::Int
+    thread::Thread
 end
 
 @cache mutable struct ABM54ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
@@ -288,7 +295,8 @@ function alg_cache(alg::ABM54, u, rate_prototype, ::Type{uEltypeNoUnits},
     t6 = zero(rate_prototype)
     t7 = zero(rate_prototype)
     t8 = zero(rate_prototype)
-    ABM54Cache(u, uprev, fsalfirst, k2, k3, k4, k5, k, tmp, t2, t3, t4, t5, t6, t7, t8, 1)
+    ABM54Cache(u, uprev, fsalfirst, k2, k3, k4, k5, k, tmp,
+        t2, t3, t4, t5, t6, t7, t8, 1, alg.thread)
 end
 
 function alg_cache(alg::ABM54, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -317,7 +325,7 @@ end
 end
 
 @cache mutable struct VCAB3Cache{uType, rateType, TabType, bs3Type, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -337,6 +345,7 @@ end
     utilde::uType
     tab::TabType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCAB3, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -395,7 +404,7 @@ function alg_cache(alg::VCAB3, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp = zero(u)
     utilde = zero(u)
     VCAB3Cache(u, uprev, fsalfirst, bs3cache, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕstar_n, β,
-        order, atmp, tmp, utilde, tab, 1)
+        order, atmp, tmp, utilde, tab, 1, alg.thread)
 end
 
 @cache mutable struct VCAB4ConstantCache{rk4constcache, tArrayType, rArrayType, cArrayType,
@@ -413,7 +422,7 @@ end
 end
 
 @cache mutable struct VCAB4Cache{uType, rateType, rk4cacheType, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -432,6 +441,7 @@ end
     tmp::uType
     utilde::uType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCAB4, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -489,7 +499,7 @@ function alg_cache(alg::VCAB4, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp = zero(u)
     utilde = zero(u)
     VCAB4Cache(u, uprev, fsalfirst, rk4cache, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕstar_n, β,
-        order, atmp, tmp, utilde, 1)
+        order, atmp, tmp, utilde, 1, alg.thread)
 end
 
 # VCAB5
@@ -509,7 +519,7 @@ end
 end
 
 @cache mutable struct VCAB5Cache{uType, rateType, rk4cacheType, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -528,6 +538,7 @@ end
     tmp::uType
     utilde::uType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCAB5, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -585,7 +596,7 @@ function alg_cache(alg::VCAB5, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp = zero(u)
     utilde = zero(u)
     VCAB5Cache(u, uprev, fsalfirst, rk4cache, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕstar_n, β,
-        order, atmp, tmp, utilde, 1)
+        order, atmp, tmp, utilde, 1, alg.thread)
 end
 
 # VCABM3
@@ -607,7 +618,7 @@ end
 
 @cache mutable struct VCABM3Cache{
     uType, rateType, TabType, bs3Type, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -628,6 +639,7 @@ end
     utilde::uType
     tab::TabType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCABM3, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -691,7 +703,7 @@ function alg_cache(alg::VCABM3, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp = zero(u)
     utilde = zero(u)
     VCABM3Cache(u, uprev, fsalfirst, bs3cache, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕ_np1,
-        ϕstar_n, β, order, atmp, tmp, utilde, tab, 1)
+        ϕstar_n, β, order, atmp, tmp, utilde, tab, 1, alg.thread)
 end
 
 # VCABM4
@@ -713,7 +725,7 @@ end
 end
 
 @cache mutable struct VCABM4Cache{uType, rateType, rk4cacheType, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -733,6 +745,7 @@ end
     tmp::uType
     utilde::uType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCABM4, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -796,7 +809,7 @@ function alg_cache(alg::VCABM4, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp = zero(u)
     utilde = zero(u)
     VCABM4Cache(u, uprev, fsalfirst, rk4cache, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕ_np1,
-        ϕstar_n, β, order, atmp, tmp, utilde, 1)
+        ϕstar_n, β, order, atmp, tmp, utilde, 1, alg.thread)
 end
 
 # VCABM5
@@ -818,7 +831,7 @@ end
 end
 
 @cache mutable struct VCABM5Cache{uType, rateType, rk4cacheType, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -838,6 +851,7 @@ end
     tmp::uType
     utilde::uType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCABM5, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -901,7 +915,7 @@ function alg_cache(alg::VCABM5, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp = zero(u)
     utilde = zero(u)
     VCABM5Cache(u, uprev, fsalfirst, rk4cache, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕ_np1,
-        ϕstar_n, β, order, atmp, tmp, utilde, 1)
+        ϕstar_n, β, order, atmp, tmp, utilde, 1, alg.thread)
 end
 
 # VCABM
@@ -924,7 +938,7 @@ end
 end
 
 @cache mutable struct VCABMCache{uType, rateType, dtType, tArrayType, cArrayType,
-    uNoUnitsType, coefType, dtArrayType} <:
+    uNoUnitsType, coefType, dtArrayType, Thread} <:
                       ABMVariableCoefficientMutableCache
     u::uType
     uprev::uType
@@ -952,6 +966,7 @@ end
     atmpm2::uNoUnitsType
     atmpp1::uNoUnitsType
     step::Int
+    thread::Thread
 end
 
 function alg_cache(alg::VCABM, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -1023,5 +1038,5 @@ function alg_cache(alg::VCABM, u, rate_prototype, ::Type{uEltypeNoUnits},
     VCABMCache(
         u, uprev, fsalfirst, k4, ϕstar_nm1, dts, c, g, ϕ_n, ϕ_np1, ϕstar_n, β, order,
         max_order, atmp, tmp, ξ, ξ0, utilde, utildem1, utildem2, utildep1, atmpm1,
-        atmpm2, atmpp1, 1)
+        atmpm2, atmpp1, 1, alg.thread)
 end
diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl
index 077fdfb4ec..5e160b2b02 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl
@@ -68,7 +68,7 @@ end
 
 @muladd function perform_step!(integrator, cache::AB3Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack tmp, fsalfirst, k2, k3, ralk2, k = cache
+    @unpack tmp, fsalfirst, k2, k3, ralk2, k, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -77,17 +77,17 @@ end
     if cache.step <= 2
         cache.step += 1
         ttmp = t + (2 / 3) * dt
-        @.. broadcast=false tmp=uprev + (2 / 3) * dt * k1
+        @.. broadcast=false thread=thread tmp=uprev + (2 / 3) * dt * k1
         f(ralk2, tmp, p, ttmp)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-        @.. broadcast=false u=uprev + (dt / 4) * (k1 + 3 * ralk2)        #Ralston Method
+        @.. broadcast=false thread=thread u=uprev + (dt / 4) * (k1 + 3 * ralk2)        #Ralston Method
         if cnt == 1
             cache.k3 .= k1
         else
             cache.k2 .= k1
         end
     else
-        @.. broadcast=false u=uprev + (dt / 12) * (23 * k1 - 16 * k2 + 5 * k3)
+        @.. broadcast=false thread=thread u=uprev + (dt / 12) * (23 * k1 - 16 * k2 + 5 * k3)
         cache.k2, cache.k3 = k3, k2
         cache.k2 .= k1
     end
@@ -128,7 +128,7 @@ end
 
 @muladd function perform_step!(integrator, cache::ABM32Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack tmp, fsalfirst, k2, k3, ralk2, k = cache
+    @unpack tmp, fsalfirst, k2, k3, ralk2, k, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -137,21 +137,21 @@ end
     if cache.step == 1
         cache.step += 1
         ttmp = t + (2 / 3) * dt
-        @.. broadcast=false tmp=uprev + (2 / 3) * dt * k1
+        @.. broadcast=false thread=thread tmp=uprev + (2 / 3) * dt * k1
         f(ralk2, tmp, p, ttmp)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-        @.. broadcast=false u=uprev + (dt / 4) * (k1 + 3 * ralk2)       #Ralston Method
+        @.. broadcast=false thread=thread u=uprev + (dt / 4) * (k1 + 3 * ralk2)       #Ralston Method
         cache.k2 .= k1
     else
         if cnt == 2
             perform_step!(integrator,
-                AB3Cache(u, uprev, fsalfirst, copy(k2), k3, ralk2, k, tmp, cnt))  #Here passing copy of k2, otherwise it will change in AB3()
+                AB3Cache(u, uprev, fsalfirst, copy(k2), k3, ralk2, k, tmp, cnt, thread))  #Here passing copy of k2, otherwise it will change in AB3()
         else
             perform_step!(integrator,
-                AB3Cache(u, uprev, fsalfirst, k2, k3, ralk2, k, tmp, cnt))
+                AB3Cache(u, uprev, fsalfirst, k2, k3, ralk2, k, tmp, cnt, thread))
         end
         k = integrator.fsallast
-        @.. broadcast=false u=uprev + (dt / 12) * (5 * k + 8 * k1 - k2)
+        @.. broadcast=false thread=thread u=uprev + (dt / 12) * (5 * k + 8 * k1 - k2)
         cache.k2, cache.k3 = k3, k2
         cache.k2 .= k1
     end
@@ -198,7 +198,7 @@ end
 
 @muladd function perform_step!(integrator, cache::AB4Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack tmp, fsalfirst, k2, k3, k4, ralk2, k, t2, t3, t4 = cache
+    @unpack tmp, fsalfirst, k2, k3, k4, ralk2, k, t2, t3, t4, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -208,14 +208,14 @@ end
         cache.step += 1
         halfdt = dt / 2
         ttmp = t + halfdt
-        @.. broadcast=false tmp=uprev + halfdt * k1
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * k1
         f(t2, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + halfdt * t2
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * t2
         f(t3, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + dt * t3
+        @.. broadcast=false thread=thread tmp=uprev + dt * t3
         f(t4, tmp, p, t + dt)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3)
-        @.. broadcast=false u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
+        @.. broadcast=false thread=thread u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
         if cnt == 1
             cache.k4 .= k1
         elseif cnt == 2
@@ -224,7 +224,9 @@ end
             cache.k2 .= k1
         end
     else
-        @.. broadcast=false u=uprev + (dt / 24) * (55 * k1 - 59 * k2 + 37 * k3 - 9 * k4)
+        @.. broadcast=false thread=thread u=uprev +
+                                            (dt / 24) *
+                                            (55 * k1 - 59 * k2 + 37 * k3 - 9 * k4)
         cache.k4, cache.k3 = k3, k4
         cache.k3 .= k2
         cache.k2 .= k1
@@ -272,7 +274,7 @@ end
 
 @muladd function perform_step!(integrator, cache::ABM43Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack tmp, fsalfirst, k2, k3, k4, ralk2, k, t2, t3, t4, t5, t6, t7 = cache
+    @unpack tmp, fsalfirst, k2, k3, k4, ralk2, k, t2, t3, t4, t5, t6, t7, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -282,14 +284,14 @@ end
         cache.step += 1
         halfdt = dt / 2
         ttmp = t + halfdt
-        @.. broadcast=false tmp=uprev + halfdt * k1
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * k1
         f(t2, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + halfdt * t2
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * t2
         f(t3, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + dt * t3
+        @.. broadcast=false thread=thread tmp=uprev + dt * t3
         f(t4, tmp, p, t + dt)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3)
-        @.. broadcast=false u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
+        @.. broadcast=false thread=thread u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
         if cnt == 1
             cache.k3 .= k1
         else
@@ -301,9 +303,10 @@ end
         t4 .= k4
         perform_step!(integrator,
             AB4Cache(u, uprev, fsalfirst, t2, t3, t4, ralk2, k, tmp, t5, t6, t7,
-                cnt))
+                cnt,thread))
         k = integrator.fsallast
-        @.. broadcast=false u=uprev + (dt / 24) * (9 * k + 19 * k1 - 5 * k2 + k3)
+        @.. broadcast=false thread=thread u=uprev +
+                                            (dt / 24) * (9 * k + 19 * k1 - 5 * k2 + k3)
         cache.k4, cache.k3 = k3, k4
         cache.k3 .= k2
         cache.k2 .= k1
@@ -354,7 +357,7 @@ end
 
 @muladd function perform_step!(integrator, cache::AB5Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack tmp, fsalfirst, k2, k3, k4, k5, k, t2, t3, t4 = cache
+    @unpack tmp, fsalfirst, k2, k3, k4, k5, k, t2, t3, t4, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -364,14 +367,14 @@ end
         cache.step += 1
         halfdt = dt / 2
         ttmp = t + halfdt
-        @.. broadcast=false tmp=uprev + halfdt * k1
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * k1
         f(t2, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + halfdt * t2
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * t2
         f(t3, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + dt * t3
+        @.. broadcast=false thread=thread tmp=uprev + dt * t3
         f(t4, tmp, p, t + dt)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-        @.. broadcast=false u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
+        @.. broadcast=false thread=thread u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
         if cnt == 1
             cache.k5 .= k1
         elseif cnt == 2
@@ -382,9 +385,10 @@ end
             cache.k2 .= k1
         end
     else
-        @.. broadcast=false u=uprev +
-                              (dt / 720) *
-                              (1901 * k1 - 2774 * k2 + 2616 * k3 - 1274 * k4 + 251 * k5)
+        @.. broadcast=false thread=thread u=uprev +
+                                            (dt / 720) *
+                                            (1901 * k1 - 2774 * k2 + 2616 * k3 - 1274 * k4 +
+                                             251 * k5)
         cache.k5, cache.k4 = k4, k5
         cache.k4 .= k3
         cache.k3 .= k2
@@ -436,7 +440,7 @@ end
 
 @muladd function perform_step!(integrator, cache::ABM54Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack tmp, fsalfirst, k2, k3, k4, k5, k, t2, t3, t4, t5, t6, t7, t8 = cache
+    @unpack tmp, fsalfirst, k2, k3, k4, k5, k, t2, t3, t4, t5, t6, t7, t8, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -446,13 +450,13 @@ end
         cache.step += 1
         halfdt = dt / 2
         ttmp = t + halfdt
-        @.. broadcast=false tmp=uprev + halfdt * k1
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * k1
         f(t2, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + halfdt * t2
+        @.. broadcast=false thread=thread tmp=uprev + halfdt * t2
         f(t3, tmp, p, ttmp)
-        @.. broadcast=false tmp=uprev + dt * t3
+        @.. broadcast=false thread=thread tmp=uprev + dt * t3
         f(t4, tmp, p, t + dt)
-        @.. broadcast=false u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
+        @.. broadcast=false thread=thread u=uprev + (dt / 6) * (2 * (t2 + t3) + (k1 + t4))   #RK4
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3)
         if cnt == 1
             cache.k4 .= k1
@@ -468,11 +472,12 @@ end
         t5 .= k5
         perform_step!(integrator,
             AB5Cache(u, uprev, fsalfirst, t2, t3, t4, t5, k, tmp, t6, t7, t8,
-                cnt))
+                cnt, thread))
         k = integrator.fsallast
-        @.. broadcast=false u=uprev +
-                              (dt / 720) *
-                              (251 * k + 646 * k1 - 264 * k2 + 106 * k3 - 19 * k4)
+        @.. broadcast=false thread=thread u=uprev +
+                                            (dt / 720) *
+                                            (251 * k + 646 * k1 - 264 * k2 + 106 * k3 -
+                                             19 * k4)
         cache.k5, cache.k4 = k4, k5
         cache.k4 .= k3
         cache.k3 .= k2
@@ -561,7 +566,7 @@ end
 
 @muladd function perform_step!(integrator, cache::VCAB3Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k4, dts, g, ϕstar_n, ϕstar_nm1, order, atmp, utilde, bs3cache = cache
+    @unpack k4, dts, g, ϕstar_n, ϕstar_nm1, order, atmp, utilde, bs3cache, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -587,12 +592,12 @@ end
         integrator.fsallast .= k4
     else
         g_coefs!(cache, k)
-        @.. broadcast=false u=uprev
+        @.. broadcast=false thread=thread u=uprev
         for i in 1:k
-            @.. broadcast=false u+=g[i] * ϕstar_n[i]
+            @.. broadcast=false thread=thread u+=g[i] * ϕstar_n[i]
         end
         if integrator.opts.adaptive
-            @.. broadcast=false utilde=g[k] * ϕstar_n[k]    # Using lower order AB from subset of coefficients
+            @.. broadcast=false thread=thread utilde=g[k] * ϕstar_n[k]    # Using lower order AB from subset of coefficients
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             integrator.EEst = integrator.opts.internalnorm(atmp, t)
@@ -694,7 +699,7 @@ end
 
 @muladd function perform_step!(integrator, cache::VCAB4Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k4, dts, g, ϕ_n, ϕstar_n, ϕstar_nm1, order, atmp, utilde, rk4cache = cache
+    @unpack k4, dts, g, ϕ_n, ϕstar_n, ϕstar_nm1, order, atmp, utilde, rk4cache, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -726,12 +731,12 @@ end
         integrator.fsallast .= rk4cache.k
     else
         g_coefs!(cache, k)
-        @.. broadcast=false u=uprev
+        @.. broadcast=false thread=thread u=uprev
         for i in 1:k
-            @.. broadcast=false u+=g[i] * ϕstar_n[i]
+            @.. broadcast=false thread=thread u+=g[i] * ϕstar_n[i]
         end
         if integrator.opts.adaptive
-            @.. broadcast=false utilde=g[k] * ϕstar_n[k]
+            @.. broadcast=false thread=thread utilde=g[k] * ϕstar_n[k]
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             integrator.EEst = integrator.opts.internalnorm(atmp, t)
@@ -842,7 +847,7 @@ end
 
 @muladd function perform_step!(integrator, cache::VCAB5Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k4, dts, g, ϕ_n, ϕstar_n, ϕstar_nm1, order, atmp, utilde, rk4cache = cache
+    @unpack k4, dts, g, ϕ_n, ϕstar_n, ϕstar_nm1, order, atmp, utilde, rk4cache, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -881,12 +886,12 @@ end
         integrator.fsallast .= rk4cache.k
     else
         g_coefs!(cache, k)
-        @.. broadcast=false u=uprev
+        @.. broadcast=false thread=thread u=uprev
         for i in 1:k
-            @.. broadcast=false u+=g[i] * ϕstar_n[i]
+            @.. broadcast=false thread=thread u+=g[i] * ϕstar_n[i]
         end
         if integrator.opts.adaptive
-            @.. broadcast=false utilde=g[k] * ϕstar_n[k]
+            @.. broadcast=false thread=thread utilde=g[k] * ϕstar_n[k]
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             integrator.EEst = integrator.opts.internalnorm(atmp, t)
@@ -988,7 +993,7 @@ end
 
 @muladd function perform_step!(integrator, cache::VCABM3Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k4, dts, g, ϕstar_n, ϕ_np1, ϕstar_nm1, order, atmp, utilde, bs3cache = cache
+    @unpack k4, dts, g, ϕstar_n, ϕ_np1, ϕstar_nm1, order, atmp, utilde, bs3cache, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -1014,16 +1019,16 @@ end
         integrator.fsallast .= k4
     else
         g_coefs!(cache, k + 1)
-        @.. broadcast=false u=uprev
+        @.. broadcast=false thread=thread u=uprev
         for i in 1:(k - 1)
-            @.. broadcast=false u+=g[i] * ϕstar_n[i]
+            @.. broadcast=false thread=thread u+=g[i] * ϕstar_n[i]
         end
         f(k4, u, p, t + dt)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
         ϕ_np1!(cache, k4, k + 1)
-        @.. broadcast=false u+=g[end - 1] * ϕ_np1[end - 1]
+        @.. broadcast=false thread=thread u+=g[end - 1] * ϕ_np1[end - 1]
         if integrator.opts.adaptive
-            @.. broadcast=false utilde=(g[end] - g[end - 1]) * ϕ_np1[end]
+            @.. broadcast=false thread=thread utilde=(g[end] - g[end - 1]) * ϕ_np1[end]
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             integrator.EEst = integrator.opts.internalnorm(atmp, t)
@@ -1131,7 +1136,7 @@ end
 
 @muladd function perform_step!(integrator, cache::VCABM4Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack k4, dts, g, ϕstar_n, ϕ_np1, ϕstar_nm1, order, atmp, utilde, rk4cache = cache
+    @unpack k4, dts, g, ϕstar_n, ϕ_np1, ϕstar_nm1, order, atmp, utilde, rk4cache, thread = cache
     k1 = integrator.fsalfirst
     if integrator.u_modified
         cache.step = 1
@@ -1163,16 +1168,16 @@ end
         integrator.fsallast .= rk4cache.k
     else
         g_coefs!(cache, k + 1)
-        @.. broadcast=false u=uprev
+        @.. broadcast=false thread=thread u=uprev
         for i in 1:(k - 1)
-            @.. broadcast=false u+=g[i] * ϕstar_n[i]
+            @.. broadcast=false thread=thread u+=g[i] * ϕstar_n[i]
         end
         f(k4, u, p, t + dt)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
         ϕ_np1!(cache, k4, k + 1)
-        @.. broadcast=false u+=g[end - 1] * ϕ_np1[end - 1]
+        @.. broadcast=false thread=thread u+=g[end - 1] * ϕ_np1[end - 1]
         if integrator.opts.adaptive
-            @.. broadcast=false utilde=(g[end] - g[end - 1]) * ϕ_np1[end]
+            @.. broadcast=false thread=thread utilde=(g[end] - g[end - 1]) * ϕ_np1[end]
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             integrator.EEst = integrator.opts.internalnorm(atmp, t)
@@ -1288,7 +1293,7 @@ end
 @muladd function perform_step!(integrator, cache::VCABM5Cache, repeat_step = false)
     @inbounds begin
         @unpack t, dt, uprev, u, f, p = integrator
-        @unpack k4, dts, g, ϕ_n, ϕ_np1, ϕstar_n, ϕstar_nm1, order, atmp, utilde, rk4cache = cache
+        @unpack k4, dts, g, ϕ_n, ϕ_np1, ϕstar_n, ϕstar_nm1, order, atmp, utilde, rk4cache, thread = cache
         k1 = integrator.fsalfirst
         if integrator.u_modified
             cache.step = 1
@@ -1327,16 +1332,16 @@ end
             integrator.fsallast .= rk4cache.k
         else
             g_coefs!(cache, 6)
-            @.. broadcast=false u=muladd(g[1], ϕstar_n[1], uprev)
-            @.. broadcast=false u=muladd(g[2], ϕstar_n[2], u)
-            @.. broadcast=false u=muladd(g[3], ϕstar_n[3], u)
-            @.. broadcast=false u=muladd(g[4], ϕstar_n[4], u)
+            @.. broadcast=false thread=thread u=muladd(g[1], ϕstar_n[1], uprev)
+            @.. broadcast=false thread=thread u=muladd(g[2], ϕstar_n[2], u)
+            @.. broadcast=false thread=thread u=muladd(g[3], ϕstar_n[3], u)
+            @.. broadcast=false thread=thread u=muladd(g[4], ϕstar_n[4], u)
             f(k4, u, p, t + dt)
             OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
             ϕ_np1!(cache, k4, 6)
-            @.. broadcast=false u=muladd(g[6 - 1], ϕ_np1[6 - 1], u)
+            @.. broadcast=false thread=thread u=muladd(g[6 - 1], ϕ_np1[6 - 1], u)
             if integrator.opts.adaptive
-                @.. broadcast=false utilde=(g[6] - g[6 - 1]) * ϕ_np1[end]
+                @.. broadcast=false thread=thread utilde=(g[6] - g[6 - 1]) * ϕ_np1[end]
                 calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                     integrator.opts.reltol, integrator.opts.internalnorm,
                     t)
@@ -1464,7 +1469,7 @@ end
 @muladd function perform_step!(integrator, cache::VCABMCache, repeat_step = false)
     @inbounds begin
         @unpack t, dt, uprev, u, f, p = integrator
-        @unpack k4, dts, g, ϕ_n, ϕ_np1, ϕstar_n, ϕstar_nm1, order, max_order, utilde, utildem2, utildem1, utildep1, atmp, atmpm1, atmpm2, atmpp1 = cache
+        @unpack k4, dts, g, ϕ_n, ϕ_np1, ϕstar_n, ϕstar_nm1, order, max_order, utilde, utildem2, utildem1, utildep1, atmp, atmpm1, atmpm2, atmpp1, thread = cache
         k1 = integrator.fsalfirst
         step = integrator.iter
         k = order
@@ -1476,16 +1481,16 @@ end
         ϕ_and_ϕstar!(cache, k1, k)
         g_coefs!(cache, k + 1)
         # unroll the predictor once
-        @.. broadcast=false u=muladd(g[1], ϕstar_n[1], uprev)
+        @.. broadcast=false thread=thread u=muladd(g[1], ϕstar_n[1], uprev)
         for i in 2:(k - 1)
-            @.. broadcast=false u=muladd(g[i], ϕstar_n[i], u)
+            @.. broadcast=false thread=thread u=muladd(g[i], ϕstar_n[i], u)
         end
         f(k4, u, p, t + dt)
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
         ϕ_np1!(cache, k4, k + 1)
-        @.. broadcast=false u=muladd(g[k], ϕ_np1[k], u)
+        @.. broadcast=false thread=thread u=muladd(g[k], ϕ_np1[k], u)
         if integrator.opts.adaptive
-            @.. broadcast=false utilde=(g[k + 1] - g[k]) * ϕ_np1[k + 1]
+            @.. broadcast=false thread=thread utilde=(g[k + 1] - g[k]) * ϕ_np1[k + 1]
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             integrator.EEst = integrator.opts.internalnorm(atmp, t)
@@ -1501,13 +1506,15 @@ end
             if step <= 4 || order < 3
                 cache.order = min(order + 1, 3)
             else
-                # @.. broadcast=false utildem2 = dt * γstar[(k-2)+1] * ϕ_np1[k-1]
-                @.. broadcast=false utildem2=(g[k - 1] - g[k - 2]) * ϕ_np1[k - 1]
-                # @.. broadcast=false utildem1 = dt * γstar[(k-1)+1] * ϕ_np1[k]
-                @.. broadcast=false utildem1=(g[k] - g[k - 1]) * ϕ_np1[k]
+                # @.. broadcast=false thread=thread utildem2 = dt * γstar[(k-2)+1] * ϕ_np1[k-1]
+                @.. broadcast=false thread=thread utildem2=(g[k - 1] - g[k - 2]) *
+                                                           ϕ_np1[k - 1]
+                # @.. broadcast=false thread=thread utildem1 = dt * γstar[(k-1)+1] * ϕ_np1[k]
+                @.. broadcast=false thread=thread utildem1=(g[k] - g[k - 1]) * ϕ_np1[k]
                 expand_ϕ_and_ϕstar!(cache, k + 1)
                 ϕ_np1!(cache, k4, k + 2)
-                @.. broadcast=false utildep1=dt * γstar[(k + 1) + 1] * ϕ_np1[k + 2]
+                @.. broadcast=false thread=thread utildep1=dt * γstar[(k + 1) + 1] *
+                                                           ϕ_np1[k + 2]
                 calculate_residuals!(atmpm2, utildem2, uprev, u, integrator.opts.abstol,
                     integrator.opts.reltol, integrator.opts.internalnorm,
                     t)
diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/algorithms.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/algorithms.jl
index af35c61140..cf3b3613dc 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/algorithms.jl
@@ -10,7 +10,13 @@ reference = """E. Hairer, S. P. Norsett, G. Wanner, Solving Ordinary Differentia
     reference,
     "",
     "")
-struct AB3 <: OrdinaryDiffEqAlgorithm end
+struct AB3{Thread} <: OrdinaryDiffEqAlgorithm
+    thread::Thread
+end
+
+function AB3()
+    AB3(False())
+end
 
 @doc generic_solver_docstring("The 4-step fourth order multistep method.
     Runge-Kutta method of order 4 is used to calculate starting values.",
@@ -19,7 +25,13 @@ struct AB3 <: OrdinaryDiffEqAlgorithm end
     reference,
     "",
     "")
-struct AB4 <: OrdinaryDiffEqAlgorithm end
+struct AB4{Thread} <: OrdinaryDiffEqAlgorithm
+    thread::Thread
+end
+function AB4()
+    AB4(False())
+end
+
 @doc generic_solver_docstring("The 5-step fifth order multistep method.
     Ralston's 3rd order Runge-Kutta method is used to calculate starting values.",
     "AB5",
@@ -27,7 +39,12 @@ struct AB4 <: OrdinaryDiffEqAlgorithm end
     reference,
     "",
     "")
-struct AB5 <: OrdinaryDiffEqAlgorithm end
+struct AB5{Thread} <: OrdinaryDiffEqAlgorithm
+    thread::Thread
+end
+function AB5()
+    AB5(False())
+end
 
 @doc generic_solver_docstring("It is third order method.
     In ABM32, AB3 works as predictor and Adams Moulton 2-steps method works as Corrector.
@@ -37,7 +54,12 @@ struct AB5 <: OrdinaryDiffEqAlgorithm end
     reference,
     "",
     "")
-struct ABM32 <: OrdinaryDiffEqAlgorithm end
+struct ABM32{Thread} <: OrdinaryDiffEqAlgorithm
+    thread::Thread
+end
+function ABM32()
+    ABM32(False())
+end
 
 @doc generic_solver_docstring("It is fourth order method.
     In ABM43, AB4 works as predictor and Adams Moulton 3-steps method works as Corrector.
@@ -47,7 +69,12 @@ struct ABM32 <: OrdinaryDiffEqAlgorithm end
     reference,
     "",
     "")
-struct ABM43 <: OrdinaryDiffEqAlgorithm end
+struct ABM43{Thread} <: OrdinaryDiffEqAlgorithm
+    thread::Thread
+end
+function ABM43()
+    ABM43(False())
+end
 
 @doc generic_solver_docstring("It is fifth order method.
     In ABM54, AB5 works as predictor and Adams Moulton 4-steps method works as Corrector.
@@ -57,7 +84,12 @@ struct ABM43 <: OrdinaryDiffEqAlgorithm end
     reference,
     "",
     "")
-struct ABM54 <: OrdinaryDiffEqAlgorithm end
+struct ABM54{Thread} <: OrdinaryDiffEqAlgorithm
+    thread::Thread
+end
+function ABM54()
+    ABM54(False())
+end
 
 # Variable Step Size Adams methods
 
@@ -68,7 +100,12 @@ struct ABM54 <: OrdinaryDiffEqAlgorithm end
     reference,
     "",
     "")
-struct VCAB3 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct VCAB3{Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    thread::Thread
+end
+function VCAB3()
+    VCAB3(False())
+end
 
 @doc generic_solver_docstring("The 4th order Adams method.
     Runge-Kutta 4 is used to calculate starting values.",
@@ -77,7 +114,12 @@ struct VCAB3 <: OrdinaryDiffEqAdaptiveAlgorithm end
     reference,
     "",
     "")
-struct VCAB4 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct VCAB4{Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    thread::Thread
+end
+function VCAB4()
+    VCAB4(False())
+end
 
 @doc generic_solver_docstring("The 5th order Adams method.
     Runge-Kutta 4 is used to calculate starting values.",
@@ -86,7 +128,12 @@ struct VCAB4 <: OrdinaryDiffEqAdaptiveAlgorithm end
     reference,
     "",
     "")
-struct VCAB5 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct VCAB5{Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    thread::Thread
+end
+function VCAB5()
+    VCAB5(False())
+end
 
 @doc generic_solver_docstring("The 3rd order Adams-Moulton method.
     Bogacki-Shampine 3/2 method is used to calculate starting values.",
@@ -95,7 +142,12 @@ struct VCAB5 <: OrdinaryDiffEqAdaptiveAlgorithm end
     reference,
     "",
     "")
-struct VCABM3 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct VCABM3{Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    thread::Thread
+end
+function VCABM3()
+    VCABM3(False())
+end
 
 @doc generic_solver_docstring("The 4th order Adams-Moulton method.
     Runge-Kutta 4 is used to calculate starting values.",
@@ -104,7 +156,12 @@ struct VCABM3 <: OrdinaryDiffEqAdaptiveAlgorithm end
     reference,
     "",
     "")
-struct VCABM4 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct VCABM4{Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    thread::Thread
+end
+function VCABM4()
+    VCABM4(False())
+end
 
 @doc generic_solver_docstring("The 5th order Adams-Moulton method.
     Runge-Kutta 4 is used to calculate starting values.",
@@ -113,7 +170,12 @@ struct VCABM4 <: OrdinaryDiffEqAdaptiveAlgorithm end
     reference,
     "",
     "")
-struct VCABM5 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct VCABM5{Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    thread::Thread
+end
+function VCABM5()
+    VCABM5(False())
+end
 
 # Variable Order and Variable Step Size Adams methods
 
@@ -124,4 +186,9 @@ struct VCABM5 <: OrdinaryDiffEqAdaptiveAlgorithm end
     reference,
     "",
     "")
-struct VCABM <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm end
+struct VCABM{Thread} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm
+    thread::Thread
+end
+function VCABM()
+    VCABM(False())
+end

From 23c054b86faeda01d9323f10afac8659b0e75a06 Mon Sep 17 00:00:00 2001
From: Michele Mesiti <michele.mesiti@kit.edu>
Date: Wed, 18 Dec 2024 11:55:51 +0100
Subject: [PATCH 2/3] Add regression tests for multithreaded ABM algs

---
 .../test/regression_test_threading.jl         | 27 +++++++++++++++++++
 .../test/runtests.jl                          |  1 +
 2 files changed, 28 insertions(+)
 create mode 100644 lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl

diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl
new file mode 100644
index 0000000000..bfc0954c57
--- /dev/null
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl
@@ -0,0 +1,27 @@
+using OrdinaryDiffEqAdamsBashforthMoulton, ODEProblemLibrary
+using Test
+using Static
+
+free_timestep_algorithms = [ VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5, VCABM]
+
+fixed_timestep_algorithms = [AB3,AB4, AB5, ABM32, ABM43, ABM54]
+
+problem = ODEProblemLibrary.prob_ode_linear
+
+function test_alg(ALG; kwargs...)
+   sol_thread = solve(problem, ALG(Static.True()); kwargs...)
+   sol_nothread = solve(problem, ALG(Static.False()); kwargs...)
+
+   @test all(sol_nothread .== sol_thread)
+end
+
+
+@testset "Regression test for threading versions vs non threading versions" begin
+    @testset "$ALG" for ALG in fixed_timestep_algorithms
+        test_alg(ALG, dt=1//2^9)
+    end
+    @testset "$ALG" for ALG in free_timestep_algorithms
+        test_alg(ALG)
+    end
+
+end
diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl
index a0518e6f70..75b03c8ca4 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl
@@ -2,3 +2,4 @@ using SafeTestsets
 
 @time @safetestset "ABM Convergence Tests" include("abm_convergence_tests.jl")
 @time @safetestset "Adams Variable Coefficients Tests" include("adams_tests.jl")
+@time @safetestset "Regression test for threading versions vs non threading versions"  include("regression_test_threading.jl")

From 69b7d167b7f9f82bb460dc683d03df0c2cd977bd Mon Sep 17 00:00:00 2001
From: Michele Mesiti <michele.mesiti@kit.edu>
Date: Wed, 18 Dec 2024 12:04:29 +0100
Subject: [PATCH 3/3] Format ABM subdirectory

Run the reformatting command
suggested in CONTRIBUTING.md,
but keeping only the changes on the ABM directory
(which my current PR is affecting)
---
 .../src/adams_bashforth_moulton_perform_step.jl    |  2 +-
 .../test/regression_test_threading.jl              | 14 ++++++--------
 .../test/runtests.jl                               |  2 +-
 3 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl
index 5e160b2b02..e22e68c517 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_perform_step.jl
@@ -303,7 +303,7 @@ end
         t4 .= k4
         perform_step!(integrator,
             AB4Cache(u, uprev, fsalfirst, t2, t3, t4, ralk2, k, tmp, t5, t6, t7,
-                cnt,thread))
+                cnt, thread))
         k = integrator.fsallast
         @.. broadcast=false thread=thread u=uprev +
                                             (dt / 24) * (9 * k + 19 * k1 - 5 * k2 + k3)
diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl
index bfc0954c57..2a2da48b9c 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/regression_test_threading.jl
@@ -2,26 +2,24 @@ using OrdinaryDiffEqAdamsBashforthMoulton, ODEProblemLibrary
 using Test
 using Static
 
-free_timestep_algorithms = [ VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5, VCABM]
+free_timestep_algorithms = [VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5, VCABM]
 
-fixed_timestep_algorithms = [AB3,AB4, AB5, ABM32, ABM43, ABM54]
+fixed_timestep_algorithms = [AB3, AB4, AB5, ABM32, ABM43, ABM54]
 
 problem = ODEProblemLibrary.prob_ode_linear
 
 function test_alg(ALG; kwargs...)
-   sol_thread = solve(problem, ALG(Static.True()); kwargs...)
-   sol_nothread = solve(problem, ALG(Static.False()); kwargs...)
+    sol_thread = solve(problem, ALG(Static.True()); kwargs...)
+    sol_nothread = solve(problem, ALG(Static.False()); kwargs...)
 
-   @test all(sol_nothread .== sol_thread)
+    @test all(sol_nothread .== sol_thread)
 end
 
-
 @testset "Regression test for threading versions vs non threading versions" begin
     @testset "$ALG" for ALG in fixed_timestep_algorithms
-        test_alg(ALG, dt=1//2^9)
+        test_alg(ALG, dt = 1 // 2^9)
     end
     @testset "$ALG" for ALG in free_timestep_algorithms
         test_alg(ALG)
     end
-
 end
diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl
index 75b03c8ca4..fd3e5ee9a4 100644
--- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl
+++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/runtests.jl
@@ -2,4 +2,4 @@ using SafeTestsets
 
 @time @safetestset "ABM Convergence Tests" include("abm_convergence_tests.jl")
 @time @safetestset "Adams Variable Coefficients Tests" include("adams_tests.jl")
-@time @safetestset "Regression test for threading versions vs non threading versions"  include("regression_test_threading.jl")
+@time @safetestset "Regression test for threading versions vs non threading versions" include("regression_test_threading.jl")