addsteps fix, and muladds

SciML · Feb 2, 2017 · 16639fb · 16639fb
1 parent 70582a1
commit 16639fb
Show file tree

Hide file tree

Showing 8 changed files with 155 additions and 129 deletions.
diff --git a/src/dense/low_order_rk_addsteps.jl b/src/dense/low_order_rk_addsteps.jl
@@ -150,7 +150,7 @@ function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Tsit5C
     end
     f(t+dt,tmp,k6)
     for i in uidx
-      u[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
+      tmp[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
     end
     f(t+dt,u,k7)
     copyat_or_push!(k,1,k1)

diff --git a/src/integrators/explicit_rk_integrator.jl b/src/integrators/explicit_rk_integrator.jl
@@ -14,28 +14,28 @@ end
   for i = 2:stages-1
     utilde = zero(kk[1])
     for j = 1:i-1
-      utilde += A[j,i]*kk[j]
+      utilde = @muladd utilde + A[j,i]*kk[j]
     end
-    kk[i] = f(t+c[i]*dt,uprev+dt*utilde);
+    kk[i] = f(@muladd(t+c[i]*dt),@muladd(uprev+dt*utilde));
   end
   #Calc Last
   utilde = zero(kk[1])
   for j = 1:stages-1
-    utilde += A[j,end]*kk[j]
+    utilde = @muladd utilde + A[j,end]*kk[j]
   end
-  kk[end] = f(t+c[end]*dt,uprev+dt*utilde); integrator.fsallast = kk[end] # Uses fsallast as temp even if not fsal
+  kk[end] = f(@muladd(t+c[end]*dt),@muladd(uprev+dt*utilde)); integrator.fsallast = kk[end] # Uses fsallast as temp even if not fsal
   # Accumulate Result
   utilde = α[1]*kk[1]
   for i = 2:stages
-    utilde += α[i]*kk[i]
+    utilde = @muladd utilde + α[i]*kk[i]
   end
-  u = uprev + dt*utilde
+  u = @muladd uprev + dt*utilde
   if integrator.opts.adaptive
     uEEst = αEEst[1]*kk[1]
     for i = 2:stages
-      uEEst += αEEst[i]*kk[i]
+      uEEst = @muladd uEEst + αEEst[i]*kk[i]
     end
-    integrator.EEst = integrator.opts.internalnorm( dt*(utilde-uEEst)/(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
+    integrator.EEst = integrator.opts.internalnorm( dt*(utilde-uEEst)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
   end
   if isfsal(integrator.alg.tableau)
     integrator.fsallast = kk[end]
@@ -69,39 +69,39 @@ end
     end
     for j = 1:i-1
       for l in uidx
-        utilde[l] += A[j,i]*kk[j][l]
+        utilde[l] = @muladd utilde[l] + A[j,i]*kk[j][l]
       end
     end
     for l in uidx
-      tmp[l] = uprev[l]+dt*utilde[l]
+      tmp[l] = @muladd uprev[l]+dt*utilde[l]
     end
-    f(t+c[i]*dt,tmp,kk[i])
+    f(@muladd(t+c[i]*dt),tmp,kk[i])
   end
   #Last
   for l in uidx
     utilde[l] = zero(kk[1][1])
   end
   for j = 1:stages-1
     for l in uidx
-      utilde[l] += A[j,end]*kk[j][l]
+      utilde[l] = @muladd utilde[l] + A[j,end]*kk[j][l]
     end
   end
   for l in uidx
-    u[l] = uprev[l]+dt*utilde[l]
+    u[l] = @muladd uprev[l]+dt*utilde[l]
   end
-  f(t+c[end]*dt,u,kk[end]) #fsallast is tmp even if not fsal
+  f(@muladd(t+c[end]*dt),u,kk[end]) #fsallast is tmp even if not fsal
   #Accumulate
   if !isfsal(integrator.alg.tableau)
     for i in uidx
       utilde[i] = α[1]*kk[1][i]
     end
     for i = 2:stages
       for l in uidx
-        utilde[l] += α[i]*kk[i][l]
+        utilde[l] = @muladd utilde[l] + α[i]*kk[i][l]
       end
     end
     for i in uidx
-      u[i] = uprev[i] + dt*utilde[i]
+      u[i] = @muladd uprev[i] + dt*utilde[i]
     end
   end
   if integrator.opts.adaptive
@@ -110,11 +110,11 @@ end
     end
     for i = 2:stages
       for j in uidx
-        uEEst[j] += αEEst[i]*kk[i][j]
+        uEEst[j] = @muladd uEEst[j] + αEEst[i]*kk[i][j]
       end
     end
     for i in uidx
-      atmp[i] = (dt*(utilde[i]-uEEst[i])/(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
+      atmp[i] = (dt*(utilde[i]-uEEst[i])/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
     end
     integrator.EEst = integrator.opts.internalnorm(atmp)
   end

diff --git a/src/integrators/feagin_rk_integrators.jl b/src/integrators/feagin_rk_integrators.jl
@@ -9,7 +9,7 @@ end
   @unpack adaptiveConst,a0100,a0200,a0201,a0300,a0302,a0400,a0402,a0403,a0500,a0503,a0504,a0600,a0603,a0604,a0605,a0700,a0704,a0705,a0706,a0800,a0805,a0806,a0807,a0900,a0905,a0906,a0907,a0908,a1000,a1005,a1006,a1007,a1008,a1009,a1100,a1105,a1106,a1107,a1108,a1109,a1110,a1200,a1203,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1300,a1302,a1303,a1305,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1400,a1401,a1404,a1406,a1412,a1413,a1500,a1502,a1514,a1600,a1601,a1602,a1604,a1605,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613,a1614,a1615,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16 = cache
   k1  = integrator.fsalfirst
   k2  = f(t + c1*dt,uprev  + dt*(a0100*k1))
-  k3  = f(t + c2*dt ,uprev + dt*(a0200*k1 + a0201*k2))
+  k3  = f(t + c2*dt,uprev + dt*(a0200*k1 + a0201*k2))
   k4  = f(t + c3*dt,uprev  + dt*(a0300*k1              + a0302*k3))
   k5  = f(t + c4*dt,uprev  + dt*(a0400*k1              + a0402*k3 + a0403*k4))
   k6  = f(t + c5*dt,uprev  + dt*(a0500*k1                           + a0503*k4 + a0504*k5))
@@ -24,8 +24,7 @@ end
   k15 = f(t + c14*dt,uprev + dt*(a1400*k1 + a1401*k2                           + a1404*k5              + a1406*k7 +                                                                     a1412*k13 + a1413*k14))
   k16 = f(t + c15*dt,uprev + dt*(a1500*k1              + a1502*k3                                                                                                                                                     + a1514*k15))
   k17 = f(t + c16*dt,uprev + dt*(a1600*k1 + a1601*k2 + a1602*k3              + a1604*k5 + a1605*k6 + a1606*k7 + a1607*k8 + a1608*k9 + a1609*k10 + a1610*k11 + a1611*k12 + a1612*k13 + a1613*k14 + a1614*k15 + a1615*k16))
-  tmp = dt*((b1*k1 + b2*k2 + b3*k3 + b5*k5) + (b7*k7 + b9*k9 + b10*k10 + b11*k11) + (b12*k12 + b13*k13 + b14*k14 + b15*k15) + (b16*k16 + b17*k17))
-  u = uprev + tmp
+  u = uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b5*k5 + b7*k7 + b9*k9 + b10*k10 + b11*k11 + b12*k12 + b13*k13 + b14*k14 + b15*k15 + b16*k16 + b17*k17)
   if integrator.opts.adaptive
     integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k16) * adaptiveConst)./(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
   end

diff --git a/src/integrators/implicit_integrators.jl b/src/integrators/implicit_integrators.jl
@@ -6,7 +6,7 @@ type RHS_IE_Scalar{F,uType,tType} <: Function
 end
 
 function (p::RHS_IE_Scalar)(u,resid)
-  resid[1] = u[1] - p.u_old[1] - p.dt*p.f(p.t+p.dt,u)[1]
+  resid[1] = u[1] - p.u_old[1] - p.dt*p.f(p.t+p.dt,u[1])[1]
 end
 
 @inline function initialize!(integrator,cache::ImplicitEulerConstantCache,f=integrator.f)
@@ -103,7 +103,7 @@ function (p::RHS_Trap)(uprev,resid)
   du1 = get_du(p.dual_cache, eltype(uprev))
   p.f(p.t+p.dt,reshape(uprev,p.sizeu),du1)
   for i in p.uidx
-    resid[i] = uprev[i] - p.u_old[i] - (p.dt/2)*(du1[i]+p.f_old[i])
+    resid[i] = @muladd uprev[i] - p.u_old[i] - (p.dt/2)*(du1[i]+p.f_old[i])
   end
 end
 
@@ -150,7 +150,7 @@ type RHS_Trap_Scalar{F,uType,rateType,tType} <: Function
 end
 
 function (p::RHS_Trap_Scalar)(uprev,resid)
-  resid[1] = uprev[1] - p.u_old[1] - (p.dt/2)*(p.f_old + p.f(p.t+p.dt,uprev)[1])
+  resid[1] = @muladd uprev[1] - p.u_old[1] - (p.dt/2)*(p.f_old + p.f(p.t+p.dt,uprev[1])[1])
 end
 
 @inline function initialize!(integrator,cache::TrapezoidConstantCache,f=integrator.f)