# Add specific elements of a CUDA matrix

**URL:** <https://discourse.julialang.org/t/add-specific-elements-of-a-cuda-matrix/111909>\
**Category:** GPU\
**Tags:** question, indexing, cuda, arithmetic\
**Created:** [March 20, 2024, 11:13pm UTC](https://discourse.julialang.org/t/add-specific-elements-of-a-cuda-matrix/111909 "2024-03-20T23:13:03Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![NonDairyNeutrino](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nondairyneutrino/32/221496_2.png) [@NonDairyNeutrino](https://discourse.julialang.org/u/NonDairyNeutrino)\
**Post date:** [March 20, 2024, 11:13pm UTC](https://discourse.julialang.org/t/add-specific-elements-of-a-cuda-matrix/111909/1 "2024-03-20T23:13:03Z")

</div>

# Question

I have an N \times N matrix A on the GPU (e.g. `CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}` from `CUDA.jl`), and two N-element “index vectors” `iVector` and `jVector` also on the GPU. I want to add the elements of A located at indices \{(i\_1, j\_1), (i\_2, j\_2), \ldots, (i\_N, j\_N)\}. With everything on the CPU, this is of course trivial with something like

```Julia
sum(A[iVector[k], jVector[k]] for k in 1:N)

```

There’s a similar method with everything on the GPU with

```Julia
sum(@inbounds view(A, view(iVector, index), view(jVector, index)) for index in eachindex(iVector))

```

but my gut says this isn’t the best way to do it and is more just a workaround. **Is there a better or more correct way to do this? What about doing the summation in parallel on the GPU as well? Possibly writing a custom CUDA kernal?**

# Example

```Julia
using Random: randperm
using CUDA
using BenchmarkTools

const dimension = 10^4
const adjacencyMatrix = CUDA.rand(dimension, dimension)
const tour = convert.(Int32, randperm(dimension)) |> cu

function fitness(
    adjacencyMatrix :: Union{Matrix{T}, CuArray}, 
    position :: Union{Vector, CuArray}
) where T <: Real
    head = view(position, 1:length(position)-1)
    tail = view(position, 2:length(position))
    return sum(
        @inbounds view(
            adjacencyMatrix, 
            view(head, index), 
            view(tail, index)
         ) for index in eachindex(head)
    )
end

display(@btime fitness(adjacencyMatrix, tour);)

```

Output:

> 119.928 ms (240377 allocations: 8.36 MiB)

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [March 21, 2024, 2:48pm UTC](https://discourse.julialang.org/t/add-specific-elements-of-a-cuda-matrix/111909/2 "2024-03-21T14:48:48Z")

</div>

You can take the `sum` of a `view` with indices:

```julia
function fitness(adjacencyMatrix::Union{Matrix{T}, CuArray},
                 position::Union{Vector, CuArray}) where T <: Real
    head = view(position, 1:length(position)-1)
    tail = view(position, 2:length(position))
    idx = CartesianIndex.(head, tail)
    sum(view(adjacencyMatrix, idx))
end

```

That’s slightly faster 🙂 100us instead of 270ms on my system.