module de_opt implicit none private integer, parameter :: prec = selected_real_kind(8), max_iter = 1000, step = 10 real(prec), parameter :: c = 0.7, r = 0.1 public :: prec, initialize, optimize contains subroutine initialize(f, pop, cost, low, up) ! generates initial "populations" implicit none interface function f(x) result(res) import real(prec), dimension(:), intent(IN) :: x real(prec) :: res end function f end interface real(kind=prec), dimension(:,:), intent(OUT) :: pop real(prec), dimension(:), intent(OUT) :: cost real(prec), intent(IN) :: low, up integer :: pop_size, i pop_size = size(pop,2) if (size(cost) /= pop_size) stop call random_number(pop) pop = low + (up - low)*pop do i = 1, pop_size cost(i) = f(pop(:,i)) end do end subroutine initialize subroutine optimize(f, pop, cost) ! the main optimization algorithm implicit none interface function f(x) result(res) import real(prec), dimension(:), intent(IN) :: x real(prec) :: res end function f end interface real(prec), dimension(:,:), intent(INOUT) :: pop real(prec), dimension(:), intent(INOUT) :: cost integer :: n, pop_size, allocstat, iter, i, idx1, idx2, idx3 real(KIND=prec), dimension(:,:), allocatable :: new_pop real(KIND=prec), dimension(size(pop,1)) :: trial real(KIND=prec) :: score logical, dimension(size(pop,1)) :: mask n = size(pop,1) pop_size = size(pop,2) write (*,'(A,I2,/,A,I4,/,A,F5.3,/,A,F5.3)') & 'Dimension: ', n, & 'Population size: ', pop_size, & 'c=', c, & 'r=', r if (size(cost) /= pop_size) stop 'Error in population size' allocate(new_pop(n,pop_size), STAT=allocstat) if (allocstat /= 0) stop 'Failed to allocate memory' write (*,'(A,I4,G14.6)') & 'Initial values, cost', 0, minval(cost) do iter = 1, max_iter do i = 1, pop_size call triple(i, pop_size, idx1, idx2, idx3) mask = rnd(n) < r mask(idx(n)) = .true. where (mask) trial = pop(:,idx3) + c*(pop(:,idx1)-pop(:,idx2)) elsewhere trial = pop(:,i) end where score = f(trial) if (score < cost(i)) then new_pop(:,i) = trial cost(i) = score else new_pop(:,i) = pop(:,i) end if end do pop = new_pop if (mod(iter,step) == 0) then write (*,'(A,I4,G14.6)') 'Iteration, minimum: ', & iter, minval(cost) end if end do end subroutine optimize subroutine triple(i, n, i1, i2, i3) ! generates three indices implicit none integer, intent(IN) :: i, n integer, intent(OUT) :: i1, i2, i3 do i1 = idx(n) if (i1 /= i) exit end do do i2 = idx(n) if (i2 /= i .and. i2 /= i1) exit end do do i3 = idx(n) if (all(i3 /= (/i,i1,i2/))) exit end do end subroutine triple function rnd(n) ! returns an array (size n) of random numbers between [0,1] implicit none integer, intent(IN) :: n real(KIND=prec), dimension(n) :: rnd call random_number(rnd) end function rnd integer function idx(n) ! returns one random integer between 1 and n implicit none integer, intent(IN) :: n real(KIND=prec) :: x call random_number(x) idx = n*x + 1 end function idx end module de_opt