# Need fast Array operation!

**URL:** <https://discourse.julialang.org/t/need-fast-array-operation/132964>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [October 8, 2025, 4:58am UTC](https://discourse.julialang.org/t/need-fast-array-operation/132964 "2025-10-08T04:58:11Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [October 8, 2025, 7:27am UTC](https://discourse.julialang.org/t/need-fast-array-operation/132964/22 "2025-10-08T07:27:07Z")

</div>

My time 822ms is close enough to 1s that I’d chalk it up to our different machines. As far as I can tell from my naive translation, the few extra allocations in the `np.exp` line had a negligible impact on runtimes relative to the total number-crunching, and `bicubic` being on par with NumPy indicates that the secret sauce was vectorization (in the SIMD sense, not the array operation API sense).

This was a reasonable outcome. NumPy/SciPy and its dependencies implement many optimizations, some top-of-the-line, so straightforward usage of its functions on NumPy arrays is usually as fast as practical. The issues only really come up when throwing pure Python elementwise functions into the mix e.g. `numpy.vectorize`, in which case there are still Python-world compiler options like Numba and Cython.

To speak more broadly, many languages have an optimizing compiler like Julia, and although it’s essential for good performance, it isn’t the only factor. An optimizing compiler can transform code quite a bit, but there are sensible limits when user intent is ambiguous. NumPy’s optimizations, whether it’s vectorization or the various multithreading schemes, are manually implemented for supported platforms in various compiled dependencies and applied to carefully designed C loops over NumPy arrays; `@turbo` asks you to adhere to a list of restrictions on your loop, otherwise it’ll blindly run into mistakes like the out-of-bounds overwriting. I’d say that NumPy was actually more convenient in this case, and Julia’s tendency to put optimizations in our hands warrants improvements to the ecosystem and discoverability.

---

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [October 8, 2025, 12:59pm UTC](https://discourse.julialang.org/t/need-fast-array-operation/132964/23 "2025-10-08T12:59:05Z")

</div>

> function bicubic(  
> N\_layer::Integer, kernel\_x::AbstractMatrix{F}, kernel\_y::AbstractMatrix{F},

Your current Julia code is a pure-Julia implementation, and it’s already near the performance ceiling. In your specific use case, even routing the current Julia kernel through BLAS is unlikely to help.

Why BLAS won’t move the needle here:

- The kernel is a fixed 4×4 stencil (16 fused multiply-adds) per output plus an `exp10`. That granularity is too small—BLAS call overhead would dominate.
- The access pattern is “gather 4×4 patches at (ix0[n]+i, iy0[n]+j) and accumulate over (k,l),” which isn’t GEMM-shaped or contiguous. Forcing it into AXPY/GEMM would hit strided, cache-unfriendly slices.
- The real costs are scattered memory traffic and the exponential; BLAS doesn’t accelerate either. A tight vectorized microkernel (e.g., LoopVectorization) already emits near-roofline FMA code on CPU.

Given that, achieving parity with NumPy’s C ufunc path under **strict single-thread** settings on both sides is already an excellent outcome.

---

<div class="post-metadata">

**Author:** ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)\
**Post date:** [October 8, 2025, 3:15pm UTC](https://discourse.julialang.org/t/need-fast-array-operation/132964/24 "2025-10-08T15:15:58Z")

</div>

I tried to write a kernel just to see on gpu

```julia-auto
using KernelAbstractions
using BenchmarkTools

@kernel inbounds=true unsafe_indices=true function bicubic_kernel!(
    buf, N_layer, kernel_x, kernel_y, ix0, iy0, kcoeff
)
    k, l, n = @index(Global, NTuple)

    klim, llim = size(kcoeff,1), size(kcoeff,2)

    @fastmath if k <= klim && l <= llim && n <= N_layer
        ix0n = ix0[n]
        iy0n = iy0[n]

        acc = zero(eltype(buf))
        @inbounds for i in 1:4, j in 1:4
            acc += kernel_x[i,n] * kernel_y[j,n] *
                   kcoeff[k, l, ix0n + i, iy0n + j]
        end

        buf[k, l, n] = exp10(acc)
    end
end

function bicubic2( N_layer, kernel_x, kernel_y, ix0, iy0, kcoeff)
    klim, llim = size(kcoeff, 1), size(kcoeff, 2)
    buf = KernelAbstractions.allocate(get_backend(kcoeff),eltype(kcoeff), klim, llim,N_layer)
    bicubic_kernel!(get_backend(kcoeff))(buf, N_layer, kernel_x, kernel_y, ix0, iy0, kcoeff;ndrange=size(buf))
    buf
end
kcoeff_intp = zeros(Float64, 100, 5011, 10)
kernel_x = rand(4,100)
kernel_y = rand(4,100)
ix0 = rand(0:36, 100) # subtract 1 if you want 0-based like NumPy
iy0 = rand(0:16, 100)
kcoeff = rand( 5011, 10,40, 20)
N_layer = 100
@btime bicubic2($N_layer,$kernel_x,$kernel_y,$ix0,$iy0,$kcoeff);

```

and suprisingly I get the best time on cpu and gpu I got yet  
`39.512 ms (198 allocations: 38.25 MiB)` on cpu and `2.641 ms (108 allocations: 3.19 KiB)` on gpu, for single threaded (ie `groupsize = size(buf)`) I get `109.986 ms (19 allocations: 38.23 MiB)` which isn’t so bad yet maybe beatable

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 8, 2025, 11:56pm UTC](https://discourse.julialang.org/t/need-fast-array-operation/132964/25 "2025-10-08T23:56:50Z")

</div>

> [@avnsh9](#):
>
> ```julia-auto
> function my_bicubic_implementation!(cache::Array{Float64}, N_layer::Int64, kernel_x::Matrix{Float64},kernel_y::Matrix{Float64},ix0::Vector{Int64},iy0::Vector{Int64},kcoeff::Array{Float64})
>     
> 
> ```

As an aside, all of these type declarations do nothing for performance. They just make your code less generic (e.g. it doesn’t work in single precision).

---

<div class="post-metadata">

**Author:** ![avnsh9](https://avatars.discourse-cdn.com/v4/letter/a/f07891/32.png) [@avnsh9](https://discourse.julialang.org/u/avnsh9)\
**Post date:** [October 9, 2025, 10:06am UTC](https://discourse.julialang.org/t/need-fast-array-operation/132964/26 "2025-10-09T10:06:55Z")

</div>

Sorry but i was placed on hold being a new user for 24 hours 😁.There is no multithreading in my version of julia implementation, you can verify this yourself by running it and open a htop window side by side, you won’t see more then one thread being used

[Previous page](https://discourse.julialang.org/t/need-fast-array-operation/132964.md?page=1)