Skip to content

Commit

Permalink
fixing more
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Jul 17, 2018
1 parent f096cef commit dc0114c
Show file tree
Hide file tree
Showing 12 changed files with 79 additions and 79 deletions.
8 changes: 4 additions & 4 deletions src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -212,15 +212,15 @@ end
@inbounds i = searchsortedfirst(ts,tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
@inbounds if ts[i] == tval
if idxs == nothing
copy!(out,timeseries[i])
copyto!(out,timeseries[i])
else
copy!(out,timeseries[i][idxs])
copyto!(out,timeseries[i][idxs])
end
elseif ts[i-1] == tval # Can happen if it's the first value!
if idxs == nothing
copy!(out,timeseries[i-1])
copyto!(out,timeseries[i-1])
else
copy!(out,timeseries[i-1][idxs])
copyto!(out,timeseries[i-1][idxs])
end
else
dt = ts[i] - ts[i-1]
Expand Down
4 changes: 2 additions & 2 deletions src/integrators/integrator_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -266,10 +266,10 @@ end
@inline function update_noise!(integrator,scaling_factor=integrator.sqdt)
if isinplace(integrator.noise)
integrator.noise(integrator.ΔW,integrator)
scale!(integrator.ΔW,scaling_factor)
rmul!(integrator.ΔW,scaling_factor)
if alg_needs_extra_process(integrator.alg)
integrator.noise(integrator.ΔZ,integrator)
scale!(integrator.ΔZ,scaling_factor)
rmul!(integrator.ΔZ,scaling_factor)
end
else
if typeof(integrator.u) <: AbstractArray
Expand Down
26 changes: 13 additions & 13 deletions src/perform_step/iif.jl
Original file line number Diff line number Diff line change
Expand Up @@ -94,9 +94,9 @@ end
integrator.g(rtmp2,uprev,p,t)

if is_diagonal_noise(integrator.sol.prob)
scale!(rtmp2,W.dW) # rtmp2 === rtmp3
rmul!(rtmp2,W.dW) # rtmp2 === rtmp3
else
A_mul_B!(rtmp3,rtmp2,W.dW)
mul!(rtmp3,rtmp2,W.dW)
end

rtmp3 .+= uprev
Expand All @@ -108,7 +108,7 @@ end

A = integrator.f[1](rtmp1,uprev,p,t)
M = expm(A*dt)
A_mul_B!(tmp,M,rtmp3)
mul!(tmp,M,rtmp3)

if integrator.iter > 1 && !integrator.u_modified
current_extrapolant!(uhold,t+dt,integrator)
Expand All @@ -120,7 +120,7 @@ end
rhs.sizeu = size(u)
nlres = alg.nlsolve(nl_rhs,uhold)

copy!(uhold,nlres)
copyto!(uhold,nlres)


end
Expand All @@ -141,15 +141,15 @@ end
integrator.g(rtmp2,uprev,p,t)
if typeof(cache) <: Union{IIF1MCache,IIF2MCache}
if is_diagonal_noise(integrator.sol.prob)
scale!(rtmp2,W.dW) # rtmp2 === rtmp3
rmul!(rtmp2,W.dW) # rtmp2 === rtmp3
else
A_mul_B!(rtmp3,rtmp2,W.dW)
mul!(rtmp3,rtmp2,W.dW)
end
else #Milstein correction
rtmp2 = M*rtmp2 # A_mul_B!(rtmp2,M,gtmp)
rtmp2 = M*rtmp2 # mul!(rtmp2,M,gtmp)
@unpack gtmp,gtmp2 = cache
#error("Milstein correction does not work.")
A_mul_B!(rtmp3,rtmp2,W.dW)
mul!(rtmp3,rtmp2,W.dW)
I = zeros(length(dW),length(dW));
Dg = zeros(length(dW),length(dW)); mil_correction = zeros(length(dW))
mil_correction .= 0.0
Expand All @@ -160,7 +160,7 @@ end
for j = 1:length(uprev)
#Kj = uprev .+ dt.*du1 + sqdt*rtmp2[:,j] # This works too
Kj = uprev .+ sqdt*rtmp2[:,j]
g(gtmp,Kj,p,t); A_mul_B!(gtmp2,M,gtmp)
g(gtmp,Kj,p,t); mul!(gtmp2,M,gtmp)
Dgj = (gtmp2 - rtmp2)/sqdt
mil_correction .+= Dgj*I[:,j]
end
Expand All @@ -170,12 +170,12 @@ end
if typeof(cache) <: IIF2MCache
integrator.f[2](t,uprev,rtmp1)
@. rtmp1 = @muladd 0.5dt*rtmp1 + uprev + rtmp3
A_mul_B!(tmp,M,rtmp1)
mul!(tmp,M,rtmp1)
elseif !(typeof(cache) <: IIF1MilCache)
@. rtmp1 = uprev + rtmp3
A_mul_B!(tmp,M,rtmp1)
mul!(tmp,M,rtmp1)
else
A_mul_B!(tmp,M,uprev)
mul!(tmp,M,uprev)
tmp .+= rtmp3
end

Expand All @@ -189,5 +189,5 @@ end
rhs.sizeu = size(u)
nlres = alg.nlsolve(nl_rhs,uhold)

copy!(uhold,nlres)
copyto!(uhold,nlres)
end
20 changes: 10 additions & 10 deletions src/perform_step/implicit_split_step.jl
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ end
elseif alg.extrapolant == :linear
@. u = uprev + integrator.fsalfirst*dt
else # :constant
copy!(u,uprev)
copyto!(u,uprev)
end

new_W = calc_W!(integrator, cache, dt, repeat_step)
Expand All @@ -157,15 +157,15 @@ end
end
iter += 1
f(k,u,p,t+a)
scale!(k,dt)
rmul!(k,dt)
if mass_matrix == I
k .-= z
else
A_mul_B!(vec(du1),mass_matrix,vec(z))
mul!(vec(du1),mass_matrix,vec(z))
k .-= du1
end
if has_invW(f)
A_mul_B!(vec(dz),W,vec(k)) # Here W is actually invW
mul!(vec(dz),W,vec(k)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(k),new_W)
end
Expand All @@ -188,15 +188,15 @@ end
@. u = uprev + dt*(1-theta)*tmp + theta*z
end
f(k,u,p,t+a)
scale!(k,dt)
rmul!(k,dt)
if mass_matrix == I
k .-= z
else
A_mul_B!(du1,mass_matrix,z)
mul!(du1,mass_matrix,z)
k .-= du1
end
if has_invW(f)
A_mul_B!(dz,W,k) # Here W is actually invW
mul!(dz,W,k) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(k),false)
end
Expand Down Expand Up @@ -236,7 +236,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. gtmp2 = gtmp*dW
else
A_mul_B!(gtmp2,gtmp,dW)
mul!(gtmp2,gtmp,dW)
end

