tem_balance_sparta Subroutine

  subroutine tem_balance_sparta(weight, myPart, nParts, comm, myElems, offset, &
    &                          sparta )
    ! ---------------------------------------------------------------------------
    !> Sorted list of weights corresponding to treeID order
    real(kind=rk),intent(in)  :: weight(:)
    integer,intent(in) :: myPart !< Rank of the calling process
    !> Number of procs the distribution should span
    integer, intent(in) :: nParts
    !> MPI Communicator
    integer, intent(in) :: comm
    !> number of elements
    integer, intent(inout) :: myElems
    !> Array of offsets with the size nParts. Offset index starts at 0.
    !! This Array needs to be allocate and deallocated outside
    integer(kind=long_k), intent(out) :: offset
    ! Count variables that state which rank gets how many elements
    ! *_count(rank0, ...)
    ! Right now the size is nParts but that will be changed in
    ! the near future to avoid O(p) allocations
    type( tem_sparta_type ), intent(inout) :: sparta
    ! ---------------------------------------------------------------------------
    integer :: iErr  ! MPI error variable
    integer :: iElem, iProc
    integer(kind=long_k) :: myElems_long
    real(kind=rk) :: w_sum, w_opt
    real(kind=rk) :: send, recv ! Send and receive buffers for MPI calls
    ! boundary values of the elements in which we search for splitters
    real(kind=rk) :: lower_boundary, upper_boundary
    ! local prefix sum array of myElems
    real(kind=rk), allocatable :: presum(:)
    integer :: rmin, rmax, lb, ub, left_off, mid
    real(kind=rk) :: opt_split, wsplit
    integer :: send_count(0:nParts-1)
    ! ---------------------------------------------------------------------------

    write(logUnit(5),*) "Balance by SpartA algorithm."

    send_count = 0

    ! Allocate Array for Prefix sum
    allocate(presum(myElems))
    ! Prefix sum over local weights. later on we will look for the splitter in
    ! this prefix sum
    presum(1) = weight(1)
    do iElem = 2,myElems
      presum(iElem) = presum(iElem-1) + weight(iElem)
    end do
    send = presum(myElems)
    ! sum up global total weight
    call MPI_ALLREDUCE(send, recv, 1, rk_mpi, mpi_sum, comm, iErr)

    w_sum = recv
    ! Calculate global optimum
    w_opt = w_sum / dble(nParts)

    ! Global prefix sum for weights
    call MPI_EXSCAN(send, recv, 1, rk_mpi, mpi_sum, comm, iErr)

    ! initialize splitter search
    lower_boundary=recv
    if (myPart == 0) lower_boundary = 0
    upper_boundary = lower_boundary + presum(myElems)

    rmin = max(floor(lower_boundary/w_opt),0)
    rmax = min(ceiling(upper_boundary / w_opt),nParts-1)

    ! Do splitter search
    left_off = 1
    do iProc = rmin,rmax
      lb = left_off
      ub = myelems
      opt_split = (iProc+1)*w_opt
      if (iProc*w_opt < upper_boundary) then
      do
        mid = (lb+ub)/2
        wsplit = presum(mid) + lower_boundary
        if (wsplit .feq. opt_split) exit
        if (wsplit < opt_split) then
           lb = mid
        else
           ub = mid
        end if
        ! exit if a single element was found, need to do this
        if (lb >= ub-1) exit
        !                    here, to have mid and wsplit set.
      end do
      if (ABS(wsplit - opt_split) > ABS(wsplit - opt_split - weight(mid))) then
        mid = mid - 1 ! return 0 if the splitter is left of the lower boundary
      else
        if (mid+1 <= myElems) then
          if (ABS(wsplit - opt_split)                                          &
            & > ABS(wsplit - opt_split + weight(mid+1))) then
            mid = mid + 1 ! return myElems at most
          end if
        else
          if (opt_split > upper_boundary) mid = myElems
        end if
      end if
      send_count(iProc) = mid - left_off + 1
      left_off = mid + 1
      end if
    end do
    ! finished splitter search. Communciate results.
    ! Each process needs to know how many elements to receive from which process
    call tem_set_sparta( sparta, comm, nParts, send_count )

    ! Calculate myElems and offset -----------------------------------
    ! total number of my elements after exchanging elements
    myElems = sparta%new_size
    myElems_long = int(myElems, kind=long_k)
    call mpi_exscan(myElems_long, offset, 1, long_k_mpi, mpi_sum, comm, ierr)
    if (myPart == 0) offset = 0
    ! write(*,"(3(A,I0))") 'myPart ', myPart, ' nElems: ', myElems, ' offset: ', offset
    ! Calculate myElems and offset -----------------------------------

  end subroutine tem_balance_sparta