accessing boundaries of array without duplicating lots of code

Question

When I try to compile my code using -fcheck=all I get a runtime error since it seems I step out of bounds of my array dimension size. It comes from the part of my code shown below. I think it is because my loops over i,j only run from -ny to ny, -nx to nx but I try to use points at i+1,j+1,i-1,j-1 which takes me out of bounds in my arrays. When the loop over j starts at -ny, it needs j-1, so it immediately takes me out of bounds since I'm trying to access -ny-1. Similarly when j=ny, i=-nx,nx.

My question is, how can I fix this problem efficiently using minimal code?

I need the array grad(1,i,j) correctly defined on the boundary, and it needs to be defined exactly as on the right hand side of the equality below, I just don't know an efficient way of doing this. I can explicitly define grad(1,nx,j), grad(1,-nx,j), etc, separately and only loop over i=-nx+1,nx-1,j=-ny+1,ny-1 but this causes lots of duplicated code and I have many of these arrays so I don't think this is the logical/efficient approach. If I do this, I just end up with hundreds of lines of duplicated code that makes it very hard to debug. Thanks.

integer :: i,j
integer, parameter :: nx = 50, ny = 50
complex, dimension (3,-nx:nx,-ny:ny) :: grad,psi
real, parameter :: h = 0.1

do j = -ny,ny
do i = -nx,nx

    psi(1,i,j) = sin(i*h)+sin(j*h)
    psi(2,i,j) = sin(i*h)+sin(j*h)    
    psi(3,i,j) = sin(i*h)+sin(j*h)

end do
end do

do j = -ny,ny
do i = -nx,nx

    grad(1,i,j) = (psi(1,i+1,j)+psi(1,i-1,j)+psi(1,i,j+1)+psi(1,i,j-1)-4*psi(1,i,j))/h**2 & 

                - (psi(2,i+1,j)-psi(2,i,j))*psi(1,i,j)/h & 

                - (psi(3,i,j+1)-psi(3,i,j))*psi(1,i,j)/h &

                - psi(2,i,j)*(psi(1,i+1,j)-psi(1,i,j))/h &

                - psi(3,i,j)*(psi(1,i,j+1)-psi(1,i,j))/h

end do
end do

If I was to do this directly for grad(1,nx,j), grad(1,-nx,j), it would be given by

   do j = -ny+1,ny-1

       grad(1,nx,j) = (psi(1,nx,j)+psi(1,nx-2,j)+psi(1,nx,j+1)+psi(1,nx,j-1)-2*psi(1,nx-1,j)-2*psi(1,nx,j))/h**2 & 

                - (psi(2,nx,j)-psi(2,nx-1,j))*psi(1,nx,j)/h & 

                - (psi(3,nx,j+1)-psi(3,nx,j))*psi(1,nx,j)/h &

                - psi(2,nx,j)*(psi(1,nx,j)-psi(1,nx-1,j))/h &

                - psi(3,nx,j)*(psi(1,nx,j+1)-psi(1,nx,j))/h

       grad(1,-nx,j) = (psi(1,-nx+2,j)+psi(1,-nx,j)+psi(1,-nx,j+1)+psi(1,-nx,j-1)-2*psi(1,-nx+1,j)-2*psi(1,-nx,j))/h**2 & 

                - (psi(2,-nx+1,j)-psi(2,-nx,j))*psi(1,-nx,j)/h & 

                - (psi(3,-nx,j+1)-psi(3,-nx,j))*psi(1,-nx,j)/h &

                - psi(2,-nx,j)*(psi(1,-nx+1,j)-psi(1,-nx,j))/h &

                - psi(3,-nx,j)*(psi(1,-nx,j+1)-psi(1,-nx,j))/h

   end do

Answer 1

One possible way for you could be using an additional index variable for the boundaries, modified from the original index to avoid getting out-of-bounds. I mean something like this:

do j = -ny,ny
  jj = max(min(j, ny-1), -ny+1)
  do i = -nx,nx
    ii = max(min(i, nx-1), -nx+1)
    grad(1,i,j) = (psi(1,ii+1,j)+psi(1,ii-1,j)+psi(1,i,jj+1)+psi(1,i,jj-1)-4*psi(1,i,j))/h**2 &
                - (psi(2,ii+1,j)-psi(2,ii,j))*psi(1,i,j)/h &
                - (psi(3,i,jj+1)-psi(3,i,jj))*psi(1,i,j)/h &
                - psi(2,i,j)*(psi(1,ii+1,j)-psi(1,ii,j))/h &
                - psi(3,i,j)*(psi(1,i,jj+1)-psi(1,i,jj))/h
  end do
end do

It's hard for me to write a proper code because it seems you trimmed part of the original expression in the code you presented in the question, but I hope you understand the idea and apply it correctly for your logic.

Opinions:

1. Even though this is what you are asking for (as far as I understand), I would not recommend doing this before profiling and checking if assigning the boundary conditions manually after a whole array operation wouldn't be more efficient, instead. Maybe those extra calculations on the indices on each iteration could impact on performance (arguably less than if conditionals or function calls). Using "ghost cells", as suggested by @evets, could be even more performant. You should profile and compare.
2. I'd recommend you declaring your arrays as dimension(-nx:nx,-ny:ny,3) instead. Fortran stores arrays in column-major order and, as you are accessing values on the neighborhood of the "x" and "y", they would be non-contiguous memory locations for a fixed "other" dimension is the leftest, and that could mean less cache-hits.

Answer 2

In somewhat pseudo-code, you can do

do j = -ny, ny

   if (j == -ny) then
      p1jm1 = XXXXX    ! Some boundary condition
   else
      p1jm1 = psi(1,i,j-1)
   end if

   if (j == ny) then
      p1jp1 = YYYYY   ! Some other boundary condition
   else
      p1jp1 = psi(1,i,j+1)
   end if

   do i = -nx, ny
      grad(1,i,j) = ... term involving p1jm1 ...  term involving p1jp1 ...
      ...
   end do
end do

The j-loop isn't bad in that you are adding 2*2*ny conditionals. The inner i-loop is adding 2*2*nx conditionals for each j iteration (or 2*2*ny * 2*2*nx conditional). Note, you need a temporary for each psi with the triplet indices are unique, ie., psi(1,i,j+1), psi(1,i,j-1), and psi(3,i,j+1).