Skip to content

Commit e7f7271

Browse files
ChrisRackauckasChrisRackauckas
committed
Merge branch 'broad_methods' into broadcast_changes
2 parents d903711 + 495d91a commit e7f7271

24 files changed

+913
-1649
lines changed

‎REQUIRE

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
julia 0.6.0-pre
2-
DiffEqBase 1.5.0
2+
DiffEqBase 1.10.0
33
Parameters 0.5.0
44
ForwardDiff 0.5.0
55
GenericSVD 0.0.2

‎src/caches.jl

Lines changed: 22 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -925,7 +925,7 @@ end
925925
alg_cache(alg::Feagin14,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Feagin14ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
926926

927927

928-
type Rosenbrock23Cache{uType,uArrayType,rateType,du2Type,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
928+
type Rosenbrock23Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
929929
u::uType
930930
uprev::uType
931931
k₁::rateType
@@ -946,7 +946,8 @@ type Rosenbrock23Cache{uType,uArrayType,rateType,du2Type,vecuType,JType,TabType,
946946
tab::TabType
947947
tf::TFType
948948
uf::UFType
949-
linsolve_tmp::vecuType
949+
linsolve_tmp::LinuType
950+
linsolve_tmp_vec::vecuType
950951
linsolve::F
951952
jac_config::JCType
952953
end
@@ -956,7 +957,7 @@ du_cache(c::Rosenbrock23Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsal
956957
jac_cache(c::Rosenbrock23Cache) = (c.J,c.W)
957958
vecu_cache(c::Rosenbrock23Cache) = (c.vectmp,c.vectmp2,c.vectmp3)
958959

959-
type Rosenbrock32Cache{uType,uArrayType,rateType,du2Type,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
960+
type Rosenbrock32Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
960961
u::uType
961962
uprev::uType
962963
k₁::rateType
@@ -977,7 +978,8 @@ type Rosenbrock32Cache{uType,uArrayType,rateType,du2Type,vecuType,JType,TabType,
977978
tab::TabType
978979
tf::TFType
979980
uf::UFType
980-
linsolve_tmp::vecuType
981+
linsolve_tmp::LinuType
982+
linsolve_tmp_vec::vecuType
981983
linsolve::F
982984
jac_config::JCType
983985
end
@@ -1009,11 +1011,12 @@ function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,tTypeNoUnit
10091011
vfr = VectorFReturn(f,size(u))
10101012
tf = TimeGradientWrapper(vf,uprev)
10111013
uf = UJacobianWrapper(vfr,t)
1012-
linsolve_tmp = vec(similar(u,indices(u)))
1014+
linsolve_tmp = similar(u,indices(u))
1015+
linsolve_tmp_vec = vec(linsolve_tmp)
10131016
jac_config = ForwardDiff.JacobianConfig(uf,vec(du1),vec(uprev),ForwardDiff.Chunk{determine_chunksize(u,alg)}())
10141017
Rosenbrock23Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,
1015-
fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,alg.linsolve,
1016-
jac_config)
1018+
fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,linsolve_tmp_vec,
1019+
alg.linsolve,jac_config)
10171020
end
10181021

10191022
function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
@@ -1037,9 +1040,10 @@ function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,tTypeNoUnit
10371040
vfr = VectorFReturn(f,size(u))
10381041
tf = TimeGradientWrapper(vf,uprev)
10391042
uf = UJacobianWrapper(vfr,t)
1040-
linsolve_tmp = vec(similar(u,indices(u)))
1043+
linsolve_tmp = similar(u,indices(u))
1044+
linsolve_tmp_vec = vec(linsolve_tmp)
10411045
jac_config = ForwardDiff.JacobianConfig(uf,vec(du1),vec(uprev),ForwardDiff.Chunk{determine_chunksize(u,alg)}())
1042-
Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,alg.linsolve,jac_config)
1046+
Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,linsolve_tmp_vec,alg.linsolve,jac_config)
10431047
end
10441048

10451049
immutable Rosenbrock23ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
@@ -1100,15 +1104,16 @@ vecu_cache(c::ImplicitEulerCache) = (c.uhold,)
11001104
dual_cache(c::ImplicitEulerCache) = (c.dual_cache,)
11011105

11021106
function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
1103-
u_old = similar(u,indices(u)); k = zeros(rate_prototype)
1107+
tmp = similar(u)
1108+
u_old = vec(tmp); k = zeros(rate_prototype)
11041109
dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
11051110
uhold = vec(u) # this makes uhold the same values as integrator.u
11061111
rhs = RHS_IE(f,u_old,t,t,dual_cache,size(u),eachindex(u))
11071112
fsalfirst = zeros(rate_prototype)
11081113
nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
1109-
tmp = u_old
1110-
ImplicitEulerCache{typeof(u),typeof(u_old),typeof(uhold),typeof(dual_cache),typeof(k),
1111-
typeof(rhs),typeof(nl_rhs)}(
1114+
1115+
ImplicitEulerCache{typeof(u),typeof(u_old),typeof(uhold),typeof(dual_cache),
1116+
typeof(k),typeof(rhs),typeof(nl_rhs)}(
11121117
u,uprev,uprev2,uhold,dual_cache,u_old,tmp,k,fsalfirst,rhs,nl_rhs)
11131118
end
11141119

@@ -1147,12 +1152,13 @@ vecu_cache(c::TrapezoidCache) = (c.uhold,)
11471152
dual_cache(c::TrapezoidCache) = (c.dual_cache,)
11481153

11491154
function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
1150-
u_old = similar(u,indices(u)); k = zeros(rate_prototype)
1155+
tmp = similar(u)
1156+
u_old = vec(tmp); k = zeros(rate_prototype)
11511157
uhold = vec(u); fsalfirst = zeros(rate_prototype)
1158+
f_old = vec(fsalfirst)
11521159
dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
1153-
rhs = RHS_Trap(f,u_old,fsalfirst,t,t,size(u),dual_cache,eachindex(u))
1160+
rhs = RHS_Trap(f,u_old,f_old,t,t,size(u),dual_cache,eachindex(u))
11541161
nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
1155-
tmp = u_old
11561162
TrapezoidCache{typeof(u),typeof(u_old),typeof(uhold),typeof(dual_cache),typeof(k),
11571163
typeof(rhs),typeof(nl_rhs)}(u,uprev,uprev2,uhold,u_old,fsalfirst,dual_cache,tmp,k,rhs,nl_rhs)
11581164
end
@@ -1206,7 +1212,6 @@ function alg_cache(alg::IIF1,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,
12061212
end
12071213

12081214
function alg_cache(alg::IIF1,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
1209-
12101215
tmp = similar(u,indices(u)); rtmp1 = zeros(rate_prototype)
12111216
dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
12121217
uhold = vec(u) # this makes uhold the same values as integrator.u

‎src/dense/generic_dense.jl

Lines changed: 8 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -105,18 +105,12 @@ times ts (sorted), with values timeseries and derivatives ks
105105
i = 2 # Start the search thinking it's between ts[1] and ts[2]
106106
tvals[idx[end]] > ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
107107
tvals[idx[1]] < ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
108-
if idxs == nothing
109-
if (eltype(timeseries) <: AbstractArray) && !(eltype(timeseries) <: Union{StaticArray,Array})
110-
vals = Vector{Vector{eltype(first(timeseries))}}(length(tvals))
111-
else
112-
vals = Vector{eltype(timeseries)}(length(tvals))
113-
end
114-
elseif typeof(idxs) <: Number
108+
if typeof(idxs) <: Number
115109
vals = Vector{eltype(first(timeseries))}(length(tvals))
116-
elseif eltype(timeseries) <: ArrayPartition
117-
vals = Vector{eltype(timeseries)}(length(tvals))
110+
elseif typeof(idxs) <: AbstractVector
111+
vals = Vector{Vector{eltype(first(timeseries))}}(length(tvals))
118112
else
119-
vals = Vector{Vector{eltype(first(timeseries))}}(length(tvals))
113+
vals = Vector{eltype(timeseries)}(length(tvals))
120114
end
121115
@inbounds for j in idx
122116
t = tvals[j]
@@ -189,27 +183,27 @@ times ts (sorted), with values timeseries and derivatives ks
189183
dt = ts[notsaveat_idxs[i]] - ts[notsaveat_idxs[i-1]]
190184
Θ = (t-ts[notsaveat_idxs[i-1]])/dt
191185
if typeof(cache) <: (DiscreteCache) || typeof(cache) <: DiscreteConstantCache
192-
if eltype(timeseries) <: Union{AbstractArray,ArrayPartition}
186+
if eltype(timeseries) <: AbstractArray
193187
ode_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs,deriv)
194188
else
195189
vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],0,cache,idxs,deriv)
196190
end
197191
elseif !id.dense
198-
if eltype(timeseries) <: Union{AbstractArray,ArrayPartition}
192+
if eltype(timeseries) <: AbstractArray
199193
linear_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],idxs,deriv)
200194
else
201195
vals[j] = linear_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],idxs,deriv)
202196
end
203197
elseif typeof(cache) <: CompositeCache
204198
ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache.caches[id.alg_choice[notsaveat_idxs[i-1]]]) # update the kcurrent
205-
if eltype(timeseries) <: Union{AbstractArray,ArrayPartition}
199+
if eltype(timeseries) <: AbstractArray
206200
ode_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs,deriv)
207201
else
208202
vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache.caches[id.alg_choice[notsaveat_idxs[i-1]]],idxs,deriv)
209203
end
210204
else
211205
ode_addsteps!(ks[i],ts[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],dt,f,cache) # update the kcurrent
212-
if eltype(timeseries) <: Union{AbstractArray,ArrayPartition}
206+
if eltype(timeseries) <: AbstractArray
213207
ode_interpolant!(vals[j],Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs,deriv)
214208
else
215209
vals[j] = ode_interpolant(Θ,dt,timeseries[notsaveat_idxs[i-1]],timeseries[notsaveat_idxs[i]],ks[i],cache,idxs,deriv)

