Хорошо, я думаю, что перепробовал все, что мог, и, к сожалению, я пришел к выводу, что не так уж много возможностей для оптимизаций, если только вы не готовы пойти на распараллеливание. Давайте посмотрим, почему, давайте посмотрим, что вы можете и не можете сделать.
Оптимизация компилятора
В настоящее время компиляторы чрезвычайно хороши в оптимизации кода, гораздо больше, чем люди. Опора на оптимизацию, выполняемую компиляторами, также дает дополнительное преимущество, заключающееся в том, что они не ухудшают читабельность вашего исходного кода. Что бы вы ни делали (при оптимизации по скорости), всегда пробуйте это с любой разумной комбинацией флагов компилятора. Вы даже можете попробовать несколько компиляторов. Лично я использовал только gfortran (входит в GCC) (операционная система - 64-битная Windows), который, как я полагаю, имеет эффективные и правильные методы оптимизации.
-O2 почти всегда резко улучшают скорость, но даже -O3 - безопасная ставка (среди прочего, она включает в себя восхитительную раскрутку петли). Для этой проблемы я также попробовал -ffast-math и -fexpensive-optimizations, они не имели какого-либо измеримого эффекта, но -march-corei7 (специфичная для архитектуры настройка процессора, специфичная для Core i7) имела, поэтому я сделал измерения с -O3 -march-corei7
Так как быстро это на самом деле?
Я написал следующий код для проверки вашего решения и скомпилировал его с -O3 -march-corei7. Обычно он работал менее 0,78-0,82 секунды.
program benchmark
implicit none
real :: start, finish
integer :: I, J, K
real :: SIDEBACK, SIDEOWN, SIDEFRONT
integer, parameter :: NX = 600
integer, parameter :: NY = 600
integer, parameter :: NZ = 600
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: STEN
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: RES
call random_number(STEN)
call cpu_time(start)
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
SIDEBACK = STEN(I-1,J-1,K-1) + STEN(I-1,J,K-1) + STEN(I-1,J+1,K-1) + &
STEN(I ,J-1,K-1) + STEN(I ,J,K-1) + STEN(I ,J+1,K-1) + &
STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1)
SIDEOWN = STEN(I-1,J-1,K) + STEN(I-1,J,K) + STEN(I-1,J+1,K) + &
STEN(I ,J-1,K) + STEN(I ,J,K) + STEN(I ,J+1,K) + &
STEN(I+1,J-1,K) + STEN(I+1,J,K) + STEN(I+1,J+1,K)
SIDEFRONT = STEN(I-1,J-1,K+1) + STEN(I-1,J,K+1) + STEN(I-1,J+1,K+1) + &
STEN(I ,J-1,K+1) + STEN(I ,J,K+1) + STEN(I ,J+1,K+1) + &
STEN(I+1,J-1,K+1) + STEN(I+1,J,K+1) + STEN(I+1,J+1,K+1)
RES(I,J,K) = ( SIDEBACK + SIDEOWN + SIDEFRONT ) / 27.0
END DO
END DO
END DO
call cpu_time(finish)
!Use the calculated value, so the compiler doesn't optimize away everything.
!Print the original value as well, because one can never be too paranoid.
print *, STEN(1,1,1), RES(1,1,1)
print '(f6.3," seconds.")',finish-start
end program
Хорошо, это то, что компилятор может взять нас. Что дальше?
Сохранить промежуточные результаты?
Как вы можете догадаться из вопросительного знака, этот на самом деле не работал. Сожалею. Но давайте не будем торопиться.
Как упоминалось в комментариях, ваш текущий код вычисляет каждую частичную сумму несколько раз, то есть STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1) одной итерации будет STEN(I,J-1,K-1) + STEN(I,J,K-1) + STEN(I,J+1,K-1) следующей итерации, поэтому нет необходимости извлекать и вычислять снова, вы можете сохранить эти частичные результаты.
Проблема в том, что мы не можем хранить слишком много частичных результатов. Как вы сказали, ваш код уже достаточно дружествен к кешу, каждая сохраненная вами частичная сумма означает на один элемент массива меньше, который вы можете сохранить в кеше L1. Мы можем хранить несколько значений из последних нескольких итераций I (значения для индекса I-2, I-3 и т. Д.), Но компилятор почти наверняка уже это сделает. У меня есть 2 доказательства этого подозрения. Во-первых, мое ручное развертывание цикла замедлило программу примерно на 5%
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX, 8
SIDEBACK(0) = STEN(I-1,J-1,K-1) + STEN(I-1,J,K-1) + STEN(I-1,J+1,K-1)
SIDEBACK(1) = STEN(I ,J-1,K-1) + STEN(I ,J,K-1) + STEN(I ,J+1,K-1)
SIDEBACK(2) = STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1)
SIDEBACK(3) = STEN(I+2,J-1,K-1) + STEN(I+2,J,K-1) + STEN(I+2,J+1,K-1)
SIDEBACK(4) = STEN(I+3,J-1,K-1) + STEN(I+3,J,K-1) + STEN(I+3,J+1,K-1)
SIDEBACK(5) = STEN(I+4,J-1,K-1) + STEN(I+4,J,K-1) + STEN(I+4,J+1,K-1)
SIDEBACK(6) = STEN(I+5,J-1,K-1) + STEN(I+5,J,K-1) + STEN(I+5,J+1,K-1)
SIDEBACK(7) = STEN(I+6,J-1,K-1) + STEN(I+6,J,K-1) + STEN(I+6,J+1,K-1)
SIDEBACK(8) = STEN(I+7,J-1,K-1) + STEN(I+7,J,K-1) + STEN(I+7,J+1,K-1)
SIDEBACK(9) = STEN(I+8,J-1,K-1) + STEN(I+8,J,K-1) + STEN(I+8,J+1,K-1)
SIDEOWN(0) = STEN(I-1,J-1,K) + STEN(I-1,J,K) + STEN(I-1,J+1,K)
SIDEOWN(1) = STEN(I ,J-1,K) + STEN(I ,J,K) + STEN(I ,J+1,K)
SIDEOWN(2) = STEN(I+1,J-1,K) + STEN(I+1,J,K) + STEN(I+1,J+1,K)
SIDEOWN(3) = STEN(I+2,J-1,K) + STEN(I+2,J,K) + STEN(I+2,J+1,K)
SIDEOWN(4) = STEN(I+3,J-1,K) + STEN(I+3,J,K) + STEN(I+3,J+1,K)
SIDEOWN(5) = STEN(I+4,J-1,K) + STEN(I+4,J,K) + STEN(I+4,J+1,K)
SIDEOWN(6) = STEN(I+5,J-1,K) + STEN(I+5,J,K) + STEN(I+5,J+1,K)
SIDEOWN(7) = STEN(I+6,J-1,K) + STEN(I+6,J,K) + STEN(I+6,J+1,K)
SIDEOWN(8) = STEN(I+7,J-1,K) + STEN(I+7,J,K) + STEN(I+7,J+1,K)
SIDEOWN(9) = STEN(I+8,J-1,K) + STEN(I+8,J,K) + STEN(I+8,J+1,K)
SIDEFRONT(0) = STEN(I-1,J-1,K+1) + STEN(I-1,J,K+1) + STEN(I-1,J+1,K+1)
SIDEFRONT(1) = STEN(I ,J-1,K+1) + STEN(I ,J,K+1) + STEN(I ,J+1,K+1)
SIDEFRONT(2) = STEN(I+1,J-1,K+1) + STEN(I+1,J,K+1) + STEN(I+1,J+1,K+1)
SIDEFRONT(3) = STEN(I+2,J-1,K+1) + STEN(I+2,J,K+1) + STEN(I+2,J+1,K+1)
SIDEFRONT(4) = STEN(I+3,J-1,K+1) + STEN(I+3,J,K+1) + STEN(I+3,J+1,K+1)
SIDEFRONT(5) = STEN(I+4,J-1,K+1) + STEN(I+4,J,K+1) + STEN(I+4,J+1,K+1)
SIDEFRONT(6) = STEN(I+5,J-1,K+1) + STEN(I+5,J,K+1) + STEN(I+5,J+1,K+1)
SIDEFRONT(7) = STEN(I+6,J-1,K+1) + STEN(I+6,J,K+1) + STEN(I+6,J+1,K+1)
SIDEFRONT(8) = STEN(I+7,J-1,K+1) + STEN(I+7,J,K+1) + STEN(I+7,J+1,K+1)
SIDEFRONT(9) = STEN(I+8,J-1,K+1) + STEN(I+8,J,K+1) + STEN(I+8,J+1,K+1)
RES(I ,J,K) = ( SIDEBACK(0) + SIDEOWN(0) + SIDEFRONT(0) + &
SIDEBACK(1) + SIDEOWN(1) + SIDEFRONT(1) + &
SIDEBACK(2) + SIDEOWN(2) + SIDEFRONT(2) ) / 27.0
RES(I + 1,J,K) = ( SIDEBACK(1) + SIDEOWN(1) + SIDEFRONT(1) + &
SIDEBACK(2) + SIDEOWN(2) + SIDEFRONT(2) + &
SIDEBACK(3) + SIDEOWN(3) + SIDEFRONT(3) ) / 27.0
RES(I + 2,J,K) = ( SIDEBACK(2) + SIDEOWN(2) + SIDEFRONT(2) + &
SIDEBACK(3) + SIDEOWN(3) + SIDEFRONT(3) + &
SIDEBACK(4) + SIDEOWN(4) + SIDEFRONT(4) ) / 27.0
RES(I + 3,J,K) = ( SIDEBACK(3) + SIDEOWN(3) + SIDEFRONT(3) + &
SIDEBACK(4) + SIDEOWN(4) + SIDEFRONT(4) + &
SIDEBACK(5) + SIDEOWN(5) + SIDEFRONT(5) ) / 27.0
RES(I + 4,J,K) = ( SIDEBACK(4) + SIDEOWN(4) + SIDEFRONT(4) + &
SIDEBACK(5) + SIDEOWN(5) + SIDEFRONT(5) + &
SIDEBACK(6) + SIDEOWN(6) + SIDEFRONT(6) ) / 27.0
RES(I + 5,J,K) = ( SIDEBACK(5) + SIDEOWN(5) + SIDEFRONT(5) + &
SIDEBACK(6) + SIDEOWN(6) + SIDEFRONT(6) + &
SIDEBACK(7) + SIDEOWN(7) + SIDEFRONT(7) ) / 27.0
RES(I + 6,J,K) = ( SIDEBACK(6) + SIDEOWN(6) + SIDEFRONT(6) + &
SIDEBACK(7) + SIDEOWN(7) + SIDEFRONT(7) + &
SIDEBACK(8) + SIDEOWN(8) + SIDEFRONT(8) ) / 27.0
RES(I + 7,J,K) = ( SIDEBACK(7) + SIDEOWN(7) + SIDEFRONT(7) + &
SIDEBACK(8) + SIDEOWN(8) + SIDEFRONT(8) + &
SIDEBACK(9) + SIDEOWN(9) + SIDEFRONT(9) ) / 27.0
END DO
END DO
END DO
И что еще хуже, легко показать, что мы уже достаточно близки к теоретически минимально возможному времени выполнения. Чтобы вычислить все эти средние значения, абсолютный минимум, который нам нужно сделать, это получить доступ к каждому элементу хотя бы один раз и разделить их на 27,0. Таким образом, вы никогда не сможете получить быстрее, чем следующий код, который выполняется на моей машине менее 0,48-0,5 секунд.
program benchmark
implicit none
real :: start, finish
integer :: I, J, K
integer, parameter :: NX = 600
integer, parameter :: NY = 600
integer, parameter :: NZ = 600
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: STEN
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: RES
call random_number(STEN)
call cpu_time(start)
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
!This of course does not do what you want to do,
!this is just an example of a speed limit we can never surpass.
RES(I, J, K) = STEN(I, J, K) / 27.0
END DO
END DO
END DO
call cpu_time(finish)
!Use the calculated value, so the compiler doesn't optimize away everything.
print *, STEN(1,1,1), RES(1,1,1)
print '(f6.3," seconds.")',finish-start
end program
Но даже отрицательный результат есть результат. Если простой доступ к каждому элементу один раз (и деление на 27,0) занимает более половины времени выполнения, это просто означает, что доступ к памяти является узким местом. Тогда, может быть, вы можете оптимизировать , что .
Меньше данных
Если вам не нужна полная точность 64-битных двойных чисел, вы можете объявить свой массив с типом real(kind=4). Но, возможно, ваши реалы уже 4 байта. В этом случае я считаю, что некоторые реализации на Фортране поддерживают нестандартные 16-битные двойные числа, или, в зависимости от ваших данных, вы можете просто использовать целые числа (возможно, числа с плавающей точкой, умноженные на число, затем округленные до целых). Чем меньше ваш базовый тип, тем больше элементов вы можете поместить в кеш. Самым идеальным было бы integer(kind=1), конечно, это вызвало увеличение скорости моей машины более чем в 2 раза по сравнению с real(kind=4). Но это зависит от точности, которая вам нужна.
Лучшая местность
Основные массивы столбцов работают медленно, когда вам нужны данные из соседнего столбца, а основные массивы столбцов работают медленно для соседних строк.
К счастью, существует причудливый способ хранения данных, называемый Кривая Z-порядка , который имеет приложений, аналогичных вашему сценарию использования в компьютерной графике.Я не могу обещать, что это поможет, может быть, это будет ужасно контрпродуктивно, но, возможно, нет. Извините, мне самому не хотелось реализовывать это, если честно.
Распараллеливание
Говоря о компьютерной графике, эта проблема тривиально и очень хорошо распараллеливается, может быть, даже на GPU, но если вы не хотите заходить так далеко, вы можете просто использовать обычный многоядерный процессор. Fortran Wiki кажется хорошим местом для поиска библиотек распараллеливания Fortran.