init_AliasTable_t Subroutine

    subroutine init_AliasTable_t(this, rank_with_info, arr)
        class(AliasTable_t), intent(inout) :: this
        integer, intent(in) :: rank_with_info
            !! The **intra-node** rank that contains the weights to be used
            !! all other arr of all other ranks are ignored (and can be allocated with size 0).
        real(dp), intent(in) :: arr(:)

        ! as this is shared memory, only node-root has to do this
        if (iProcIndex_intra == rank_with_info) then
        block
            integer :: i, j, cV, cU
            integer(int64) :: arrSize
            integer, allocatable :: overfull(:), underfull(:)

            arrSize = size(arr, kind=int64)
            ! initialize the probabilities
            this%bias%ptr = arr / sum(arr) * arrSize

            ! indices of subarrays
            allocate(overfull(arrSize), underfull(arrSize))

            cV = 0
            cU = 0
            do i = 1, int(arrSize)
                call assignLabel(i, overfull, cV, underfull, cU)
            end do
            ! we now labeled each entry

            ! it is more efficient to start with the largest biases
            ! -> reverse overfull
            overfull(1:cV) = overfull(cV:1:-1)
            do
                if ((cV == 0) .or. (cU == 0)) then
                    exit
                end if
                ! pick one overfull and one underfull index
                i = overfull(cV)
                j = underfull(cU)
                ! set the alias of the underfull to be the other
                this%alias%ptr(j) = i
                ! correct the bias
                this%bias%ptr(i) = this%bias%ptr(i) + this%bias%ptr(j) - 1.0_dp

                ! unmark j
                cU = cU - 1
                ! reassign i based on the new bias
                cV = cV - 1
                call assignLabel(i, overfull, cV, underfull, cU)
            end do

            ! make sure we do not leave anything unfilled
            call roundTo1(overfull, cV)
            call roundTo1(underfull, cU)

        end block
        end if

    contains

        subroutine assignLabel(i, overfull, cV, underfull, cU)
            integer, intent(in) :: i
            integer, intent(inout) :: overfull(:), cV, underfull(:), cU

            if (this%bias%ptr(i) > 1) then
                cV = cV + 1
                overfull(cV) = i
            else
                cU = cU + 1
                underfull(cU) = i
            end if
        end subroutine assignLabel

        subroutine roundTo1(labels, cI)
            integer, intent(in) :: labels(:)
            integer, intent(in) :: cI
            integer :: i

            ! if, due to floating point errors, one of the categories is not empty, empty it
            ! (error is negligible then)
            if (cI > 0) then
                do i = 1, cI
                    this%bias%ptr(labels(i)) = 1.0_dp
                    this%alias%ptr(labels(i)) = labels(i)
                end do
            end if

        end subroutine roundTo1

    end subroutine init_AliasTable_t