# Optimization of nonlocal stencil-like operator on subset of regular grid

I am trying to optimize the execution time for this particular piece of fortran code.

Details: i_gc is a (ngpts, 3) array of containing (i,j,k) indices for each grid point. This is a subset of the full NX x NY x NZ grid. The subset is connected in the sense that density of sampling varies, but the points aren't scattered randomly throughout the domain.

Essentially I have a non-local operator on the grid which is a tensor product of 3 one dimensional operators. If I had the data on the full 3d grid it would be pretty easy to write this code without all the if statements in the inner loop. Also, I would like to avoid storing a set {i},{j},{k} indices for each grid point, but maybe this is the way to go?

I am wondering if there is a way to do much better than this? I do have some freedom to reorder the points if needed. Maybe I can exploit that somehow? I suppose I could try something like a Hilbert curve, but I'm not sure how well that could apply to my case where I have an oddly shaped subset of a regular grid.

edit: mathematically the operation is like this

$Y_{ijk} = \sum^{N^x}_{q=1} M^x_{qi} X_{qjk} + \sum^{N^y}_{q=1} M^y_{qj} X_{iqk} + \sum^{N^z}_{q=1} M^z_{qk} X_{ijq}$

where the sums are only over points in the subset (i_gc)

subroutine apply_lf4(ngpts, i_gc, mx,my,mz,x,y)
implicit none
integer, intent(in) :: ngpts
integer, intent(in), dimension(:,:) :: i_gc
real(8), intent(in), dimension(:,:) :: mx,my,mz
real(8), intent(in), dimension(ngpts) :: x
real(8), intent(out), dimension(ngpts) :: y

integer :: ig,jg,i1,i2,i3,j1,j2,j3

y=0.0d0

do ig=1, ngpts
i1=i_gc(ig,1); i2=i_gc(ig,2); i3=i_gc(ig,3)
inner: do jg=1, ngpts
j1=i_gc(jg,1)
if (j1==i1) then
j2=i_gc(jg,2)
if (j2==i2) then
j3 = i_gc(jg,3)
y(ig) = y(ig) + mz(j3,i3)*x(jg)
cycle inner
end if
j3=i_gc(jg,3)
if (j3==i3) then
y(ig) = y(ig) + my(j2,i2)*x(jg)
cycle inner
end if
else
j2 = i_gc(jg,2)
if (j2 == i2) then
j3 = i_gc(jg,3)
if (j3 == i3) then
y(ig) = y(ig) + mx(j1,i1)*x(jg)
cycle inner
end if
end if
end if
end do inner
end do

end subroutine apply_lf4

• How large is ngpts relative to NX, NY, NZ? – Jesse Chan May 6 '15 at 0:21
• What is this computing in mathematical terms? I'm a bit lost. – Bill Barth May 6 '15 at 0:22
• @JesseChan a typical NX=NY=NZ=400, Ngpts=10^6 – user1984528 May 6 '15 at 2:15
• @BillBarth I have updated the question to include some mathematical terms. Hopefully that is clear enough? – user1984528 May 6 '15 at 2:22