if typeof(cache) <: ISSEulerHeunCache
Expand All @@ -247,7 +247,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. gtmp2 = gtmp*dW
else
A_mul_B!(gtmp2,gtmp,dW)
mul!(gtmp2,gtmp,dW)
end
end

Expand All @@ -262,7 +262,7 @@ end
f(Val{:jac},J,uprev,p,t)
end

A_mul_B!(vec(z),J,vec(tmp))
mul!(vec(z),J,vec(tmp))
@. k = dt*dt*z/2

# k is Ed
Expand Down
28 changes: 14 additions & 14 deletions src/perform_step/kencarp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -283,7 +283,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. z₄ = chi2*g1 # use z₄ as storage for the g1*chi2
else
A_mul_B!(z₄,g1,chi2) # use z₄ as storage for the g1*chi2
mul!(z₄,g1,chi2) # use z₄ as storage for the g1*chi2
end

@. tmp = uprev + γ*z₁ + nb021*z₄
Expand All @@ -309,7 +309,7 @@ end
f(k,u,p,tstep)
@. b = dt*k - z₂
if has_invW(f)
A_mul_B!(vec(dz),W,vec(b)) # Here W is actually invW
mul!(vec(dz),W,vec(b)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(b),new_W)
end
Expand All @@ -327,7 +327,7 @@ end
f(k,u,p,tstep)
@. b = dt*k - z₂
if has_invW(f)
A_mul_B!(vec(dz),W,vec(b)) # Here W is actually invW
mul!(vec(dz),W,vec(b)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(b),false)
end
Expand Down Expand Up @@ -375,7 +375,7 @@ end
f(k,u,p,tstep)
@. b = dt*k - z₃
if has_invW(f)
A_mul_B!(vec(dz),W,vec(b)) # Here W is actually invW
mul!(vec(dz),W,vec(b)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(b),false)
end
Expand All @@ -393,7 +393,7 @@ end
f(k,u,p,tstep)
@. b = dt*k - z₃
if has_invW(f)
A_mul_B!(vec(dz),W,vec(b)) # Here W is actually invW
mul!(vec(dz),W,vec(b)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(b),false)
end
Expand Down Expand Up @@ -443,7 +443,7 @@ end
f(k,u,p,tstep)
@. b = dt*k - z₄
if has_invW(f)
A_mul_B!(vec(dz),W,vec(b)) # Here W is actually invW
mul!(vec(dz),W,vec(b)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(b),false)
end
Expand All @@ -461,7 +461,7 @@ end
f(k,u,p,tstep)
@. b = dt*k - z₄
if has_invW(f)
A_mul_B!(vec(dz),W,vec(b)) # Here W is actually invW
mul!(vec(dz),W,vec(b)) # Here W is actually invW
else
cache.linsolve(vec(dz),W,vec(b),false)
end
Expand Down Expand Up @@ -496,8 +496,8 @@ end
end
else
g1 .-= g4
A_mul_B!(E₂,g1,chi2)
A_mul_B!(tmp,g4,integrator.W.dW)
mul!(E₂,g1,chi2)
mul!(tmp,g4,integrator.W.dW)
for i in eachindex(u)
@inbounds u[i] = uprev[i] + a41*z₁[i] + a42*z₂[i] + a43*z₃[i] + γ*z₄[i] + eb1*k1[i] + eb2*k2[i] + eb3*k3[i] + eb4*k4[i] + tmp[i] + E₂[i]
end
Expand All @@ -511,8 +511,8 @@ end
end
else
g1 .-= g4
A_mul_B!(E₂,g1,chi2)
A_mul_B!(tmp,g4,integrator.W.dW)
mul!(E₂,g1,chi2)
mul!(tmp,g4,integrator.W.dW)
for i in eachindex(u)
@inbounds u[i] = uprev[i] + a41*z₁[i] + a42*z₂[i] + a43*z₃[i] + γ*z₄[i] + tmp[i] + E₂[i]
end
Expand All @@ -534,7 +534,7 @@ end
@inbounds dz[i] = btilde1*z₁[i] + btilde2*z₂[i] + btilde3*z₃[i] + btilde4*z₄[i] + ebtilde1*k1[i] + ebtilde2*k2[i] + ebtilde3*k3[i] + ebtilde4*k4[i] + chi2[i]*(g1[i]-g4[i])
end
else
# dz already holds A_mul_B!(dz,g1,chi2)!
# dz already holds mul!(dz,g1,chi2)!
#@. dz += btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄ + ebtilde1*k1 + ebtilde2*k2 + ebtilde3*k3 + ebtilde4*k4
for i in eachindex(dz)
@inbounds dz[i] += btilde1*z₁[i] + btilde2*z₂[i] + btilde3*z₃[i] + btilde4*z₄[i] + ebtilde1*k1[i] + ebtilde2*k2[i] + ebtilde3*k3[i] + ebtilde4*k4[i]
Expand All @@ -547,14 +547,14 @@ end
@inbounds dz[i] = btilde1*z₁[i] + btilde2*z₂[i] + btilde3*z₃[i] + btilde4*z₄[i] + chi2[i]*(g1[i]-g4[i])
end
else
# dz already holds A_mul_B!(dz,g1,chi2)!
# dz already holds mul!(dz,g1,chi2)!
@. dz += btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄
end
end
if alg.smooth_est # From Shampine
if has_invW(f)
A_mul_B!(vec(tmp),W,vec(dz))
mul!(vec(tmp),W,vec(dz))
else
cache.linsolve(vec(tmp),W,vec(dz),false)
end
Expand Down
6 changes: 3 additions & 3 deletions src/perform_step/lamba.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. tmp=L*W.dW
else
A_mul_B!(tmp,L,W.dW)
mul!(tmp,L,W.dW)
end

