futhark/Lobby

Hi, I'm new to GPGPU programming and Futhark. It's been interesting, but I've hit a roadblock trying to optimize an algorithm. The problem is somewhat involved, so apologies for the length of this post.

Seam carving is an algorithm for "content-aware image resizing". I'm trying to implement "non-cumulative seam carving", an adaptation of the algorithm for GPUs described in this paper. Here's the relevant description: (pp. 3-4)

To summarize, there are two steps:

1. Create the index map. For each pixel, the map's value is the index of the adjacent pixel in the next row with the least energy.
2. Following the paths created by the indexes, sum up the energy for each seam (which runs from the top to the bottom of the image).

Creating the index map is just a map over the energy array. But their algorithm for summing seam energy uses a tree reduction that modifies SumMap and OffsetMap in place. It doesn't seem like Futhark can do this, as reduce does not operate on unique arrays. So, my current implementation creates temporary arrays: (energy_sum and next_index)

-- Note that I use a 2D array while the paper uses a 1D array, so my
-- indexes are x-coordinates.
let create_index_map [h][w] (energy: [h][w]u32): [h][w]i32 =
  map (\y ->
    let row = energy[y + 1]
    let mn i1 i2 = unsafe if row[i1] < row[i2] then i1 else i2
    in map (\x ->
      let min = x
      let min = if x != 0     then mn min (x - 1) else min
      let min = if x != w - 1 then mn min (x + 1) else min
      in min
    ) (iota w)
  ) (iota (h - 1)) ++ [replicate w 0]

let sum_seam_energy [h][w] (energy: [h][w]u32) (index: [h][w]i32): [w]u32 =
  let op (e1, i1) (e2, i2) =
    let energy_sum = map2 (\e i -> e + unsafe e2[i]) e1 i1
    let next_index = map (\i -> unsafe i2[i]) i1
    in (energy_sum, next_index)
  let ne = ((replicate w 0), (iota w))
  let as = zip energy index
  in (reduce op ne as).1

entry main [h][w] (energy: [h][w]u32): [w]u32 =
  sum_seam_energy energy (create_index_map energy)

Is there some way to use reduce but modify the array in place? I've looked at SOACs like reduce_by_index and scatter, but neither does what I want. (They each operate on arrays of destination indexes, but the indexes here are source indexes.) I ask this because on my laptop (Intel i5-5350U @ 4x 2.9GHz, Intel HD Graphics 6000), benchmarking OpenCL with a [2000][2000]u32 energy array gives a throughput of 0.1GB/s--4x slower than the CPU. My only guesses for the bad performance are poor hardware or too much allocation. Though the former is true, it seems like the latter is a bigger culprit. Also, if it's helpful, running with -D shows a peak memory usage of 48MB.

Thank you!