# Clamped View

**URL:** <https://discourse.julialang.org/t/clamped-view/96717>\
**Category:** Performance\
**Tags:** arrays, views\
**Created:** [March 28, 2023, 3:21pm UTC](https://discourse.julialang.org/t/clamped-view/96717 "2023-03-28T15:21:19Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2023, 3:21pm UTC](https://discourse.julialang.org/t/clamped-view/96717/1 "2023-03-28T15:21:19Z")

</div>

I am trying to define a **performant** read-only view into an array where the index is clamped to lie within the array axis. The motivation for this is to apply stencils like finite differences and convolutions without needing padding or offsets or subarrays. Here’s a quick version, but it’s much slower than using an OffsetArray. (PaddedViews has similar performance to mine.)

```julia
struct ClampedView{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    parent :: A
    @inline ClampedView(parent::AbstractArray{T,N}) where {T,N} = new{T,N,typeof(parent)}(parent)
end
Base.parent(A::ClampedView) = A.parent
Base.size(A::ClampedView) = size(parent(A))
Base.@propagate_inbounds @inline Base.getindex(A::ClampedView{<:Any,N}, I::Vararg{Int,N}) where N = 
    @inbounds parent(A)[clamp.(I,axes(parent(A)))...]

@inline δ(i) = CartesianIndex(ntuple(j -> j==i ? 1 : 0, 2))
laplace!(b,a) = @inbounds @simd for I ∈ CartesianIndices(b)
    b[I] = a[I-δ(1)]+a[I+δ(1)]+a[I-δ(2)]+a[I+δ(2)]-4a[I]
end

N = (2^10,2^10)
a = ClampedView(rand(Float32,N));
a[CartesianIndex(-1,1)] # works
b = zeros(Float32,N);
laplace!(b,a); # works

using OffsetArrays: Origin
aPadded = Origin(0)(rand(Float32,N.+2));

using BenchmarkTools
@btime laplace!(b,a) #1.163 ms (0 allocations: 0 bytes)
@btime laplace!(b,aPadded) #135.400 μs (0 allocations: 0 bytes)

```

(Edited based on @mikmoore 's suggestion)

A view of a subarray is essentially the same cost as the original array. Is it possible to get similar performance for this clamped view idea? What is the secret-sauce that allows views to be so fast?

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [March 28, 2023, 3:51pm UTC](https://discourse.julialang.org/t/clamped-view/96717/2 "2023-03-28T15:51:16Z")

</div>

> [@weymouth](#):
>
> ```julia
> I2 = CartesianIndex(0,1); I1 = CartesianIndex(1,0)
> laplace!(b,a) = @inbounds @simd for I ∈ CartesianIndices(b)
> b[I] = a[I-I1]+a[I+I1]+a[I-I2]+a[I+I2]-4a[I]
> end
> 
> ```

`I2` and `I1` are global non-constants. Move them inside the function definition.

Your use of `@inbounds` in `laplace!` is unsafe because you haven’t validated the bounds of `a`. Instead, your `ClampedView` is always technically `@inbounds` in its access. However, you haven’t quite set up its `getindex` correctly, so I don’t believe it’s actually eliding the inner bounds checks. Further, internally it can still access out-of-bounds if the parent array does not have `1:N` axes. I’m trying to find a good documentation of its use but the best I’ve found is [this bit](https://docs.julialang.org/en/v1/base/base/#Base.@boundscheck).

In any case, try something like

```julia
Base.getindex(A::ClampedView{<:Any,N}, I::Vararg{Int,N}) where N = @inbounds parent(A)[clamp.(I,axes(parent(A)))...]

```

You will never need to `@inbounds` an access to a `ClampedView` because it always checks the parent array and ensures it’s accessing the closest valid index. Meanwhile, with the above definition it is always using `@inbounds` to access the parent array at the clamped index.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [March 28, 2023, 4:10pm UTC](https://discourse.julialang.org/t/clamped-view/96717/3 "2023-03-28T16:10:14Z")

</div>

Hoping you get something simple and performant. I can remember a couple of times where I wanted something like this.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2023, 4:19pm UTC](https://discourse.julialang.org/t/clamped-view/96717/4 "2023-03-28T16:19:18Z")

</div>

Thanks for the feedback! The globals was just from writing the MWE too quickly. I’ve replaced them with inlined functions in the code above.

Your version of the clamping is certainly better and I like that it will work with general `axes`. The timing using ClampedView and OffsetArray now match when using `laplace!` without `@inbounds`. But with that macro applied (which is safe in the context I will be using the function), the OffsetArray becomes 8x faster, but the ClampedView has no speed up.

```julia
@btime laplace!(b,a) #1.163 ms (0 allocations: 0 bytes)
@btime laplace!(b,aPadded) #135.400 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [March 28, 2023, 5:10pm UTC](https://discourse.julialang.org/t/clamped-view/96717/5 "2023-03-28T17:10:19Z")

</div>

Indeed, in a sense the `ClampedView` always does bounds checks because it has to move the indices inward if they aren’t legal. The `OffsetArray` is padded out to excess size, so (in this setting) it is always accessed inbounds and no checking of any kind is necessary.

I’m surprised at the 8x performance difference you see, however. I’m suspecting this has something to do with SIMD but don’t have a great idea what’s going wrong or how to fix it. You might add a `@fastmath` to your `b[I] = ...` line to allow it to re-order the summation, but I don’t expect that to change things significantly.

At the tedious end of your options, you could replace your `CartesianIndices(b)` loop with a double-loop over the rows/columns where you handle the first and last row/column specially (outside of the loop) and can simply use a normal `Array` for `a` with inbounds access for the interior of the array.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2023, 5:39pm UTC](https://discourse.julialang.org/t/clamped-view/96717/6 "2023-03-28T17:39:19Z")

</div>

Actually, to my surprise, @simd is not helping here at all! Same timings if I take it off. `@fastmath` speeds up both loops a little, but doesn’t address the question.

So it must be the `clamp` that’s slowing us down. Either directly or because of it blocking some compiler optimizations.

But this brings me back to thinking about how views work. How is the shift in indices accomplished so well with views and can it be replicated here?

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2023, 6:41pm UTC](https://discourse.julialang.org/t/clamped-view/96717/7 "2023-03-28T18:41:01Z")

</div>

Small aside: KernelAbstractions.jl is pretty good at closing the gap between `ClampedView` and `OffsetArray`

```julia
using KernelAbstractions,CUDA,CUDA.CUDAKernels,Adapt
CUDA.allowscalar(false)

@kernel function lap_kernel(b, a)
    I = @index(Global, Cartesian)
    @fastmath @inbounds b[I] = a[I-δ(1)]+a[I+δ(1)]+a[I-δ(2)]+a[I+δ(2)]-4a[I]
end
lap_CUDA!(b,a) = lap_kernel(CUDABackend(), 64)(b, a, ndrange=size(b))
lap_CPU!(b,a) = lap_kernel(CPU(), 64)(b, a, ndrange=size(b))

@btime lap_CPU!(b,a) #215.500 μs (0 allocations: 0 bytes)
@btime lap_CPU!(b,aPadded) # 60.000 μs (0 allocations: 0 bytes) ! 3.5x faster !

b = adapt(CuArray,b);
a = adapt(CuArray,a);
aPadded = adapt(CuArray,aPadded);

@btime lap_CUDA!(b,a) #2.067 μs (0 allocations: 0 bytes)
@btime lap_CUDA!(b,aPadded) #2.178 μs (0 allocations: 0 bytes)

```

CUDA is happily burning through both at the same speed.

---

<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:** [March 28, 2023, 7:06pm UTC](https://discourse.julialang.org/t/clamped-view/96717/8 "2023-03-28T19:06:50Z")

</div>

> [@weymouth](#):
>
> The motivation for this is to apply stencils like finite differences and convolutions without needing padding or offsets or subarrays.

My experiences is that boundary checks inside inner loops are always going to be slow, and are inflexible as well. What if you want to change to periodic boundary conditions, for example? Better to use padding/ghost elements, which allows you to completely decouple the boundary conditions from the inner loops.

See also [Arrays with periodic boundaries - #4 by stevengj](https://discourse.julialang.org/t/arrays-with-periodic-boundaries/4015/4) or [Seemingly unnecessary allocation within for loop - #8 by stevengj](https://discourse.julialang.org/t/seemingly-unnecessary-allocation-within-for-loop/72556/8)

> [@weymouth](#):
>
> How is the shift in indices accomplished so well with views and can it be replicated here?

The shift in indices for views is not conditional — it’s just a constant offset. Indeed, the compiler can probably pull the offset completely out of an inner loop by simply offsetting the data pointer.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 28, 2023, 7:24pm UTC](https://discourse.julialang.org/t/clamped-view/96717/9 "2023-03-28T19:24:47Z")

</div>

> [@stevengj](#):
>
> What if you want to change to periodic boundary conditions, for example?

I get what you mean, but  
`getindex(A,I)=parent(A)[mod.(I,axes(A))...]`  
is a _gorgeous_ way to do periodic conditions. It also saves me from going through and updating ghost values every iteration, and it works regardless of the stencil size. Not for every use-case, but it would clearly have its uses.

> [@stevengj](#):
>
> Indeed, the compiler can probably pull the offset completely out of an inner loop by simply offsetting the data pointer.

I suspect you are right. Applying a shift or skip to an index doesn’t seem like much less work than a min/max or mod, so the compiler might be avoiding that math completely somehow with views.

However, I’m happy to be proven wrong… Anyone have ideas on how to help the compiler figure this out?

---

<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:** [March 28, 2023, 11:00pm UTC](https://discourse.julialang.org/t/clamped-view/96717/10 "2023-03-28T23:00:17Z")

</div>

> [@weymouth](#):
>
> Applying a shift or skip to an index doesn’t seem like much less work than a min/max or mod

[Would that it were so simple](https://youtu.be/-rDw2YBUz6A?t=119). Conditionals are expensive. Mod is expensive. And both inhibit SIMD vectorization.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 29, 2023, 11:38am UTC](https://discourse.julialang.org/t/clamped-view/96717/11 "2023-03-29T11:38:25Z")

</div>

But it’s just as fast on the GPU as OffsetArrays. Are GPUs better at this?

---

<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:** [March 29, 2023, 12:24pm UTC](https://discourse.julialang.org/t/clamped-view/96717/12 "2023-03-29T12:24:32Z")

</div>

> [@weymouth](#):
>
> But it’s just as fast on the GPU as OffsetArrays.

Maybe your GPU calculation is limited by the speed of memory access.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [March 29, 2023, 2:54pm UTC](https://discourse.julialang.org/t/clamped-view/96717/13 "2023-03-29T14:54:39Z")

</div>

> [@stevengj](#):
>
> Conditionals are expensive. Mod is expensive. And both inhibit SIMD vectorization.

Yes, branches are expensive and prohibit SIMD – but not all conditionals are branches. Most `min`/`max` operations can be computed branchlessly via `cmov`/`blend`-style instructions (if not native `min`/`max`) and do not intrinsically restrict SIMD. `mod` appears to be branchless (but still expensive, maybe complicated enough to not SIMD) for fixed divisors but is branchy for variable divisors, so would likely prohibit SIMD for compile-time unknown divisors.

However, I imagine the more fundamental issue might be that the index modifications break SIMD by producing irregular memory access patterns.

Dividing the kernel into a boundary component and an interior component would at least allow the interior a good shot at vectorization.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 29, 2023, 3:05pm UTC](https://discourse.julialang.org/t/clamped-view/96717/14 "2023-03-29T15:05:28Z")

</div>

True, and this is what I’m doing now. It’s just a lot of code which distracts from the main point and introduces potential errors.