         vid_t* u_nodes = hornet.vertex(u).neighbor_ptr();
         vid_t* v_nodes = hornet.vertex(v).neighbor_ptr();
         
-        int work_per_thread = std::max(total_work/threads_per_union, (unsigned long)1);
-        int diag_id;
-        diag_id = thread_id*work_per_thread;
-
+        int work_per_thread = total_work/threads_per_union;
+        int remainder_work = total_work % threads_per_union;
+        int diag_id, next_diag_id;
+        diag_id = thread_union_id*work_per_thread + std::min(thread_union_id, remainder_work);
+        next_diag_id = (thread_union_id+1)*work_per_thread + std::min(thread_union_id+1, remainder_work);
+        //printf("u=%d, v=%d, diag_id=%d, union_id=%d, total_work=%d, work_per_thread=%d, remainder_work=%d\n",
+        //        u, v, diag_id, thread_union_id, total_work, work_per_thread, remainder_work);
         vid_t low_ui, low_vi, high_vi, high_ui, ui_curr, vi_curr;
-        if ((diag_id > 0) && (diag_id < total_work-1)) {
+        if (diag_id > 0 && diag_id < total_work) {
             // For the binary search, we are figuring out the initial poT of search.
             if (diag_id < u_len) {
                 low_ui = diag_id-1;


    // Initialize queue sizes to zero
    for (auto i = 0; i < MAX_ADJ_UNIONS_BINS; i++)
        hd_queue_info().queue_sizes[i] = 0;

    // Phase 1: determine and bin all edges based on edge neighbor properties
    // First, count the number to avoid creating excessive queues
    timer::Timer<timer::DEVICE> TM(5);
    TM.start();
    forAllEdgesSrcDst(hornet, bin_edges {hd_queue_info, true}, load_balancing);
    TM.stop();
    TM.print("counting queues");
    TM.reset();
    hd_queue_info.sync();

    for (auto i = 0; i < MAX_ADJ_UNIONS_BINS; i++)

    int       id = blockIdx.x * blockDim.x + threadIdx.x;
    int queue_id = id / threads_per_union;
    int thread_id = id % threads_per_union;
    int thread_union_id = threadIdx.x % threads_per_union;
    int block_local_id = threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    int queue_stride = stride / threads_per_union;