‎src/dense/high_order_rk_addsteps.jl

Lines changed: 45 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -4,33 +4,33 @@ function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP8Con
44
@unpack c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = cache
55
@unpack d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = cache
66
k1 = f(t,uprev)
7-
k2 = f(t+c2*dt,uprev+dt*(a0201*k1))
8-
k3 = f(t+c3*dt,uprev+dt*(a0301*k1+a0302*k2))
9-
k4 = f(t+c4*dt,uprev+dt*(a0401*k1 +a0403*k3))
10-
k5 = f(t+c5*dt,uprev+dt*(a0501*k1 +a0503*k3+a0504*k4))
11-
k6 = f(t+c6*dt,uprev+dt*(a0601*k1 +a0604*k4+a0605*k5))
12-
k7 = f(t+c7*dt,uprev+dt*(a0701*k1 +a0704*k4+a0705*k5+a0706*k6))
13-
k8 = f(t+c8*dt,uprev+dt*(a0801*k1 +a0804*k4+a0805*k5+a0806*k6+a0807*k7))
14-
k9 = f(t+c9*dt,uprev+dt*(a0901*k1 +a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8))
15-
k10= f(t+c10*dt,uprev+dt*(a1001*k1 +a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9))
16-
k11= f(t+c11*dt,uprev+dt*(a1101*k1 +a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10))
17-
k12= f(t+dt,uprev+dt*(a1201*k1 +a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11))
18-
kupdate= b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12
19-
update = dt*kupdate
7+
k2 = f(t+c2*dt,@.(@muladd(uprev+dt*(a0201*k1))))
8+
k3 = f(t+c3*dt,@.(@muladd(uprev+dt*(a0301*k1+a0302*k2))))
9+
k4 = f(t+c4*dt,@.(@muladd(uprev+dt*(a0401*k1 +a0403*k3))))
10+
k5 = f(t+c5*dt,@.(@muladd(uprev+dt*(a0501*k1 +a0503*k3+a0504*k4))))
11+
k6 = f(t+c6*dt,@.(@muladd(uprev+dt*(a0601*k1 +a0604*k4+a0605*k5))))
12+
k7 = f(t+c7*dt,@.(@muladd(uprev+dt*(a0701*k1 +a0704*k4+a0705*k5+a0706*k6))))
13+
k8 = f(t+c8*dt,@.(@muladd(uprev+dt*(a0801*k1 +a0804*k4+a0805*k5+a0806*k6+a0807*k7))))
14+
k9 = f(t+c9*dt,@.(@muladd(uprev+dt*(a0901*k1 +a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8))))
15+
k10= f(t+c10*dt,@.(@muladd(uprev+dt*(a1001*k1 +a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9))))
16+
k11= f(t+c11*dt,@.(@muladd(uprev+dt*(a1101*k1 +a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10))))
17+
k12= f(t+dt,@.(muladd(uprev+dt*(a1201*k1 +a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11))))
18+
kupdate= @. @muladd b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12
19+
update = @. dt*kupdate
2020
utmp = uprev + update
2121
k13 = f(t+dt,utmp)
22-
k14 = f(t+c14*dt,uprev+dt*(a1401*k1 +a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13))
23-
k15 = f(t+c15*dt,uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8 +a1511*k11+a1512*k12+a1513*k13+a1514*k14))
24-
k16 = f(t+c16*dt,uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9 +a1613*k13+a1614*k14+a1615*k15))
22+
k14 = f(t+c14*dt,@.(@muladd(uprev+dt*(a1401*k1 +a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13))))
23+
k15 = f(t+c15*dt,@.(@muladd(uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8 +a1511*k11+a1512*k12+a1513*k13+a1514*k14))))
24+
k16 = f(t+c16*dt,@.(@muladd(uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9 +a1613*k13+a1614*k14+a1615*k15))))
2525
udiff = kupdate
2626
copyat_or_push!(k,1,udiff)
2727
bspl = k1 - udiff
2828
copyat_or_push!(k,2,bspl)
29-
copyat_or_push!(k,3,udiff - k13 - bspl)
30-
copyat_or_push!(k,4,(d401*k1+d406*k6+d407*k7+d408*k8+d409*k9+d410*k10+d411*k11+d412*k12+d413*k13+d414*k14+d415*k15+d416*k16))
31-
copyat_or_push!(k,5,(d501*k1+d506*k6+d507*k7+d508*k8+d509*k9+d510*k10+d511*k11+d512*k12+d513*k13+d514*k14+d515*k15+d516*k16))
32-
copyat_or_push!(k,6,(d601*k1+d606*k6+d607*k7+d608*k8+d609*k9+d610*k10+d611*k11+d612*k12+d613*k13+d614*k14+d615*k15+d616*k16))
33-
copyat_or_push!(k,7,(d701*k1+d706*k6+d707*k7+d708*k8+d709*k9+d710*k10+d711*k11+d712*k12+d713*k13+d714*k14+d715*k15+d716*k16))
29+
copyat_or_push!(k,3,@.(udiff - k13 - bspl))
30+
copyat_or_push!(k,4,@.(@muladd((d401*k1+d406*k6+d407*k7+d408*k8+d409*k9+d410*k10+d411*k11+d412*k12+d413*k13+d414*k14+d415*k15+d416*k16))))
31+
copyat_or_push!(k,5,@.(@muladd((d501*k1+d506*k6+d507*k7+d508*k8+d509*k9+d510*k10+d511*k11+d512*k12+d513*k13+d514*k14+d515*k15+d516*k16))))
32+
copyat_or_push!(k,6,@.(@muladd((d601*k1+d606*k6+d607*k7+d608*k8+d609*k9+d610*k10+d611*k11+d612*k12+d613*k13+d614*k14+d615*k15+d616*k16))))
33+
copyat_or_push!(k,7,@.(@muladd((d701*k1+d706*k6+d707*k7+d708*k8+d709*k9+d710*k10+d711*k11+d712*k12+d713*k13+d714*k14+d715*k15+d716*k16))))
3434
end
3535
end
3636