@. u = K+tmp
Expand Down Expand Up @@ -127,7 +127,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. tmp=L*W.dW
else
A_mul_B!(tmp,L,W.dW)
mul!(tmp,L,W.dW)
end

@. tmp = K+tmp
Expand All @@ -139,7 +139,7 @@ end
@. tmp=(1/2)*W.dW*(L+gtmp)
else
@. gtmp = (1/2)*(L+gtmp)
A_mul_B!(tmp,gtmp,W.dW)
mul!(tmp,gtmp,W.dW)
end

dto2 = dt*(1/2)
Expand Down
10 changes: 5 additions & 5 deletions src/perform_step/low_order.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. rtmp2 *= W.dW # rtmp2 === rtmp3
else
A_mul_B!(rtmp3,rtmp2,W.dW)
mul!(rtmp3,rtmp2,W.dW)
end

@. u += rtmp3
Expand Down Expand Up @@ -72,7 +72,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. nrtmp=gtmp1*W.dW
else
A_mul_B!(nrtmp,gtmp1,W.dW)
mul!(nrtmp,gtmp1,W.dW)
end

@. tmp = @muladd uprev + ftmp1*dt + nrtmp
Expand All @@ -84,7 +84,7 @@ end
@. nrtmp=(1/2)*W.dW*(gtmp1+gtmp2)
else
@. gtmp1 = (1/2)*(gtmp1+gtmp2)
A_mul_B!(nrtmp,gtmp1,W.dW)
mul!(nrtmp,gtmp1,W.dW)
end

