Замена циклов do-while в качестве запроса конвергенции в Fortran OpenMP - PullRequest
0 голосов
/ 24 апреля 2018

Я создал программу, которая использует дифференциальную эволюцию для оптимизации позиций атомов с точки зрения их парного потенциала, и теперь хочу распараллелить ее с OpenMP, для которого я совсем новичок. Дифференциальная эволюция использует общий цикл do-while, в котором запрос конвергенции используется в качестве условия выхода.

Это значит

  1. Я знаю, что не могу просто !$OMP PARALLEL DO цикл do-while
  2. Я не могу предсказать, в какой момент цикл прерывается
  3. следующие итерации также будут соответствовать условию. Следующее мой непараллельный код:

    implicit none
    integer, parameter :: dp = selected_real_kind(15, 307)
    real(kind = dp) :: sum_dif, sum, sum_old, final_toler, F, CR, d1, d2, d3, max_step, this_step, scale, randomtest
    integer :: pop, dim, arfa, beta, gama, delt, bi, jrand, kf, ki, kj, kg, dim_low, i, g, num_Ti, idum
    logical :: monitor_progress, history
    integer, dimension(0:) :: index
    real(kind=dp), dimension(0:,0:) :: depop, tempop 
    real(kind=dp), dimension(0:) :: best,temp_best, proj, d_fx, t_fx
    real(kind=dp) :: midpoint
    logical :: best_DE, use_maxstep
    integer*2 :: best_DE_num
    
    sum_dif = 0.0
    do while ((abs(sum_old - sum + sum_dif)) > final_toler) !UTILIZE DIFFERENTIAL EVOLUTION until convergence
        ! ↑ enclosing convergence loop
        sum_old = sum
        !**initialize DE-Parameters**
        if (best_DE) then
            F = 0.2_dp + 0.6_dp * ran2(idum)
            CR = 0.6_dp + 0.4_dp * ran2(idum)
        else
            F = 0.4_dp + 0.4_dp * ran2(idum)
            CR = 0.7_dp + 0.2_dp * ran2(idum)
        end if
        !******
        !** Create Mutant Population AND Perform Population Update**
        do i = 0, pop - 1
            do
                arfa = RNGgen(pop)
                if (arfa /= i) exit
            end do
            do
                beta = RNGgen(pop)
                if (beta /= arfa) then
                    if (beta /= i) exit
                end if
            end do
            do
                gama = RNGgen(pop)
                if (gama /= beta) then
                    if (gama /= arfa) then
                        if (gama /= i) exit
                    end if
                end if
            end do
            do
                delt = RNGgen(pop)
                if (delt /= gama) then
                    if (delt /= beta) then
                        if (delt /= arfa) then
                            if (delt /= i) exit !loops will be continued until arfa!=beta!=gama!=delt!=i
                        end if
                    end if
                end if
            end do
            jrand = RNGgen(dim + 1)/3 
            !**Create mutant population per atom not per dim
            kf = 0
    
            do while (kf < dim)
                randomtest = ran2(idum)
                if ((randomtest <= CR).or.(kf == jrand)) then 
                    !**Create hybrid child**
                d1 = F * (depop(arfa, kf) - depop(gama, kf) + best_DE_num * (depop(beta, kf) - depop(delt, kf)))
                d2 = F * (depop(arfa, kf + 1) - depop(gama, kf + 1) + best_DE_num * (depop(beta, kf + 1) - depop(delt, kf + 1)))
                d3 = F * (depop(arfa, kf + 2) - depop(gama, kf + 2) + best_DE_num * (depop(beta, kf + 2) - depop(delt, kf + 2)))
                    !******
                    if (use_maxstep) then
                        this_step = d1 * d1 + d2 * d2 + d3 * d3!norm^2
                        if (this_step > max_step) then
                            scale = sqrt(max_step/this_step)
                            d1 = d1 * scale
                            d2 = d2 * scale
                            d3 = d3 * scale
                        end if
                    end if !end if use_maxstep
                    tempop(i, kf) = best_DE_num * best(kf)+(1 - best_DE_num) * depop(beta, kf) + d1
                    tempop(i, kf + 1) = best_DE_num * best(kf + 1)+(1 - best_DE_num) * depop(beta, kf + 1) + d2
                    tempop(i, kf + 2) = best_DE_num * best(kf + 2)+(1 - best_DE_num) * depop(beta, kf + 2) + d3
                else
                    tempop(i, kf) = depop(i, kf)
                    tempop(i, kf + 1) = depop(i, kf + 1)
                    tempop(i, kf + 2) = depop(i, kf + 2)
                end if
                kf = kf + 3
            end do !end dim do loop
            !******
            call trans_to_cent(tempop(i,:), midpoint, midpoint, midpoint, 0, dim_low)
            !******
            !**Evaluate Objective Function for Mutant
            tempop(i, dim) = Analytic_U(num_Ti, tempop(i,:))
            t_fx(i) = tempop(i, dim) !store tempop fx values
            !******
        end do !end pop do loop
        do i = 0, pop - 1
            d_fx(i) = depop(i, dim) !store depop fx values
        end do
        !******
        !**SORTED MUTANT REPLACEMENT**
        index = Max_DelF(d_fx, t_fx)
        do kg = 0, pop - 1
            if (tempop(index(kg), dim) < depop(kg, dim)) then
                depop(kg,:) = tempop(index(kg),:)
            end if
            d_fx(kg) = depop(kg, dim) 
        end do
        !******
        !**Obtain total cost function values for tolerance comparison**
        sum = 0
        do ki = 0, pop - 1
            sum = sum + depop(ki, dim)
        end do
        sum = sum/pop !calculate average energy value            
        !******
        !**Obtain global best**
        do kj = 0, pop - 1
            if (best(dim) > depop(kj, dim)) then
                best = depop(kj,:)
                bi = kj
            end if !determine "best" vector
        end do
        !******
        if (monitor_progress) then
            print*, "Progress (min = 1.0): ", (abs(sum_old - sum + sum_dif))/final_toler
        end if
        g = g + 1 !increment iteration counter
    end do !end generation (while) loop
    

