I am trying to populate a large, dense 3D matrix (3000 x 18 x 1000). I tried iterating over all positions of the 3D matrix (testing various row/column iteration orders) - see Attempt 1. I also tried array comprehensions on an array of matrices - see Attempt 2. The timing for both is about the same (see below). Does anyone have suggestions for how to further speed this up? The examples below have the first dimension at 50 instead of 3000, but I need to scale up to 3000 (see the `pointNumber`

dimension in the code).

**Attempt 1 - iterate over all positions of 3D array with nested for loop**

## Summary

```
using LinearAlgebra
using BenchmarkTools
function greensMatrix3D1(matchCoordinates,axisCoordinates,zArray)
greensMatrix = zeros(Float64,length(matchCoordinates),length(axisCoordinates),length(zArray))
for v = 1:length(matchCoordinates),q = 1:length(axisCoordinates),k=1:length(zArray)
rvec = matchCoordinates[v] - [axisCoordinates[q][1],axisCoordinates[q][2],zArray[k]]
rivec = rvec + [0, 2*axisCoordinates[q][2],0]
r = norm(rvec)
ri = norm(rivec)
greensMatrix[v,q,k] = (1/r - 1/ri)
end
return greensMatrix
end
pointNumber = 50 #target 3000
matchCoordinates = [[10,i*0.02,0] for i in 1:pointNumber]
axisCoordinates = [[i-1,10,0] for i in 1:18]
zArray = range(-100,100,length=1024)
@btime output1 = greensMatrix3D1(matchCoordinates,axisCoordinates,zArray)
```

**Attempt 2 - use comprehensions and an array of matrices**

## Summary

```
using LinearAlgebra
function oneGreens(xyz1,xyz2)
rvec = xyz1 - xyz2
rivec = rvec + [0, 2*xyz2[2],0]
r = norm(rvec)
ri = norm(rivec)
greensMatrix = (1/r - 1/ri)
return greensMatrix
end
function greensMatrix3D(matchCoordinates,axisCoordinates,zArray)
greensMatrixArray = [zeros(Float64,length(matchCoordinates),length(axisCoordinates)) for _ in 1:length(zArray)]
for i = 1:length(zArray)
axisTempCoordinates = axisCoordinates
setindex!.(axisTempCoordinates,zArray[i],3)
greensMatrixArray[i] = [oneGreens(j,k) for j in matchCoordinates, k in axisTempCoordinates]
end
return greensMatrixArray
end
pointNumber = 50 #target 3000
matchCoordinates = [[10.0,i*0.02,0.0] for i in 1:pointNumber]
axisCoordinates = [[i-1.0,10.0,0.0] for i in 1:18]
zArray = range(-100.0,100.0,length=1024)
@btime output2 = greensMatrix3D(matchCoordinates,axisCoordinates,zArray)
```

The timing results for the above are:

Attempt 1:>>> `1.826 s (8294402 allocations: 513.28 MiB)`

Attempt 2:>>> `1.673 s (7375873 allocations: 422.41 MiB)`