Merge pull request #138 from baggepinnen/patch-1

reduce allocations by two_stage_method

Author: ChrisRackauckas

src/two_stage_method.jl

@@ -47,7 +47,7 @@ function construct_t2(t,tpoints)
 end
 function construct_w(t,tpoints,h,kernel)
     W = @. calckernel((kernel,),(tpoints-t)/h)/h
-    Matrix(Diagonal(W))
+    Diagonal(W)
 end
 function construct_estimated_solution_and_derivative!(data,kernel,tpoints)
   _one = oneunit(first(data))
@@ -57,17 +57,22 @@ function construct_estimated_solution_and_derivative!(data,kernel,tpoints)
   n = length(tpoints)
   h = (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16))
 
+  Wd = similar(data, n, size(data,1))
+  WT1 = similar(data, n, 2)
+  WT2 = similar(data, n, 3)
   x = map(tpoints) do _t
       T1 = construct_t1(_t,tpoints)
       T2 = construct_t2(_t,tpoints)
       W = construct_w(_t,tpoints,h,kernel)
-      e2'*inv(T2'*W*T2)T2'*W*data',e1'*inv(T1'*W*T1)*T1'*W*data'
+      mul!(Wd,W,data')
+      mul!(WT1,W,T1)
+      mul!(WT2,W,T2)
+      (e2'*((T2'*WT2)\T2'))*Wd,(e1'*((T1'*WT1)\T1'))*Wd
   end
   estimated_derivative = reduce(hcat,transpose.(first.(x)))
   estimated_solution = reduce(hcat,transpose.(last.(x)))
   estimated_derivative,estimated_solution
 end
-
 function construct_iip_cost_function(f,du,preview_est_sol,preview_est_deriv,tpoints)
   function (p)
       _du = DiffEqBase.get_tmp(du,p)