@@ -42,78 +42,45 @@ function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP8Cac
4242
@unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,update,udiff,bspl,dense_tmp3,dense_tmp4,dense_tmp5,dense_tmp6,dense_tmp7,kupdate,utilde,tmp,atmp,atmp2 = cache
4343
utmp = utilde
4444
k = [cache.udiff,cache.bspl,cache.dense_tmp3,cache.dense_tmp4,cache.dense_tmp5,cache.dense_tmp6,cache.dense_tmp7]
45-
uidx = eachindex(u)
4645
f(t,uprev,k1)
47-
for i in uidx
48-
tmp[i] = uprev[i]+dt*(a0201*k1[i])
49-
end
46+
@. tmp = @muladd uprev+dt*(a0201*k1)
5047
f(t+c2*dt,tmp,k2)
51-
for i in uidx
52-
tmp[i] = uprev[i]+dt*(a0301*k1[i]+a0302*k2[i])
53-
end
48+
@. tmp = @muladd uprev+dt*(a0301*k1+a0302*k2)
5449
f(t+c3*dt,tmp,k3)
55-
for i in uidx
56-
tmp[i] = uprev[i]+dt*(a0401*k1[i]+a0403*k3[i])
57-
end
50+
@. tmp = @muladd uprev+dt*(a0401*k1+a0403*k3)
5851
f(t+c4*dt,tmp,k4)
59-
for i in uidx
60-
tmp[i] = uprev[i]+dt*(a0501*k1[i]+a0503*k3[i]+a0504*k4[i])
61-
end
52+
@. tmp = @muladd uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
6253
f(t+c5*dt,tmp,k5)
63-
for i in uidx
64-
tmp[i] = uprev[i]+dt*(a0601*k1[i]+a0604*k4[i]+a0605*k5[i])
65-
end
54+
@. tmp = @muladd uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
6655
f(t+c6*dt,tmp,k6)
67-
for i in uidx
68-
tmp[i] = uprev[i]+dt*(a0701*k1[i]+a0704*k4[i]+a0705*k5[i]+a0706*k6[i])
69-
end
56+
@. tmp = @muladd uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
7057
f(t+c7*dt,tmp,k7)
71-
for i in uidx
72-
tmp[i] = uprev[i]+dt*(a0801*k1[i]+a0804*k4[i]+a0805*k5[i]+a0806*k6[i]+a0807*k7[i])
73-
end
58+
@. tmp = @muladd uprev+dt*(a0801*k1+a0804*k4+a0805*k5+a0806*k6+a0807*k7)
7459
f(t+c8*dt,tmp,k8)
75-
for i in uidx
76-
tmp[i] = uprev[i]+dt*(a0901*k1[i]+a0904*k4[i]+a0905*k5[i]+a0906*k6[i]+a0907*k7[i]+a0908*k8[i])
77-
end
60+
@. tmp = @muladd uprev+dt*(a0901*k1+a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8)
7861
f(t+c9*dt,tmp,k9)
79-
for i in uidx
80-
tmp[i] = uprev[i]+dt*(a1001*k1[i]+a1004*k4[i]+a1005*k5[i]+a1006*k6[i]+a1007*k7[i]+a1008*k8[i]+a1009*k9[i])
81-
end
62+
@. tmp = @muladd uprev+dt*(a1001*k1+a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
8263
f(t+c10*dt,tmp,k10)
83-
for i in uidx
84-
tmp[i] = uprev[i]+dt*(a1101*k1[i]+a1104*k4[i]+a1105*k5[i]+a1106*k6[i]+a1107*k7[i]+a1108*k8[i]+a1109*k9[i]+a1110*k10[i])
85-
end
64+
@. tmp = @muladd uprev+dt*(a1101*k1+a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
8665
f(t+c11*dt,tmp,k11)
87-
for i in uidx
88-
tmp[i] = uprev[i]+dt*(a1201*k1[i]+a1204*k4[i]+a1205*k5[i]+a1206*k6[i]+a1207*k7[i]+a1208*k8[i]+a1209*k9[i]+a1210*k10[i]+a1211*k11[i])
89-
end
66+
@. tmp = @muladd uprev+dt*(a1201*k1+a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
9067
f(t+dt,tmp,k12)
91-
for i in uidx
92-
kupdate[i] = b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k11[i]+b12*k12[i]
93-
update[i] = dt*kupdate[i]
94-
utmp[i] = uprev[i] + update[i]
95-
end
68+
@. kupdate = @muladd b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12
69+
@. update = dt*kupdate
70+
@. utmp = uprev + update
9671
f(t+dt,utmp,k13)
97-
for i in uidx
98-
tmp[i] = uprev[i]+dt*(a1401*k1[i]+a1407*k7[i]+a1408*k8[i]+a1409*k9[i]+a1410*k10[i]+a1411*k11[i]+a1412*k12[i]+a1413*k13[i])
99-
end
72+
@. tmp = @muladd uprev+dt*(a1401*k1+a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13)
10073
f(t+c14*dt,tmp,k14)
101-
for i in uidx
102-
tmp[i] = uprev[i]+dt*(a1501*k1[i]+a1506*k6[i]+a1507*k7[i]+a1508*k8[i]+a1511*k11[i]+a1512*k12[i]+a1513*k13[i]+a1514*k14[i])
103-
end
74+
@. tmp = @muladd uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8+a1511*k11+a1512*k12+a1513*k13+a1514*k14)
10475
f(t+c15*dt,tmp,k15)
105-
for i in uidx
106-
tmp[i] = uprev[i]+dt*(a1601*k1[i]+a1606*k6[i]+a1607*k7[i]+a1608*k8[i]+a1609*k9[i]+a1613*k13[i]+a1614*k14[i]+a1615*k15[i])
107-
end
76+
@. tmp = @muladd uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9+a1613*k13+a1614*k14+a1615*k15)
10877
f(t+c16*dt,tmp,k16)
109-
for i in uidx
110-
udiff[i]= kupdate[i]
111-
bspl[i] = k1[i] - udiff[i]
112-
k[3][i] = udiff[i] - k13[i] - bspl[i]
113-
k[4][i] = (d401*k1[i]+d406*k6[i]+d407*k7[i]+d408*k8[i]+d409*k9[i]+d410*k10[i]+d411*k11[i]+d412*k12[i]+d413*k13[i]+d414*k14[i]+d415*k15[i]+d416*k16[i])
114-
k[5][i] = (d501*k1[i]+d506*k6[i]+d507*k7[i]+d508*k8[i]+d509*k9[i]+d510*k10[i]+d511*k11[i]+d512*k12[i]+d513*k13[i]+d514*k14[i]+d515*k15[i]+d516*k16[i])
115-
k[6][i] = (d601*k1[i]+d606*k6[i]+d607*k7[i]+d608*k8[i]+d609*k9[i]+d610*k10[i]+d611*k11[i]+d612*k12[i]+d613*k13[i]+d614*k14[i]+d615*k15[i]+d616*k16[i])
116-
k[7][i] = (d701*k1[i]+d706*k6[i]+d707*k7[i]+d708*k8[i]+d709*k9[i]+d710*k10[i]+d711*k11[i]+d712*k12[i]+d713*k13[i]+d714*k14[i]+d715*k15[i]+d716*k16[i])
117-
end
78+
@. udiff= kupdate
79+
@. bspl = k1 - udiff
80+
@. k[3] = udiff - k13 - bspl
81+
@. k[4] = @muladd (d401*k1+d406*k6+d407*k7+d408*k8+d409*k9+d410*k10+d411*k11+d412*k12+d413*k13+d414*k14+d415*k15+d416*k16)
82+
@. k[5] = @muladd (d501*k1+d506*k6+d507*k7+d508*k8+d509*k9+d510*k10+d511*k11+d512*k12+d513*k13+d514*k14+d515*k15+d516*k16)
83+
@. k[6] = @muladd (d601*k1+d606*k6+d607*k7+d608*k8+d609*k9+d610*k10+d611*k11+d612*k12+d613*k13+d614*k14+d615*k15+d616*k16)
84+
@. k[7] = @muladd (d701*k1+d706*k6+d707*k7+d708*k8+d709*k9+d710*k10+d711*k11+d712*k12+d713*k13+d714*k14+d715*k15+d716*k16)
11885
end
11986
end

0 commit comments

Comments
 (0)