Skip to content

Commit

Permalink
Reduce allocations in the Rosenbrock methods
Browse files Browse the repository at this point in the history
This changes uses of `eye(J) - ....` to `I - ...`. `I`, the `UniformScaling` operator, uses dispatch to perform this calculation without actually allocating the identity matrix, improving the speed in almost every case. 

Also, in `oderosenbrock` someone had `eye(J)/gamma/hs` ... that's the wrong equation. I fixed it to I/(gamma*hs). That's probably why there was a `#Fix Me!` right above it.
  • Loading branch information
ChrisRackauckas committed Nov 12, 2016
1 parent 5736f4d commit 80300ae
Showing 1 changed file with 4 additions and 8 deletions.
12 changes: 4 additions & 8 deletions src/ODE.jl
Expand Up @@ -281,14 +281,14 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
end

if size(J,1) == 1
W = one(J) - h*d*J
W = I - h*d*J
else
# note: if there is a mass matrix M on the lhs of the ODE, i.e.,
# M * dy/dt = F(t,y)
# we can simply replace eye(J) by M in the following expression
# (see Sec. 5 in [SR97])

W = lufact( eye(J) - h*d*J )
W = lufact( I - h*d*J )
end

# approximate time-derivative of F
Expand Down Expand Up @@ -363,12 +363,8 @@ function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)
hs = h[solstep]
xs = x[solstep]
dFdx = G(ts, xs)
# FIXME
if size(dFdx,1) == 1
jac = 1/gamma/hs - dFdx[1]
else
jac = eye(dFdx)/gamma/hs - dFdx
end

jac = I/(gamma*hs) - dFdx

g = Array(typeof(x0), size(a,1))
g[1] = (jac \ F(ts + b[1]*hs, xs))
Expand Down

0 comments on commit 80300ae

Please sign in to comment.