Unsure why a function isn't assigning a new value to a view

General Usage
1

I am writing code that goes through an array and does a calculation on each element, with the results previous calculations being used for future ones.

I have a function that did this one index at a time, step!, and it works fine, but it’s slow.

@views function step!(arr, U, M, dx, dt)
  (Nx, Ny) = size(arr)
  I1 = CartesianIndex(1, 1)
  for i in 2:(Nx-1)
    for j in 2:(Ny-1)
      loc = CartesianIndex(j, i)
      @inbounds calc!(arr[loc-I1:loc+I1], U, M, dx, dt)
    end
  end
end

I then wrote a function to calculate the indices that could be done all at the same time using SIMD, but were dependent on a previous calculation. However, this function doesn’t allocate to arr and I am not sure why.

@views function step_loc!(arr, inds, U, M, dx, dt)
  I1 = CartesianIndex(1, 1)
  for it in inds
    @simd for loc in it
      @inbounds calc!(arr[loc-I1:loc+I1], U, M, dx, dt)
    end
  end
end

I’m wondering what the difference between these two functions is. calc! hasn’t changed, but it assigns as expected in step!,. but not in step_loc!.

calc! has the signature

@views function calc(arr, U, M, dx, dt)

and assigns specifically to arr[5] by

arr[5] = result

Thank you in advance for the help!

2

Are you sure @views encapsulates the arr[loc-I1:loc+I1] in your inner loop? What happens if you move the @views to that indexing location?

From what you’ve given, it should work. Perhaps there’s a hidden copy in calc! somewhere? Hard to tell without a complete MWE.

3

I tried this

function step_loc!(arr, inds, U, M, dx, dt)
  I1 = CartesianIndex(1, 1)

  for it in inds
    @simd for loc in it
      @inbounds @views calc(arr[loc-I1:loc+I1], U, M, dx, dt)
    end
  end
end

Did you mean something different?

and, the calc code is the following:

@views function calc!(arr, U, M, dx, dt)
  corner_sum = arr[1] + arr[3] + arr[7] + arr[9]
  edge_sum = arr[2] + arr[4] + arr[6] + arr[8]

  ij = arr[5]

  grad = (2.0 / 3.0) * (edge_sum + 0.25 * corner_sum - 5.0 * ij) / dx / dx

  t1 = ij - ij * ij * ij
  t2 = -U * (1 - ij * ij) * (1 - ij * ij)

  dFdt = grad + t1 + t2

  arr[5] = ij + M * dFdt * dt
end
4

Hm, that looks fine - what indices are you calculating for inds? What are the exact differences you observe between the two functions - does step_loc! not assign at all? Do you perhaps have a complete example we could run on our end?

5 
function get_diag(idx, arr)
  Ny, _ = size(arr)
  I2Right = CartesianIndex(0, 2)
  I1Up = CartesianIndex(1, 0)
  result = CartesianIndex[]
  index = idx

  while true
    push!(result, index)
    index += I2Right
    index -= I1Up
    all(i -> i > 1 && i < Ny, Tuple(index)) || break
  end

  return result

end

@views function indexes!(x, out)
  Ny, Nx = size(x)
  inds = CartesianIndices(x)
  lin_inds = LinearIndices(x)

  r1 = lin_inds[inds[2, begin+1:end-1]]
  r2 = lin_inds[inds[3, begin+1:end-1]]
  interleaved = inds[[r1 r2]'[:]]

  for i in union(interleaved, 2Ny-1:Ny:(Ny*Nx)-Ny)
    push!(out, get_diag(inds[i], x))
  end

end

@views function process!(arr, inds, U, M, dx, dt)
  I = CartesianIndex(1, 1)
  for locs in inds
    @simd for loc in locs
      @inbounds calc(arr[loc-I:loc+I], U, M, dx, dt)
    end
  end
end

@views function calc(arr, U, M, dx, dt)
  corner_sum = arr[1] + arr[3] + arr[7] + arr[9]
  edge_sum = arr[2] + arr[4] + arr[6] + arr[8]

  ij = arr[5]

  grad = (2.0 / 3.0) * (edge_sum + 0.25 * corner_sum - 5.0 * ij) / dx / dx

  t1 = ij - ij * ij * ij
  t2 = -U * (1 - ij * ij) * (1 - ij * ij)

  dFdt = grad + t1 + t2
  arr[5] = ij + M * dFdt * dt
end

U = -0.5
M = 1.0
dx = 0.5
dt = 0.03
Ny, Nx = 514, 514

result = Vector{CartesianIndex}[]

test_arr = Matrix{Float64}(undef, Ny, Nx)
indexes!(test_arr, result)


  process!(test_arr, result, U, M, dx, dt)

test_arr

I made a minimal example, but this works? test_arr gets assigned in this one.

6

Looks like the issue was in my get_diag! function! Thanks for all the help, even though it was just my silly mistake.

2 Likes