--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/benchmarkTest.jl	Fri Jan 25 16:47:51 2019 +0100
@@ -0,0 +1,19 @@
+using BenchmarkTools
+using Profile
+using ProfileView
+m = 1000;
+n = 450;
+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)

--- a/diffOp.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/diffOp.jl	Fri Jan 25 16:47:51 2019 +0100
@@ -48,8 +48,9 @@
 end
 
 # u = L*v
-function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
-    u .= 0*u
+using UnsafeArrays
+@inline function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
+    fill!(u,0)
     h = Grid.spacings(L.grid)
 
     li = LinearIndices(L.grid.numberOfPointsPerDim)
@@ -59,23 +60,26 @@
     # For each x
     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])
-
-        u[li[i,:]] += temp
+        @inbounds indices = uview(li,i,:)
+        @inbounds apply!(L.op, temp, uview(v, indices), h[2])
+        for i ∈ eachindex(indices)
+            @inbounds u[indices[i]] = temp[i]
+        end
     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])
-
-        u[li[:,i]] += temp
+        @inbounds indices = uview(li,:,i)
+        @inbounds apply!(L.op, temp, uview(v, indices), h[1])
+        for i ∈ eachindex(indices)
+            @inbounds u[indices[i]] += temp[i]
+        end
     end
 
-    u .= L.a*u
+    for i ∈ eachindex(u)
+        @inbounds u[i] = L.a*u[i]
+    end
 
     return nothing
 end

--- a/sbpD2.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/sbpD2.jl	Fri Jan 25 16:47:51 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)*applybackwards(op.closureStencils[N-i+1], v, i)/h^2
     end
 
     return nothing
@@ -26,10 +26,10 @@
     even = 1
 end
 
-struct D2{T,N,M,K} <: ConstantStencilOperator
+struct D2{T} <: ConstantStencilOperator
     quadratureClosure::Vector{T}
-    innerStencil::Stencil{T,N}
-    closureStencils::NTuple{M, Stencil{T,K}}
+    innerStencil::Stencil{T}
+    closureStencils::Vector{Stencil{T}} # TBD: Should this be a tuple?
     eClosure::Vector{T}
     dClosure::Vector{T}
     parity::Parity
@@ -48,17 +48,17 @@
     width = length(innerStencilWeights)
     r = (-div(width,2), div(width,2))
 
-    innerStencil = Stencil(r, Tuple(innerStencilWeights))
+    innerStencil = Stencil(r, innerStencilWeights)
 
     # Create boundary stencils
     boundarySize = length(d["boundary_stencils"])
-    closureStencils = Vector{Stencil}()
+    closureStencils = Vector{typeof(innerStencil)}()
 
     for i ∈ 1:boundarySize
         stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
         width = length(stencilWeights)
         r = (1-i,width-i)
-        closureStencils = (closureStencils..., Stencil(r, Tuple(stencilWeights)))
+        push!(closureStencils,Stencil(r, stencilWeights))
     end
 
     d2 = D2(

--- a/stencil.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/stencil.jl	Fri Jan 25 16:47:51 2019 +0100
@@ -10,18 +10,26 @@
 
 # Provides index into the Stencil based on offset for the root element
 function Base.getindex(s::Stencil, i::Int)
-    # TBD: Rearrange to mark with @boundscheck?
-    if s.range[1] <= i <= s.range[2]
-        return s.weights[1 + i - s.range[1]]
-    else
-        return 0
+    @boundscheck if i < s.range[1] || s.range[2] < i
+        return eltype(s.weights)(0)
     end
+    return s.weights[1 + i - s.range[1]]
 end
 
-function apply(s::Stencil, v::AbstractVector, i::Int)
+Base.@propagate_inbounds 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 weight = s[j]
+        w += weight*v[i+j]
     end
     return w
 end
+
+Base.@propagate_inbounds function applybackwards(s::Stencil, v::AbstractVector, i::Int)
+    w = zero(eltype(v))
+    for j ∈ s.range[1]:s.range[2]
+        @inbounds weight = s[j]
+        w += weight*v[i-j]
+    end
+    return w
+end