perturb_and_evaluate Subroutine

   subroutine perturb_and_evaluate(me,x,vm,xp,fp,nfcnev,ier,first)

      class(simulated_annealing_type),intent(inout) :: me
      real(wp),dimension(:),intent(in)  :: x       !! input optimization variable vector to perturb
      real(wp),dimension(:),intent(in)  :: vm      !! step length vector
      real(wp),dimension(:),intent(out) :: xp      !! the perturbed `x` value
      real(wp),intent(out)              :: fp      !! the value of the user function at `xp`
      integer,intent(inout)             :: nfcnev  !! total number of function evaluations
      integer,intent(inout)             :: ier     !! status output code
      logical,intent(in),optional       :: first   !! to use the input `x` the first time

      integer :: i          !! counter
      integer :: istat      !! user function status code
      logical :: first_try  !! local copy of `first`
      integer :: n_inputs   !! number of inputs to send to the user function for parallel evaluation
      logical :: reallocate !! whether to reallocate `xp_mat` for parallel evaluation
      real(wp),dimension(:,:),allocatable :: xp_mat !! array of `xp` vectors for parallel evaluation

      if (present(first)) then
         first_try = first
      else
         first_try = .false.
      end if

      do  ! we must return at least one that doesn't fail (or we exceed the max number of function evaluations)

         if (me%parallel_mode) then
            ! parallel mode: generate and send multiple x vectors
            call me%n_inputs_to_send(n_inputs)  !! get the inputs to send to the user function for parallel evaluation
            n_inputs = max(1, n_inputs)  !! there must be at least one
            ! do we need to reallocate xp_mat? (if n_inputs has changed):
            reallocate = .true.
            if (allocated(xp_mat)) then
               if (size(xp_mat,1) == me%n .and. size(xp_mat,2) == n_inputs) then
                  reallocate = .false.
               end if
            end if
            if (reallocate) then
               if (allocated(xp_mat)) deallocate(xp_mat)
               allocate(xp_mat(me%n, n_inputs))
            end if
            do i = 1, n_inputs
               xp_mat(:,i) = get_xp()  ! get a perturbed `x` value for each input
            end do
            call me%fcn_parallel_input(xp_mat)  !! send the perturbed `x` values to the user function for parallel evaluation
            call me%fcn_parallel_output(xp, fp, istat)  !! get the function values and status code back from the user function for parallel evaluation
         else
            ! serial mode: generate and send a single x vector
            xp = get_xp()  ! single funciton evaluation
            call me%fcn(xp, fp, istat)  ! evaluate the function with the trial point xp and return as fp.
         end if

         if (istat==0 .and. (me%ireport == 1 .or. me%ireport == 3)) then ! report this value to the user
            call me%report(xp, fp, istat=1) ! convert f back to user's sign if necessary
         end if

         nfcnev = nfcnev + 1  ! function eval counter (note: for parallel runs,
                              ! only count one evaluation per returned values.
                              ! others can still be running)

         ! if too many function evaluations occur, terminate the algorithm.
         if (nfcnev > me%maxevl) then
            if (me%iprint>0) then
               write(me%iunit, '(A)') ' too many function evaluations; consider'
               write(me%iunit, '(A)') ' increasing maxevl or eps, or decreasing'
               write(me%iunit, '(A)') ' nt or rt. these results are likely to be'
               write(me%iunit, '(A)') ' poor.'
            end if
            ier = 1
            return
         end if

         select case (istat)
          case(-2)
            ! user stop
            write(me%iunit, '(A)') ' user stopped in function.'
            ier = 4
            exit
          case(-1)
            ! try another one until we get a valid evaluation
            cycle
          case default
            ! continue
         end select

         exit ! finished

      end do

   contains

      function get_xp() result(xp)
         !! get a perturbed `x` value
         real(wp),dimension(me%n) :: xp     !! the perturbed `x` value
         real(wp)                 :: lower  !! lower bound to use for random interval
         real(wp)                 :: upper  !! upper bound to use for random interval
         integer                  :: i      !! counter

         if (first_try) then
            if (me%use_initial_guess) then
               ! use the initial guess
               ! [note if this evauation fails, the subsequent ones
               !  are perturbed from this one]
               xp = x
               first_try = .false.
            else
               ! a random point in the bounds:
               ! [if it fails, a new random one is tried next time]
               do i = 1, me%n
                  xp(i) = me%perturb_variable(i, x(i), me%distribution_mode(i), &
                                              me%lb(i), me%ub(i))
               end do
            end if
         else
            ! perturb all of them using per-variable distributions:
            do i = 1, me%n
               lower = max( me%lb(i), x(i) - vm(i) )
               upper = min( me%ub(i), x(i) + vm(i) )
               xp(i) = me%perturb_variable(i, x(i), me%distribution_mode(i), &
                                           lower, upper)
            end do
         end if
      end function get_xp

   end subroutine perturb_and_evaluate