Здесь указанная петля сверху. Условие выхода срабатывает, когда разность энергий от одной итерации к следующей ниже определенного порога. Код включает в себя несколько других циклов do внутри этого цикла, но они должны быть распараллеливаемыми без особых проблем.

Теперь у меня вопрос: можно ли распараллелить вмещающий цикл, не отказываясь от большей части повышения производительности, или все еще будет повышение, если я попытаюсь исключить его из параллельной области?

Если бы я сделал это, я также мог бы исключить генерацию мутантных популяционных переменных arfa, beta, gama, delt, потому что их генерация должна была бы быть сделана с !$OMP CRITICAL со значительными накладными расходами, так как они случайным образом строятся с arfa!=beta!=gama!=delt!, верно? Я использую random_number встроенный со случайными семенами в моей функции RNGgen. Мой компилятор gfortran.

1 Ответ

0 голосов
/ 29 апреля 2018

Пример, который вы предоставили, был неполным: он не компилировался. Я решил построить небольшой (полная) программа, которая делает что-то похожее. Надеюсь, это поможет!

Программа запускает параллельный сеанс, в котором обнаруживаются новые группы населения. Всякий раз, когда найдено лучшее население, чем лучшее, лучшее население обновляется. Итерация останавливается, когда глобальные вычисления проводят слишком много итераций. между последовательными улучшениями.

В этой программе каждое следующее население построено полностью с нуля. В вашей программе более развитое поколение следующего населения. в пользу «лучшие» популяции по сравнению с «худшими».

В наивной параллеллизации каждый поток будет следовать своему поиску пространство, и он не будет «учиться» из того, что обнаружили другие темы. Для обмена поисковой информацией между потоками вам необходимо разработать метод, а затем запрограммируйте это. Казалось, что это выходит за рамки этого вопроса.

А вот и программа:

program optimize_population
use Population ! contains the target function
use omp_lib
implicit none
    integer, parameter :: really_long = 100000
    real(kind=dp) :: best, targetFun
    real(kind=dp) :: pop(2,npop), best_pop(2,npop)
    integer, allocatable :: iter(:)
    integer :: i, nt, it, ierr, last_improvement 

    call initTarget()  ! initialize the target function

    ! Allocate an iteration count for every thread
    nt = omp_get_max_threads()
    allocate(iter(nt), stat=ierr)
    if (ierr/=0) stop('allocation error')
    iter = 0


    best = -1e10
    last_improvement = 0

!$OMP PARALLEL PRIVATE(it, pop, i, targetFun)
      it = omp_get_thread_num()+1  ! thread number
      do
          iter(it) = iter(it) + 1

          ! Create a new population
          do i = 1,npop
             pop(1,i) = rand()
             pop(2,i) = rand()
          end do

          ! Evaluate target function  
          targetFun = popFun(pop)

          if (targetFun>best) then
          ! If this is the best population so far,
          !    then register this

!$OMP          CRITICAL
                  best_pop = pop
                  best = targetFun
                  print '(a,i0,a,i7,a,1p,e13.5)', '[',it,'] iteration ',sum(iter),' Best score until now: ',TargetFun
                  last_improvement = sum(iter) ! remember the global iteration count for the last time an improved population was found
!$OMP          END CRITICAL
          end if

          ! Done when further improvement takes too long
          if (last_improvement < sum(iter) - really_long) exit
      end do
!$OMP END PARALLEL

    ! Report the best population found
    targetFun = popFun(best_pop)
    print '(a,1p,e13.5)', 'Best score found: ',targetFun
    print '(a,1p,e13.5)', '    best population found:'
    do i = 1,npop
       print '(1p,10(a,e13.5))', '    (',best_pop(1,i),',',best_pop(2,i),')'
    end do

end program  optimize_population

Программа нуждается в целевой функции, предоставляемой модулем Population, ниже:

module Population
integer, parameter :: npop  = 20, ncenter = 3
integer, parameter :: dp = kind(1d0)
real(kind=dp)      :: center(2,ncenter)
contains

    subroutine initTarget()
    implicit none
       integer :: i
       do i = 1,ncenter
          center(1,i) = rand()
          center(2,i) = rand()
       end do
       print '(a,i0,a)', &
          'Looking for a population of ',npop,' points in the unit square,'
       print '(a,i0,a)', &
          'equally spread out in space, but clustered around the points'
       print '(1p,10(a,e13.5))', &
          '    (',center(1,1),',',center(2,1), &
          ('),    (',center(1,i),',',center(2,i), i=2,ncenter-1), &
          ')    and (',center(1,ncenter),',',center(2,ncenter),')'
    end subroutine initTarget


    function popFun(pop) result(targetFun)
    implicit none
        real(kind=dp), intent(in) :: pop(:,:)
        real(kind=dp) :: targetFun

        integer :: i,j
        real(kind=dp) :: sum_center, sum_dist

        sum_dist = 0
        sum_center = 0
        do i = 1,npop
           do j = i+1,npop
              sum_dist   = sum_dist + (pop(1,i)-pop(1,j))**2 + (pop(2,i)-pop(2,j))**2
           end do
           do j = 1,ncenter
              sum_center = sum_center + (pop(1,i)-center(1,j))**2 + (pop(2,i)-center(2,j))**2
           end do
        end do

        targetFun = sum_dist - sum_center
    end function popFun

end module Population
...