dto2 = dt*(1/2)
Expand Down Expand Up @@ -210,10 +210,10 @@ end
if integrator.opts.adaptive
ggprime_norm += integrator.opts.internalnorm(Dgj)
end
A_mul_B!(tmp,Dgj,@view(I[:,j]))
mul!(tmp,Dgj,@view(I[:,j]))
mil_correction .+= tmp
end
A_mul_B!(tmp,L,dW)
mul!(tmp,L,dW)
@. u .= uprev + dt*du1 + tmp + mil_correction

if integrator.opts.adaptive
Expand Down
4 changes: 2 additions & 2 deletions src/perform_step/predcorr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. gtmp *= dW
else
A_mul_B!(gdWtmp,gtmp,dW)
mul!(gdWtmp,gtmp,dW)
end

@. utmp = uprev + ftmp*dt + gdWtmp
Expand All @@ -53,7 +53,7 @@ end
if is_diagonal_noise(integrator.sol.prob)
@. gtmp *= dW
else
A_mul_B!(gdWtmp,gtmp,dW)
mul!(gdWtmp,gtmp,dW)
end

@. u += theta*(ftmp - eta*bbprimetmp)*dt + eta*gdWtmp
Expand Down
Loading

0 comments on commit dc0114c

Please sign in to comment.