julianized code

julia/leapfrog.jl

@@ -1,57 +1,50 @@
 using StaticArrays
+using LinearAlgebra: norm
 
-const dims = 3
-const PosVector = SVector{dims, Float64}
-const zerov = zero(SVector{dims})
-
-@inline function gravitate!(x::Vector{PosVector}, a::Vector{PosVector},
-                            m::Vector)
+function gravitate!(x::Vector{T}, a::Vector{T},m::AbstractVector) where T
     G = 6.6742367e-11
-    fill!(a, zerov)
-    @inbounds for i in 1:length(m)
-        for j in 1:length(m)
+    fill!(a, zero(T))
+    for i in eachindex(x)
+        for j in eachindex(x)
             if i == j
                 continue
             end
-
             d = x[i] - x[j]
             prefactor = -G/norm(d)^3 * m[j]
             a[i] += prefactor*d
         end
     end
 end
 
-function main(;dt::Real=0.08, tmax::Real=3.6525e8)
-
-    # Arrays (initialized to zero by default)
-    x = [zerov,
-         @SVector([0.0162, 6.57192058353e-15, 5.74968548652e-16])]
-
-    v = [zerov,
-         @SVector([-1.48427302304e-14, 0.0399408809121, 0.00349437429104])]
-
-    a = [zerov, zerov]
-
-    m = [0.08, 3.0e-6]   # M_Sun
-
+function main(
+    x::Vector{T},v::Vector{T},m::AbstractVector;
+    dt::Real=0.08, tmax::Real=3.6525e8
+) where T
+    a = similar(x)
     hdt = dt / 2
-
-    @inbounds for k in 1:round(Int, tmax/dt)
-        for i in 1:length(m)
-            x[i] += hdt*v[i]
-        end
+    for k in 1:round(Int,tmax/dt)
+        @. x += hdt * v
         gravitate!(x, a, m)
-        for i in 1:length(m)
-            v[i] += a[i]*dt
-            x[i] += v[i]*hdt
-        end
+        @. v += a*dt
+        @. x += v*hdt
     end
-
     println("Positions:", x)
 end
 
-println("Running up to tmax=10 to compile...")
-@time main(tmax=1e2)
 println("Running full benchmark...")
-@time main()
+
+# Dimension of particles
+T = SVector{3,Float64}
+
+# Initial positions and velocities
+x = [zero(T),T(0.0162, 6.57192058353e-15, 5.74968548652e-16)]
+v = [zero(T),T(-1.48427302304e-14, 0.0399408809121, 0.00349437429104)]
+
+# Masses
+m = @SVector [0.08, 3.0e-6]
+
+# Run simulation
+main(x,v,m)
+
+
 

julia/target/optimized/leapfrog

@@ -1,3 +1,3 @@
 #!/bin/bash
 # OSX: /Applications/Julia-0.5.app/Contents/Resources/julia/bin/julia
-julia -O3 leapfrog.jl
+julia -O3 --check-bounds=no leapfrog.jl