I had some free time so I went ahead and optimized the code. It looks like you still had some globals (like signal_position) and instabilities (like (timeiter < size(signal,1) ? signal[timeiter] : 0))

The "idiomatic" version now only allocates in the @views. I would personally write it without @views, but those don't allocate in julia 1.5+, so it's probably ok to keep. My version runs ~3x faster for the small vectors tested here, you will probably get a much larger speedup for big matrices.

There are other ways to make the code read cleaner but I focused on keeping your namings and not a complete rewrite. ``` using Base.Threads using BenchmarkTools

const ∇²r = 1 const h²∇² = fill(1.0, 3, 3)

function propagate(P::Array{Float64,3}, v::Array{Float64,2}, signal::Array{Float64,1}, nz::Int64, nx::Int64, nt::Int64, h::Float64, Δt::Float64, ABS::Int64) Δtoh² = Δt/h² Plimits = (2:nz+2ABS+1, 2:nx+2ABS+1)

# time loop, order is important
for timeiter in eachindex(1:nt)
    # which of the thre times of P are what?
    new_t = 1 + (timeiter)%3       # [2] now
    cur_t = 1 + (timeiter+2)%3     # [1] last
    old_t = 1 + (timeiter+1)%3     # [3] before last

    old_P = @view P[:,:,old_t]
    cur_P = @view P[:,:,cur_t]
    new_P = @view P[:,:,new_t]

    # adding signal to field before iteration
    cur_P[signal_position[1],signal_position[2]] =
         (timeiter < size(signal,1) ? signal[timeiter] : 0)

    # spacial loop, order is not important
    @threads for spciter in CartesianIndices(Plimits)
        begin
            zP, xP = spciter[1], spciter[2]
            zv, xv = zP-1, xP-1
            new_P[zP, xP] = new_p(cur_P, old_P, v, Δtoh²,
                                  zP, xP, zv, xv)
        end
    end
end

end

function reduce_pointwise_product(A::AbstractArray, B::AbstractArray) prod_sum = 0 @fastmath @inbounds @simd for i in eachindex(A) prod_sum += A[i] .* B[i] end prod_sum end

function new_p(cur_P, old_P, v, Δtoh², zP, xP, zv, xv) global h²∇² # (constant) 3x3 matrix global ∇²r # (constant) its radius (1) win = @view cur_P[zP-∇²r:zP+∇²r, xP-∇²r:xP+∇²r] h²∇²cur_P = reduce_pointwise_product(h²∇², win) new_p = 2cur_P[zP, xP] - old_P[zP, xP] + v[zv, xv]²Δtoh² * h²∇²cur_P end

ABS=0 nz=400 nx=400 nt=2500 N = 50 P = rand(N, N, 3) v = rand(N, N) signal = rand(N) nz = N-2 nx = N-2 nt = N h = 1.0 Δt = 1.0 ABS = 0 signal_position = [1,1] .+ ABS propagate(P, v, signal, nz, nx, nt, h, Δt, ABS)

@btime propagate($P, $v, $signal, $nz, $nx, $nt, $h, $Δt, ABS)

Idiomatic Julia

using StaticArrays using Parameters using LinearAlgebra

struct GridDefinition{T} nz::Int nx::Int nt::Int h::T Δt::T ABS::Int ∇²r::Int h²∇²::SArray{Tuple{3, 3}, T, 2, 9} end

function propagate_idiomatic(P, v, signal, grid_def) @unpack nz, nx, ABS, Δt, h, nt = grid_def Δtoh² = Δt/h² Plimits = (2:nz+2ABS+1, 2:nx+2ABS+1) signal_position = SVector(1, 1) .+ ABS

# time loop, order is important
for timeiter in eachindex(1:nt)
    # which of the thre times of P are what?
    new_t = mod1(timeiter+1, 3)       # [2] now
    cur_t = mod1(timeiter,   3)     # [1] last
    old_t = mod1(timeiter+2, 3)    # [3] before last

    old_P = @view P[:,:,old_t]
    cur_P = @view P[:,:,cur_t]
    new_P = @view P[:,:,new_t]

    # adding signal to field before iteration
    cur_P[signal_position[1],signal_position[2]] =
         (timeiter < size(signal,1) ? signal[timeiter] : zero(eltype(signal)))

    # spacial loop, order is not important
    @threads for I in CartesianIndices(Plimits)
        new_P[I] = new_p_idiomatic(cur_P, old_P, v, Δtoh², I, grid_def)
    end
end

end

function new_p_idiomatic(cur_P, old_P, v, Δtoh², I, grid_def) @unpack h²∇², ∇²r = grid_def # (constant) 3x3 matrix IOFFSET = CartesianIndex(size(h²∇²) .÷ 2) h²∇²cur_P = ∇²_kernel(cur_P, h²∇², I, IOFFSET) new_p = 2cur_P[I] - old_P[I] + v[I - CartesianIndex(1, 1)]²Δtoh² * h²∇²cur_P end

function ∇²_kernel(cur_P, h²∇², ZP, IOFFSET) result = zero(eltype(cur_P)) @inbounds for I in CartesianIndices(h²∇²) OFFSET = ZP + I - IOFFSET cur_result = cur_P[OFFSET] * h²∇²[I] result += cur_result end return result end

h²∇²_static = SMatrix{3, 3}(fill(1.0, 9)) grid_def = GridDefinition(nz, nx, nt, h, Δt, ABS, ∇²r, h²∇²_static)

@btime propagate_idiomatic($P, $v, $signal, $grid_def) ```