# Looking for elegant cartesian grid array assignment

**URL:** <https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360>\
**Category:** General Usage\
**Created:** [July 20, 2020, 5:11am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360 "2020-07-20T05:11:43Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![purplishrock](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/purplishrock/32/13451_2.png) [@purplishrock](https://discourse.julialang.org/u/purplishrock)\
**Post date:** [July 20, 2020, 5:11am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/1 "2020-07-20T05:11:43Z")

</div>

I’m building an array of cube centers for a 3-D cartesion grid, something like (ignore for the moment, for simplicity, that the cube centers are actually located at an offset of dx/2, dy/2, dz/2).

```julia
for i=1:nx
  for j=1:ny
    for k = 1:nz
        a[i,j,k,1] = i*dx # i=1:nx
        a[i,j,k,2] = j*dy # j=1:ny
        a[i,j,k,3] = k*dz # k=1:nz
    end
  end
end

```

realizing that the sequences are repeated, i.e. you can do

```julia
xs = dx * range(0,nx-1,nx)
for j=1:ny
  for k=1:nz
    a[:,j,k,1] .= xs
  end
end

```

and then you need two more loops for the 2nd and 3rd indexes.

Works great, but it’s not clever enough 😆 . It seems like this is amenable to a one-liner hack. I am seeking that one-liner hack, just for fun. Even the brute force 3-loop calculation of the values for a 64M array takes LESS than 2 seconds, so not worried about efficiency at all.

I’m just figuring some julia hacker out there can give me a valuable lesson in broadcasting 🙂

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 20, 2020, 5:16am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/2 "2020-07-20T05:16:13Z")

</div>

You are iterating the loops in an order suboptimal for the memory layout. Iterate the first dimension in the tightest loop, i.e, swap the loop order in the first snippet.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 20, 2020, 5:19am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/3 "2020-07-20T05:19:11Z")

</div>

Doesn’t

```julia
a[:,:, :,1] .= xs

```

work?

---

<div class="post-metadata">

**Author:** ![purplishrock](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/purplishrock/32/13451_2.png) [@purplishrock](https://discourse.julialang.org/u/purplishrock)\
**Post date:** [July 20, 2020, 5:47am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/4 "2020-07-20T05:47:31Z")

</div>

> [@baggepinnen](#):
>
> . Iterate the first dimension in the tightest loop, i.e, swap the loop order in the first snippet.

yowzer. a factor of 3 speed-up. I really should have known that as many discussions as i’ve seen about the column-major-ness of julia.

> [@baggepinnen](#):
>
> Doesn’t
> 
> ```julia
> a[:,:, :,1] .= xs
> 
> ```
> 
> work?

yes it does! How ?! More impressively if the first and second dimensions are the same it still works… However

```julia
ys = dy * range(0, ny-1, ny)

a1[:,:,:,2] .= ys

```

does not. But if you make ys a 1xny vector then i think it does work, but then i couldn’t figure out how to make things work for the 3rd index.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 20, 2020, 5:55am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/5 "2020-07-20T05:55:49Z")

</div>

> [@purplishrock](#):
>
> yes it does! How ?!

“Julia provides broadcast, which expands singleton dimensions in array arguments to match the corresponding dimension in the other array without using extra memory, and applies the given function elementwise”  
[https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting-1](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting-1)

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [July 20, 2020, 6:25am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/6 "2020-07-20T06:25:16Z")

</div>

> [@purplishrock](#):
>
> but then i couldn’t figure out how to make things work for the 3rd index.

You can do this, just putting the vector along the 3rd axis:

```julia
zsteps = reshape(range(0, nz-1, step=1), 1,1,:) 
a1[:,:,:,3] .= zsteps .* dz

```

---

<div class="post-metadata">

**Author:** ![purplishrock](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/purplishrock/32/13451_2.png) [@purplishrock](https://discourse.julialang.org/u/purplishrock)\
**Post date:** [July 21, 2020, 3:16am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/7 "2020-07-21T03:16:24Z")

</div>

Now I get it! LOL. Thank you both for your very helpful answers!

---

<div class="post-metadata">

**Author:** ![purplishrock](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/purplishrock/32/13451_2.png) [@purplishrock](https://discourse.julialang.org/u/purplishrock)\
**Post date:** [July 21, 2020, 4:15am UTC](https://discourse.julialang.org/t/looking-for-elegant-cartesian-grid-array-assignment/43360/8 "2020-07-21T04:15:40Z")

</div>

Doesn’t everyone like benchmarks ? 😉

Remarkably, the element-by-element assign and the precomputed vector assigns seem to be about the same speed. I’ll just include one result here, but i ran several different dimensions, as well as reversing the call order, and they track each other very closely.

```julia
  14.889 ms (3 allocations: 42.06 MiB)
[0.0005, 0.0005, 0.0005]
[0.2995, 0.17450000000000002, 0.0345]
(300, 175, 35, 3)
1837500
  13.552 ms (5 allocations: 42.06 MiB)
[0.0005, 0.0005, 0.0005]
[0.2995, 0.17450000000000002, 0.0345]
(300, 175, 35, 3)
1837500

```

```julia
using BenchmarkTools

function f0(nx, ny, nz, dx, dy, dz)
    a1 = Array{Float64,4}(undef, nx, ny, nz, 3)
    
    for k=1:nz
        for j=1:ny
            for i=1:nx
                a1[i,j,k,1] = (i-0.5)*dx
                a1[i,j,k,2] = (j-0.5)*dy
                a1[i,j,k,3] = (k-0.5)*dz
            end
        end
    end

    a1
end

                
function f1(nx, ny, nz, dx, dy, dz)
    
    a1 = Array{Float64,4}(undef, nx, ny, nz, 3)

    # println(a)

    # sequences of values that need only be calculated once
    
    xs = dx * (0.5 .+ range(0, nx-1, step=1))
    ys = reshape(dy * (0.5 .+ range(0, ny-1, step=1)), 1, :, 1)
    zs = reshape(dz * (0.5 .+ range(0, nz-1, step=1)), 1, 1, :)

    # assign x, y, z values in sequence
    
    a1[:,:,:,1] .= xs
    a1[:,:,:,2] .= ys
    a1[:,:,:,3] .= zs

    a1
end

function main()
    # electric fields are established at the node i,j,k exactly
    # the magnetic fields are offset 1/2 step

    nx = 300
    ny = 175
    nz = 35

    dx = 1E-3
    dy = 1E-3
    dz = 1E-3

    a = @btime f0($nx, $ny, $nz, $dx, $dy, $dz)
    println(a[1,1,1,:])
    println(a[end,end,end,:])
    println(size(a))
    println(nx*ny*nz)

    a = @btime f1($nx, $ny, $nz, $dx, $dy, $dz)
    println(a[1,1,1,:])
    println(a[end,end,end,:])
    println(size(a))
    println(nx*ny*nz)
end

main()

```