Как оптимизировать эту подпрограмму Фортрана с большим количеством циклов? - PullRequest
Новая
       89

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

1 голос
/ 12 апреля 2020

Я имею дело с подпрограммой, которая очень неэффективна, когда размер массива становится большим, например, NN = 1000, KK = 200, MM = 200. Но я не могу придумать идеи для его оптимизации. 

program main

  implicit none

  integer :: NN, KK, MM
  integer, allocatable, dimension(:,:) :: id
  complex*16, allocatable, dimension(:) :: phase
  complex*16 :: phase_base(3)
  real*8, allocatable, dimension(:,:) :: wave_base

  complex*16, allocatable, dimension(:,:) :: wave
  integer :: i, j, k, n

  NN = 1000
  KK = 200
  MM = 200

  allocate(id(MM,3))
  allocate(phase(KK))
  allocate(wave_base(KK, NN*(NN+1)/2 ))
  allocate(wave(NN, NN))

  id(:,:) = 2

  phase_base(:) = (1.0d0,1.0d0)

  wave_base(:,:) = 1.0d0

  phase(:) = (1.0d0,1.0d0)

  call  noise_wave(NN, KK, MM, id, phase, phase_base, wave_base, wave)

  deallocate(id)
  deallocate(phase)
  deallocate(wave_base)
  deallocate(wave)

end program main

subroutine noise_wave(NN, KK, MM, id, phase_1, phase_base, wave_base, wave)
  implicit none

  integer, intent(in) :: NN, KK, MM
  integer, intent(in), dimension(MM, 3) :: id
  complex*16, intent(in) :: phase_1(KK)
  complex*16, intent(in) :: phase_base(3)
  real*8,  intent(in) :: wave_base(KK, NN*(NN+1)/2 )

  complex*16, intent(out) :: wave(NN, NN)

  integer :: i, j, k, p, n
  integer :: x, y, z
  real :: start, finish
  complex*16 :: phase_2, phase_2_conjg

  do p = 1, MM

    x = id(p, 1)
    y = id(p, 2)
    z = id(p, 3)

    phase_2 = (phase_base(1) ** x) * (phase_base(2) ** y) * (phase_base(3) ** z)

    phase_2_conjg = conjg(phase_2)

    n = 0
    do j = 1, NN
      do i = 1, j   ! upper triangle

        n = n + 1

        do k = 1, KK

          wave(i,j) = wave(i,j) + wave_base(k,n) * phase_1(k) * phase_2_conjg

        enddo

        wave(j,i) = conjg(wave(i,j) )

      enddo
    enddo
  enddo

end subroutin

Может ли кто-нибудь дать мне подсказку? (Я выполнил предложенную оптимизацию. Кроме того, по предложению Яна я добавил небольшой тест. Таким образом, вы можете проверить его напрямую.)

Ответы [ 2 ]

3 голосов
/ 12 апреля 2020

Вы можете получить измеримое ускорение, если измените свое гнездо l oop на 

  do p = 1, MM

     x = id(p, 1)
     y = id(p, 2)
     z = id(p, 3)
     phase = (phase_base(1) ** x) * (phase_base(2) ** y) * (phase_base(3) ** z)
     conjg_phase = conjg(phase)  ! new variable, calculate here, use below

     n = 0
     do j = 1, NN
        do i = 1, j   
           n = n + 1
           do k = 1, KK
              wave(i,j) = wave(i,j) + wave_base(k,n) * conjg_phase
           enddo
        enddo
        wave(j,i) = conjg(wave(i,j) )
     enddo
  enddo

(и оно все равно может быть правильным, если я понял код!). Даже небольшие вычисления, подобные тем, которые я вытащил из нижней части гнезда l oop, - это перетаскивание, если они повторяются достаточно часто. И скорость выполнения может также выиграть от перемещения этих значений в кеш и из кеша реже.

Возможно, стоит (немного) поменять местами измерения id, тогда чтение id(1:3,p), вероятно, будет более дружественным к кэшу, чем текущая версия.

И если выполнение скорость по-прежнему не по вкусу, время изучать OpenMP (если вы его еще не знаете).

1 голос
/ 12 апреля 2020

Вот мое решение, следуя приведенным выше хорошим идеям. До OpenMP еще есть место для повышения эффективности. Например, первый kl oop в подпрограмме можно исключить с помощью функции sum. 

program main

  implicit none

  integer :: NN, KK, MM
  integer, allocatable, dimension(:,:) :: id
  complex*16, allocatable, dimension(:) :: phase
  complex*16 :: phase_base(3)
  real*8, allocatable, dimension(:,:) :: wave_base

  complex*16, allocatable, dimension(:,:) :: wave
  integer :: i, j, k, n

  NN = 1000
  KK = 200
  MM = 200

  allocate(id(MM,3))
  allocate(phase(KK))
  allocate(wave_base(KK, NN*(NN+1)/2 ))
  allocate(wave(NN, NN))

  id(:,:) = 2

  phase_base(:) = (1.0d0,1.0d0)

  wave_base(:,:) = 1.0d0

  phase(:) = (1.0d0,1.0d0)

  call  noise_wave(NN, KK, MM, id, phase, phase_base, wave_base, wave)

  deallocate(id)
  deallocate(phase)
  deallocate(wave_base)
  deallocate(wave)

end program main

subroutine noise_wave(NN, KK, MM, id, phase_1, phase_base, wave_base, wave)
  implicit none

  integer, intent(in) :: NN, KK, MM
  integer, intent(in), dimension(MM, 3) :: id
  complex*16, intent(in) :: phase_1(KK)
  complex*16, intent(in) :: phase_base(3)
  real*8,     intent(in) :: wave_base(KK, NN*(NN+1)/2 )
  complex*16, intent(out):: wave(NN, NN)

  integer :: i, j, k, p, n
  integer :: x, y, z
  real :: start, finish
  complex*16 :: phase_2, phase_2_conjg
  complex*16 :: wave_tmp(NN*(NN+1)/2)
  complex*16 :: wave_tmp_2(NN*(NN+1)/2)

  do k = 1, KK

    wave_tmp(:) = wave_tmp(:) + wave_base(k,:) * phase_1(k)

  enddo

  do p = 1, MM
    phase_2 = product(phase_base(:)**id(p,:) )
    phase_2_conjg = conjg(phase_2)

    wave_tmp2(:) = wave_tmp2(:) + wave_tmp(n) * phase_2_conjg
  enddo

  n = 0
  do j = 1, NN
    do i = 1, j
        n = n + 1
        wave(i,j) = wave_tmp2(n)
        wave(j,i) = conjg(wave_tmp2(n) )
    enddo
  enddo

end subroutine
...