--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/benchmarkTest.jl	Thu Jan 24 14:58:22 2019 +0100
@@ -0,0 +1,21 @@
+include("sbp.jl")
+using BenchmarkTools
+using Profile
+using ProfileView
+
+function benchmarkTest(m,n)
+  g = sbp.Grid.EquidistantGrid((m,n), (0.0, 0.0), (2pi, 3/2*pi))
+  op = sbp.readOperator("d2_4th.txt","h_4th.txt")
+  Laplace = sbp.Laplace(g, 1.0, op)
+
+  init(x,y) = sin(x) + cos(y)
+  v = sbp.Grid.evalOn(g,init)
+
+  u = zeros(eltype(v),length(v))
+
+  sbp.apply!(Laplace,u,v)
+
+  @time sbp.apply!(Laplace,u,v)
+
+  #@benchmark sbp.apply!(Laplace,u,v)
+end

--- a/diffOp.jl	Thu Jan 24 14:33:49 2019 +0100
+++ b/diffOp.jl	Thu Jan 24 14:58:22 2019 +0100
@@ -48,7 +48,7 @@
 end
 
 # u = L*v
-function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
+@inline function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
     u .= 0*u
     h = Grid.spacings(L.grid)
 
@@ -60,19 +60,19 @@
     temp = zeros(eltype(u), n_y)
     for i ∈ 1:n_x
 
-        v_i = view(v, li[i,:])
-        apply!(L.op, temp, v_i, h[2])
+        @inbounds v_i = view(v, li[i,:])
+        @inbounds apply!(L.op, temp, v_i, h[2])
 
-        u[li[i,:]] += temp
+        @inbounds u[li[i,:]] += temp
     end
 
     # For each y
     temp = zeros(eltype(u), n_x)
     for i ∈ 1:n_y
-        v_i = view(v, li[:,i])
-        apply!(L.op, temp, v_i, h[1])
+        @inbounds v_i = view(v, li[:,i])
+        @inbounds apply!(L.op, temp, v_i, h[1])
 
-        u[li[:,i]] += temp
+        @inbounds u[li[:,i]] += temp
     end
 
     u .= L.a*u

--- a/sbpD2.jl	Thu Jan 24 14:33:49 2019 +0100
+++ b/sbpD2.jl	Thu Jan 24 14:58:22 2019 +0100
@@ -1,21 +1,21 @@
 abstract type ConstantStencilOperator end
 
-function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real)
+@inline function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real)
     N = length(v)
     cSize = closureSize(op)
 
     for i ∈ range(1; length=cSize)
-        u[i] = apply(op.closureStencils[i], v, i)/h^2
+        @inbounds u[i] = apply(op.closureStencils[i], v, i)/h^2
     end
 
     innerStart = 1 + cSize
     innerEnd = N - cSize
     for i ∈ range(innerStart, stop=innerEnd)
-        u[i] = apply(op.innerStencil, v, i)/h^2
+        @inbounds u[i] = apply(op.innerStencil, v, i)/h^2
     end
 
     for i ∈ range(innerEnd+1, length=cSize)
-        u[i] = Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
+        @inbounds u[i] = Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
     end
 
     return nothing

--- a/stencil.jl	Thu Jan 24 14:33:49 2019 +0100
+++ b/stencil.jl	Thu Jan 24 14:58:22 2019 +0100
@@ -25,10 +25,10 @@
     end
 end
 
-function apply(s::Stencil, v::AbstractVector, i::Int)
+@inline function apply(s::Stencil, v::AbstractVector, i::Int)
     w = zero(eltype(v))
     for j ∈ s.range[1]:s.range[2]
-        w += s[j]*v[i+j] # TBD: Make this without boundschecks?
+        @inbounds w += s[j]*v[i+j] # TBD: Make this without boundschecks?
     end
     return w
 end