доступ к границам массива без дублирования большого количества кода

Question

Когда я пытаюсь скомпилировать мой код, используя -fcheck=all, я получаю ошибку во время выполнения, поскольку кажется, что я выхожу за пределы размера моего массива. Это происходит из части моего кода, показанной ниже. Я думаю, это потому, что мои циклы над i, j бегут только от -ny до ny, -nx до nx, но я стараюсь использовать точки на i+1,j+1,i-1,j-1, что выводит меня за пределы моих массивов. Когда цикл над j начинается с -ny, ему нужно j-1, поэтому он сразу выводит меня за пределы, поскольку я пытаюсь получить доступ к -ny-1. Точно так же, когда j=ny, i=-nx,nx.

У меня вопрос, как я могу эффективно решить эту проблему, используя минимальный код?

Мне нужно правильно определить array grad(1,i,j) на границе, и оно должно быть определено точно так же, как в правой части приведенного ниже равенства, я просто не знаю эффективного способа сделать это. Я могу явно определить grad(1,nx,j), grad(1,-nx,j), etc, отдельно и только зацикливаться на i=-nx+1,nx-1,j=-ny+1,ny-1, но это вызывает много дублированного кода, и у меня есть много таких массивов, поэтому я не думаю, что это логический / эффективный подход. Если я сделаю это, я просто получу сотни строк дублированного кода, что делает его очень сложным для отладки. Спасибо.

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

Если бы я должен был сделать это напрямую для grad(1,nx,j), grad(1,-nx,j), это было бы дано

   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

Одним из возможных способов для вас может быть использование дополнительной индексной переменной для границ, измененной по сравнению с исходным индексом, чтобы избежать выхода за пределы. Я имею в виду что-то вроде этого:

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

Мне трудно написать правильный код, потому что кажется, что вы обрезали часть исходного выражения в коде, который вы представили в вопросе, но я надеюсь, что вы понимаете идею и правильно применяете ее для своей логики.

Мнения:

1. Даже если это то, о чем вы просите (насколько я понимаю), я бы не рекомендовал делать это перед профилированием и проверкой, не будет ли более эффективным присваивать граничные условия вручную после операции с целым массивом. Возможно, эти дополнительные вычисления индексов на каждой итерации могут повлиять на производительность (возможно, меньше, чем if условных выражений или вызовов функций). Использование "клеток-призраков", как предлагает @evets, может быть еще более эффективным. Вы должны профиль и сравнить.
2. Я бы порекомендовал вам объявить ваши массивы как dimension(-nx:nx,-ny:ny,3). Fortran хранит массивы в главном порядке столбцов, и, поскольку вы обращаетесь к значениям в окрестности «x» и «y», они будут несмежными ячейками памяти, поскольку фиксированное «другое» измерение является самым левым, и это может значит меньше кеш-хитов.

Answer 2

В некотором псевдокоде вы можете выполнить:

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

J-цикл неплох, поскольку вы добавляете 2 * 2 * ny условия.Внутренний i-цикл добавляет 2 * 2 * nx условных выражений для каждой j-й итерации (или 2 * 2 * ny * 2 * 2 * nx условных).Обратите внимание, вам нужен временный для каждого пси с триплетными индексами являются уникальными, т. Е. psi(1,i,j+1), psi(1,i,j-1) и psi(3,i,j+1).