pulp_ve.cpp

#include <mpi.h>
#include <omp.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

#include "comms.h"
#include "pulp_data.h"
#include "pulp_util.h"
#include "pulp_ve.h"
#include "util.h"
#include "xtrapulp.h"

//#define X 1.0
//#define Y 0.25
#define CUT_CHANGE 0.95
#define BAL_CHANGE 0.95
#define BAL_CUTOFF 1.05

extern int procid, nprocs;
extern int seed;
extern bool verbose, debug, verify;
extern float X, Y;

int pulp_ve(dist_graph_t* g, mpi_data_t* comm, queue_data_t* q,
            pulp_data_t* pulp, uint64_t outer_iter, uint64_t balance_iter,
            uint64_t refine_iter, double vert_balance, double edge_balance) {
  if (debug) {
    printf("Task %d pulp_ve() start\n", procid);
  }
  double elt = 0.0;
  if (verbose) {
    MPI_Barrier(MPI_COMM_WORLD);
    elt = omp_get_wtime();
  }

  q->queue_size = 0;
  q->next_size = 0;
  q->send_size = 0;

  for (int32_t i = 0; i < nprocs; ++i) comm->sendcounts_temp[i] = 0;

  pulp->max_v = 0.0;
  pulp->max_e = 0.0;
  for (int p = 0; p < pulp->num_parts; ++p) {
    if ((double)pulp->part_sizes[p] / pulp->avg_size > pulp->max_v)
      pulp->max_v = (double)pulp->part_sizes[p] / pulp->avg_size;
    if ((double)pulp->part_edge_sizes[p] / pulp->avg_edge_size > pulp->max_e)
      pulp->max_e = (double)pulp->part_edge_sizes[p] / pulp->avg_edge_size;
  }

  double tot_iter = (double)(outer_iter * (refine_iter + balance_iter));
  double cur_iter = 0.0;
  double multiplier = (double)nprocs * ((X - Y) * (cur_iter / tot_iter) + Y);
  // double running_bal = pulp->max_e;
  // double running_cut = (double)pulp->cut_size;

  uint64_t num_swapped_1 = 0;
  uint64_t num_swapped_2 = 0;
  comm->global_queue_size = 1;
#pragma omp parallel default(shared)
  {
    thread_queue_t tq;
    thread_comm_t tc;
    thread_pulp_t tp;
    init_thread_queue(&tq);
    init_thread_comm(&tc);
    init_thread_pulp(&tp, pulp);
    xs1024star_t xs;
    xs1024star_seed((uint64_t)(seed + omp_get_thread_num()), &xs);

    for (uint64_t cur_outer_iter = 0; cur_outer_iter < outer_iter;
         ++cur_outer_iter) {
#pragma omp single
      {
        // if (procid == 0) printf("EVAL ec VE EB
        // ------------------------------\n");  part_eval(g, pulp);
        update_pulp_data(g, pulp);
        num_swapped_1 = 0;
      }

      for (uint64_t cur_bal_iter = 0; cur_bal_iter < balance_iter;
           ++cur_bal_iter) {
        for (int32_t p = 0; p < pulp->num_parts; ++p) {
          tp.part_weights[p] =
              vert_balance * pulp->avg_size / (double)pulp->part_sizes[p] - 1.0;
          tp.part_edge_weights[p] = pulp->max_e * pulp->avg_edge_size /
                                        (double)pulp->part_edge_sizes[p] -
                                    1.0;
          if (tp.part_weights[p] < 0.0) tp.part_weights[p] = 0.0;
          if (tp.part_edge_weights[p] < 0.0) tp.part_edge_weights[p] = 0.0;
        }

#pragma omp for schedule(guided) reduction(+ : num_swapped_1) nowait
        for (uint64_t vert_index = 0; vert_index < g->n_local; ++vert_index) {
          int32_t part = pulp->local_parts[vert_index];
          for (int32_t p = 0; p < pulp->num_parts; ++p) tp.part_counts[p] = 0.0;

          uint64_t out_degree = out_degree(g, vert_index);
          uint64_t* outs = out_vertices(g, vert_index);
          for (uint64_t j = 0; j < out_degree; ++j) {
            uint64_t out_index = outs[j];
            int32_t part_out = pulp->local_parts[out_index];
            tp.part_counts[part_out] += 1.0;
          }

          int32_t max_part = part;
          double max_val = 0.0;
          uint64_t num_max = 0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            if (tp.part_weights[p] > 0.0 && tp.part_edge_weights[p] > 0.0)
              tp.part_counts[p] *=
                  (tp.part_weights[p] * tp.part_edge_weights[p] *
                   pulp->weight_exponent_e);
            else
              tp.part_counts[p] = 0.0;

            if (tp.part_counts[p] == max_val) {
              tp.part_counts[num_max++] = (double)p;
            } else if (tp.part_counts[p] > max_val) {
              max_val = tp.part_counts[p];
              max_part = p;
              num_max = 0;
              tp.part_counts[num_max++] = (double)p;
            }
          }

          if (num_max > 1)
            max_part =
                (int32_t)tp.part_counts[(xs1024star_next(&xs) % num_max)];

          if (max_part != part) {
            ++num_swapped_1;
#pragma omp atomic
            --pulp->part_size_changes[part];
#pragma omp atomic
            ++pulp->part_size_changes[max_part];
#pragma omp atomic
            pulp->part_edge_size_changes[part] -= (int64_t)out_degree;
#pragma omp atomic
            pulp->part_edge_size_changes[max_part] += (int64_t)out_degree;

            tp.part_weights[part] =
                vert_balance * pulp->avg_size /
                    ((double)pulp->part_sizes[part] +
                     multiplier * (double)pulp->part_size_changes[part]) -
                1.0;
            tp.part_weights[max_part] =
                vert_balance * pulp->avg_size /
                    ((double)pulp->part_sizes[max_part] +
                     multiplier * (double)pulp->part_size_changes[max_part]) -
                1.0;

            tp.part_edge_weights[part] =
                pulp->max_e * pulp->avg_edge_size /
                    ((double)pulp->part_edge_sizes[part] +
                     multiplier * (double)pulp->part_edge_size_changes[part]) -
                1.0;
            tp.part_edge_weights[max_part] =
                pulp->max_e * pulp->avg_edge_size /
                    ((double)pulp->part_edge_sizes[max_part] +
                     multiplier *
                         (double)pulp->part_edge_size_changes[max_part]) -
                1.0;

            if (tp.part_weights[part] < 0.0) tp.part_weights[part] = 0.0;
            if (tp.part_weights[max_part] < 0.0)
              tp.part_weights[max_part] = 0.0;

            if (tp.part_edge_weights[part] < 0.0)
              tp.part_edge_weights[part] = 0.0;
            if (tp.part_edge_weights[max_part] < 0.0)
              tp.part_edge_weights[max_part] = 0.0;

            pulp->local_parts[vert_index] = max_part;
            add_vid_to_send(&tq, q, vert_index);
            // add_vid_to_queue(&tq, q, vert_index);
          }
        }

        empty_send(&tq, q);
// empty_queue(&tq, q);
#pragma omp barrier

        for (int32_t i = 0; i < nprocs; ++i) tc.sendcounts_thread[i] = 0;

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_sendcounts_thread(g, &tc, vert_index);
        }

        for (int32_t i = 0; i < nprocs; ++i) {
#pragma omp atomic
          comm->sendcounts_temp[i] += tc.sendcounts_thread[i];

          tc.sendcounts_thread[i] = 0;
        }
#pragma omp barrier

#pragma omp single
        { init_sendbuf_vid_data(comm); }

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_vid_data_queues(g, &tc, comm, vert_index,
                                 pulp->local_parts[vert_index]);
        }

        empty_vid_data(&tc, comm);
#pragma omp barrier

#pragma omp single
        { exchange_vert_data(g, comm, q); }  // end single

#pragma omp for
        for (uint64_t i = 0; i < comm->total_recv; ++i) {
          uint64_t index = get_value(g->map, comm->recvbuf_vert[i]);
          pulp->local_parts[index] = comm->recvbuf_data[i];
        }

#pragma omp single
        {
          clear_recvbuf_vid_data(comm);

          MPI_Allreduce(MPI_IN_PLACE, pulp->part_size_changes, pulp->num_parts,
                        MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);
          MPI_Allreduce(MPI_IN_PLACE, pulp->part_edge_size_changes,
                        pulp->num_parts, MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);

          pulp->max_v = 0.0;
          pulp->max_e = 0.0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            pulp->part_sizes[p] += pulp->part_size_changes[p];
            pulp->part_edge_sizes[p] += pulp->part_edge_size_changes[p];
            pulp->part_size_changes[p] = 0;
            pulp->part_edge_size_changes[p] = 0;

            if ((double)pulp->part_sizes[p] / pulp->avg_size > pulp->max_v)
              pulp->max_v = (double)pulp->part_sizes[p] / pulp->avg_size;
            if ((double)pulp->part_edge_sizes[p] / pulp->avg_edge_size >
                pulp->max_e)
              pulp->max_e =
                  (double)pulp->part_edge_sizes[p] / pulp->avg_edge_size;
          }
          pulp->weight_exponent_e *= pulp->max_e / edge_balance;
          if (pulp->max_e < edge_balance) {
            pulp->max_e = edge_balance;
            pulp->weight_exponent_e = 1.0;
          }

          cur_iter += 1.0;
          multiplier = (double)nprocs * ((X - Y) * (cur_iter / tot_iter) + Y);

          if (debug)
            printf("Task %d num_swapped_1 %lu \n", procid, num_swapped_1);
          num_swapped_1 = 0;
        }

      }  // end balance loop

#pragma omp single
      {
        // if (procid == 0) printf("EVAL ec VE ER
        // ------------------------------\n");  part_eval(g, pulp);
        update_pulp_data(g, pulp);
        num_swapped_2 = 0;
      }

      for (uint64_t cur_ref_iter = 0; cur_ref_iter < refine_iter;
           ++cur_ref_iter) {
#pragma omp for schedule(guided) reduction(+ : num_swapped_2) nowait
        for (uint64_t vert_index = 0; vert_index < g->n_local; ++vert_index) {
          int32_t part = pulp->local_parts[vert_index];
          for (int32_t p = 0; p < pulp->num_parts; ++p) tp.part_counts[p] = 0.0;

          uint64_t out_degree = out_degree(g, vert_index);
          uint64_t* outs = out_vertices(g, vert_index);
          for (uint64_t j = 0; j < out_degree; ++j) {
            uint64_t out_index = outs[j];
            int32_t part_out = pulp->local_parts[out_index];
            tp.part_counts[part_out] += 1.0;
          }

          int32_t max_part = part;
          double max_val = 0.0;
          uint64_t num_max = 0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            if (tp.part_counts[p] == max_val) {
              tp.part_counts[num_max++] = (double)p;
            } else if (tp.part_counts[p] > max_val) {
              max_val = tp.part_counts[p];
              max_part = p;
              num_max = 0;
              tp.part_counts[num_max++] = (double)p;
            }
          }

          if (num_max > 1)
            max_part =
                (int32_t)tp.part_counts[(xs1024star_next(&xs) % num_max)];

          if (max_part != part) {
            int64_t new_size = (int64_t)pulp->avg_size;
            int64_t new_edge_size = (int64_t)pulp->avg_edge_size;

            new_size =
                pulp->part_size_changes[max_part] + 1 < 0
                    ? pulp->part_sizes[max_part] +
                          pulp->part_size_changes[max_part] + 1
                    : (int64_t)((double)pulp->part_sizes[max_part] +
                                multiplier *
                                    (double)pulp->part_size_changes[max_part] +
                                1.0);

            new_edge_size =
                pulp->part_edge_size_changes[max_part] + out_degree < 0
                    ? pulp->part_edge_sizes[max_part] +
                          pulp->part_edge_size_changes[max_part] + out_degree
                    : (int64_t)(
                          (double)pulp->part_edge_sizes[max_part] +
                          multiplier *
                              (double)pulp->part_edge_size_changes[max_part] +
                          (double)out_degree);

            if (new_size < (int64_t)(pulp->avg_size * vert_balance) &&
                new_edge_size < (int64_t)(pulp->avg_edge_size * pulp->max_e)) {
              ++num_swapped_2;
#pragma omp atomic
              --pulp->part_size_changes[part];
#pragma omp atomic
              ++pulp->part_size_changes[max_part];
#pragma omp atomic
              pulp->part_edge_size_changes[part] -= (int64_t)out_degree;
#pragma omp atomic
              pulp->part_edge_size_changes[max_part] += (int64_t)out_degree;

              pulp->local_parts[vert_index] = max_part;
              add_vid_to_send(&tq, q, vert_index);
              // add_vid_to_queue(&tq, q, vert_index);
            }
          }
        }

        empty_send(&tq, q);
// empty_queue(&tq, q);
#pragma omp barrier

        for (int32_t i = 0; i < nprocs; ++i) tc.sendcounts_thread[i] = 0;

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_sendcounts_thread(g, &tc, vert_index);
        }

        for (int32_t i = 0; i < nprocs; ++i) {
#pragma omp atomic
          comm->sendcounts_temp[i] += tc.sendcounts_thread[i];

          tc.sendcounts_thread[i] = 0;
        }
#pragma omp barrier

#pragma omp single
        { init_sendbuf_vid_data(comm); }

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_vid_data_queues(g, &tc, comm, vert_index,
                                 pulp->local_parts[vert_index]);
        }

        empty_vid_data(&tc, comm);
#pragma omp barrier

#pragma omp single
        { exchange_vert_data(g, comm, q); }  // end single

#pragma omp for
        for (uint64_t i = 0; i < comm->total_recv; ++i) {
          uint64_t index = get_value(g->map, comm->recvbuf_vert[i]);
          pulp->local_parts[index] = comm->recvbuf_data[i];
        }

#pragma omp single
        {
          clear_recvbuf_vid_data(comm);

          MPI_Allreduce(MPI_IN_PLACE, pulp->part_size_changes, pulp->num_parts,
                        MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);
          MPI_Allreduce(MPI_IN_PLACE, pulp->part_edge_size_changes,
                        pulp->num_parts, MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);

          // pulp->max_e = 0.0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            pulp->part_sizes[p] += pulp->part_size_changes[p];
            pulp->part_edge_sizes[p] += pulp->part_edge_size_changes[p];
            pulp->part_size_changes[p] = 0;
            pulp->part_edge_size_changes[p] = 0;

            // if ((double)pulp->part_edge_sizes[p] / pulp->avg_edge_size >
            // pulp->max_e)
            // pulp->max_e = (double)pulp->part_edge_sizes[p] /
            // pulp->avg_edge_size;
          }

          cur_iter += 1.0;
          multiplier = (double)nprocs * ((X - Y) * (cur_iter / tot_iter) + Y);

          // int32_t* tmp = pulp->local_parts;
          // pulp->local_parts = pulp->local_parts_next;
          // pulp->local_parts_next = tmp;
          if (debug)
            printf("Task %d num_swapped_2 %lu \n", procid, num_swapped_2);
          num_swapped_2 = 0;
        }

      }  // end refine iter

      /*if (cur_outer_iter + 1 == outer_iter)
      {
    #pragma omp single
    {
        update_pulp_data(g, pulp);

        if ( (pulp->max_e > edge_balance*BAL_CUTOFF &&
              pulp->max_e*BAL_CHANGE < running_bal) ||
             (double)pulp->cut_size < running_cut*CUT_CHANGE )
        {
          outer_iter += 1;
          tot_iter += (double)(balance_iter + refine_iter);
          running_cut = (double)pulp->cut_size;
          running_bal = pulp->max_e;
          multiplier = (double)nprocs*( (X - Y)*(cur_iter/tot_iter) + Y );
        }
    }
      }*/

    }  // end outer loop

    clear_thread_queue(&tq);
    clear_thread_comm(&tc);
    clear_thread_pulp(&tp);
  }  // end parallel

  // part_eval(g, pulp);

  if (verbose) {
    elt = omp_get_wtime() - elt;
    printf("Task %d, pulp_ve() time %9.6f (s)\n", procid, elt);
  }
  if (debug) {
    printf("Task %d pulp_ve() success\n", procid);
  }

  return 0;
}

int pulp_ve_weighted(dist_graph_t* g, mpi_data_t* comm, queue_data_t* q,
                     pulp_data_t* pulp, uint64_t outer_iter,
                     uint64_t balance_iter, uint64_t refine_iter,
                     double vert_balance, double edge_balance) {
  if (debug) {
    printf("Task %d pulp_ve_weighted() start\n", procid);
  }
  double elt = 0.0;
  if (verbose) {
    MPI_Barrier(MPI_COMM_WORLD);
    elt = omp_get_wtime();
  }

  bool has_vwgts = (g->vertex_weights != NULL);
  bool has_ewgts = (g->edge_weights != NULL);

  q->queue_size = 0;
  q->next_size = 0;
  q->send_size = 0;

  for (int32_t i = 0; i < nprocs; ++i) comm->sendcounts_temp[i] = 0;

  pulp->max_v = 0.0;
  pulp->max_e = 0.0;
  for (int p = 0; p < pulp->num_parts; ++p) {
    if ((double)pulp->part_sizes[p] / pulp->avg_size > pulp->max_v)
      pulp->max_v = (double)pulp->part_sizes[p] / pulp->avg_size;
    if ((double)pulp->part_edge_sizes[p] / pulp->avg_edge_size > pulp->max_e)
      pulp->max_e = (double)pulp->part_edge_sizes[p] / pulp->avg_edge_size;
  }

  double tot_iter = (double)(outer_iter * (refine_iter + balance_iter));
  double cur_iter = 0.0;
  double multiplier = (double)nprocs * ((X - Y) * (cur_iter / tot_iter) + Y);
  // double running_bal = pulp->max_e;
  // double running_cut = (double)pulp->cut_size;

  uint64_t num_swapped_1 = 0;
  uint64_t num_swapped_2 = 0;
  comm->global_queue_size = 1;
#pragma omp parallel default(shared)
  {
    thread_queue_t tq;
    thread_comm_t tc;
    thread_pulp_t tp;
    init_thread_queue(&tq);
    init_thread_comm(&tc);
    init_thread_pulp(&tp, pulp);
    xs1024star_t xs;
    xs1024star_seed((uint64_t)(seed + omp_get_thread_num()), &xs);

    for (uint64_t cur_outer_iter = 0; cur_outer_iter < outer_iter;
         ++cur_outer_iter) {
#pragma omp single
      {
        // if (procid == 0) printf("EVAL ec VE EB
        // ------------------------------\n");  part_eval(g, pulp);
        update_pulp_data_weighted(g, pulp);
        num_swapped_1 = 0;
      }

      for (uint64_t cur_bal_iter = 0; cur_bal_iter < balance_iter;
           ++cur_bal_iter) {
        for (int32_t p = 0; p < pulp->num_parts; ++p) {
          tp.part_weights[p] =
              vert_balance * pulp->avg_size / (double)pulp->part_sizes[p] - 1.0;
          tp.part_edge_weights[p] = pulp->max_e * pulp->avg_edge_size /
                                        (double)pulp->part_edge_sizes[p] -
                                    1.0;
          if (tp.part_weights[p] < 0.0) tp.part_weights[p] = 0.0;
          if (tp.part_edge_weights[p] < 0.0) tp.part_edge_weights[p] = 0.0;
        }

#pragma omp for schedule(guided) reduction(+ : num_swapped_1) nowait
        for (uint64_t vert_index = 0; vert_index < g->n_local; ++vert_index) {
          int32_t part = pulp->local_parts[vert_index];
          int32_t vert_weight = 1;
          if (has_vwgts) vert_weight = g->vertex_weights[vert_index];

          for (int32_t p = 0; p < pulp->num_parts; ++p) tp.part_counts[p] = 0.0;

          uint64_t out_degree = out_degree(g, vert_index);
          uint64_t* outs = out_vertices(g, vert_index);
          int32_t* weights = out_weights(g, vert_index);
          for (uint64_t j = 0; j < out_degree; ++j) {
            uint64_t out_index = outs[j];
            int32_t part_out = pulp->local_parts[out_index];
            double weight_out = 1.0;
            if (has_ewgts) weight_out = (double)weights[j];
            tp.part_counts[part_out] += weight_out;
          }

          int32_t max_part = part;
          double max_val = 0.0;
          uint64_t num_max = 0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            if (tp.part_weights[p] > 0.0 && tp.part_edge_weights[p] > 0.0)
              tp.part_counts[p] *=
                  (tp.part_weights[p] * tp.part_edge_weights[p] *
                   pulp->weight_exponent_e);
            else
              tp.part_counts[p] = 0.0;

            if (tp.part_counts[p] == max_val) {
              tp.part_counts[num_max++] = (double)p;
            } else if (tp.part_counts[p] > max_val) {
              max_val = tp.part_counts[p];
              max_part = p;
              num_max = 0;
              tp.part_counts[num_max++] = (double)p;
            }
          }

          if (num_max > 1)
            max_part =
                (int32_t)tp.part_counts[(xs1024star_next(&xs) % num_max)];

          if (max_part != part) {
            ++num_swapped_1;
#pragma omp atomic
            pulp->part_size_changes[part] -= vert_weight;
#pragma omp atomic
            pulp->part_size_changes[max_part] += vert_weight;
#pragma omp atomic
            pulp->part_edge_size_changes[part] -= (int64_t)out_degree;
#pragma omp atomic
            pulp->part_edge_size_changes[max_part] += (int64_t)out_degree;

            tp.part_weights[part] =
                vert_balance * pulp->avg_size /
                    ((double)pulp->part_sizes[part] +
                     multiplier * (double)pulp->part_size_changes[part]) -
                1.0;
            tp.part_weights[max_part] =
                vert_balance * pulp->avg_size /
                    ((double)pulp->part_sizes[max_part] +
                     multiplier * (double)pulp->part_size_changes[max_part]) -
                1.0;

            tp.part_edge_weights[part] =
                pulp->max_e * pulp->avg_edge_size /
                    ((double)pulp->part_edge_sizes[part] +
                     multiplier * (double)pulp->part_edge_size_changes[part]) -
                1.0;
            tp.part_edge_weights[max_part] =
                pulp->max_e * pulp->avg_edge_size /
                    ((double)pulp->part_edge_sizes[max_part] +
                     multiplier *
                         (double)pulp->part_edge_size_changes[max_part]) -
                1.0;

            if (tp.part_weights[part] < 0.0) tp.part_weights[part] = 0.0;
            if (tp.part_weights[max_part] < 0.0)
              tp.part_weights[max_part] = 0.0;

            if (tp.part_edge_weights[part] < 0.0)
              tp.part_edge_weights[part] = 0.0;
            if (tp.part_edge_weights[max_part] < 0.0)
              tp.part_edge_weights[max_part] = 0.0;

            pulp->local_parts[vert_index] = max_part;
            add_vid_to_send(&tq, q, vert_index);
            // add_vid_to_queue(&tq, q, vert_index);
          }
        }

        empty_send(&tq, q);
// empty_queue(&tq, q);
#pragma omp barrier

        for (int32_t i = 0; i < nprocs; ++i) tc.sendcounts_thread[i] = 0;

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_sendcounts_thread(g, &tc, vert_index);
        }

        for (int32_t i = 0; i < nprocs; ++i) {
#pragma omp atomic
          comm->sendcounts_temp[i] += tc.sendcounts_thread[i];

          tc.sendcounts_thread[i] = 0;
        }
#pragma omp barrier

#pragma omp single
        { init_sendbuf_vid_data(comm); }

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_vid_data_queues(g, &tc, comm, vert_index,
                                 pulp->local_parts[vert_index]);
        }

        empty_vid_data(&tc, comm);
#pragma omp barrier

#pragma omp single
        { exchange_vert_data(g, comm, q); }  // end single

#pragma omp for
        for (uint64_t i = 0; i < comm->total_recv; ++i) {
          uint64_t index = get_value(g->map, comm->recvbuf_vert[i]);
          pulp->local_parts[index] = comm->recvbuf_data[i];
        }

#pragma omp single
        {
          clear_recvbuf_vid_data(comm);

          MPI_Allreduce(MPI_IN_PLACE, pulp->part_size_changes, pulp->num_parts,
                        MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);
          MPI_Allreduce(MPI_IN_PLACE, pulp->part_edge_size_changes,
                        pulp->num_parts, MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);

          pulp->max_v = 0.0;
          pulp->max_e = 0.0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            pulp->part_sizes[p] += pulp->part_size_changes[p];
            pulp->part_edge_sizes[p] += pulp->part_edge_size_changes[p];
            pulp->part_size_changes[p] = 0;
            pulp->part_edge_size_changes[p] = 0;

            if ((double)pulp->part_sizes[p] / pulp->avg_size > pulp->max_v)
              pulp->max_v = (double)pulp->part_sizes[p] / pulp->avg_size;
            if ((double)pulp->part_edge_sizes[p] / pulp->avg_edge_size >
                pulp->max_e)
              pulp->max_e =
                  (double)pulp->part_edge_sizes[p] / pulp->avg_edge_size;
          }
          pulp->weight_exponent_e *= pulp->max_e / edge_balance;
          if (pulp->max_e < edge_balance) {
            pulp->max_e = edge_balance;
            pulp->weight_exponent_e = 1.0;
          }

          cur_iter += 1.0;
          multiplier = (double)nprocs * ((X - Y) * (cur_iter / tot_iter) + Y);

          if (debug)
            printf("Task %d num_swapped_1 %lu \n", procid, num_swapped_1);
          num_swapped_1 = 0;
        }

      }  // end balance loop

#pragma omp single
      {
        // if (procid == 0) printf("EVAL ec VE ER
        // ------------------------------\n");  part_eval(g, pulp);
        update_pulp_data_weighted(g, pulp);
        num_swapped_2 = 0;
      }

      for (uint64_t cur_ref_iter = 0; cur_ref_iter < refine_iter;
           ++cur_ref_iter) {
#pragma omp for schedule(guided) reduction(+ : num_swapped_2) nowait
        for (uint64_t vert_index = 0; vert_index < g->n_local; ++vert_index) {
          int32_t part = pulp->local_parts[vert_index];
          int32_t vert_weight = 1;
          if (has_vwgts) vert_weight = g->vertex_weights[vert_index];

          for (int32_t p = 0; p < pulp->num_parts; ++p) tp.part_counts[p] = 0.0;

          uint64_t out_degree = out_degree(g, vert_index);
          uint64_t* outs = out_vertices(g, vert_index);
          int32_t* weights = out_weights(g, vert_index);
          for (uint64_t j = 0; j < out_degree; ++j) {
            uint64_t out_index = outs[j];
            int32_t part_out = pulp->local_parts[out_index];
            double weight_out = 1.0;
            if (has_ewgts) weight_out = (double)weights[j];
            tp.part_counts[part_out] += weight_out;
          }

          int32_t max_part = part;
          double max_val = 0.0;
          uint64_t num_max = 0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            if (tp.part_counts[p] == max_val) {
              tp.part_counts[num_max++] = (double)p;
            } else if (tp.part_counts[p] > max_val) {
              max_val = tp.part_counts[p];
              max_part = p;
              num_max = 0;
              tp.part_counts[num_max++] = (double)p;
            }
          }

          if (num_max > 1)
            max_part =
                (int32_t)tp.part_counts[(xs1024star_next(&xs) % num_max)];

          if (max_part != part) {
            int64_t new_size = (int64_t)pulp->avg_size;
            int64_t new_edge_size = (int64_t)pulp->avg_edge_size;

            new_size =
                pulp->part_size_changes[max_part] + (int64_t)vert_weight < 0
                    ? pulp->part_sizes[max_part] +
                          pulp->part_size_changes[max_part] +
                          (int64_t)vert_weight
                    : (int64_t)((double)pulp->part_sizes[max_part] +
                                multiplier *
                                    (double)pulp->part_size_changes[max_part] +
                                (double)vert_weight);

            new_edge_size =
                pulp->part_edge_size_changes[max_part] + out_degree < 0
                    ? pulp->part_edge_sizes[max_part] +
                          pulp->part_edge_size_changes[max_part] + out_degree
                    : (int64_t)(
                          (double)pulp->part_edge_sizes[max_part] +
                          multiplier *
                              (double)pulp->part_edge_size_changes[max_part] +
                          (double)out_degree);

            if (new_size < (int64_t)(pulp->avg_size * vert_balance) &&
                new_edge_size < (int64_t)(pulp->avg_edge_size * pulp->max_e)) {
              ++num_swapped_2;
#pragma omp atomic
              pulp->part_size_changes[part] -= vert_weight;
#pragma omp atomic
              pulp->part_size_changes[max_part] += vert_weight;
#pragma omp atomic
              pulp->part_edge_size_changes[part] -= (int64_t)out_degree;
#pragma omp atomic
              pulp->part_edge_size_changes[max_part] += (int64_t)out_degree;

              pulp->local_parts[vert_index] = max_part;
              add_vid_to_send(&tq, q, vert_index);
              // add_vid_to_queue(&tq, q, vert_index);
            }
          }
        }

        empty_send(&tq, q);
// empty_queue(&tq, q);
#pragma omp barrier

        for (int32_t i = 0; i < nprocs; ++i) tc.sendcounts_thread[i] = 0;

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_sendcounts_thread(g, &tc, vert_index);
        }

        for (int32_t i = 0; i < nprocs; ++i) {
#pragma omp atomic
          comm->sendcounts_temp[i] += tc.sendcounts_thread[i];

          tc.sendcounts_thread[i] = 0;
        }
#pragma omp barrier

#pragma omp single
        { init_sendbuf_vid_data(comm); }

#pragma omp for schedule(guided) nowait
        for (uint64_t i = 0; i < q->send_size; ++i) {
          uint64_t vert_index = q->queue_send[i];
          update_vid_data_queues(g, &tc, comm, vert_index,
                                 pulp->local_parts[vert_index]);
        }

        empty_vid_data(&tc, comm);
#pragma omp barrier

#pragma omp single
        { exchange_vert_data(g, comm, q); }  // end single

#pragma omp for
        for (uint64_t i = 0; i < comm->total_recv; ++i) {
          uint64_t index = get_value(g->map, comm->recvbuf_vert[i]);
          pulp->local_parts[index] = comm->recvbuf_data[i];
        }

#pragma omp single
        {
          clear_recvbuf_vid_data(comm);

          MPI_Allreduce(MPI_IN_PLACE, pulp->part_size_changes, pulp->num_parts,
                        MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);
          MPI_Allreduce(MPI_IN_PLACE, pulp->part_edge_size_changes,
                        pulp->num_parts, MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD);

          // pulp->max_e = 0.0;
          for (int32_t p = 0; p < pulp->num_parts; ++p) {
            pulp->part_sizes[p] += pulp->part_size_changes[p];
            pulp->part_edge_sizes[p] += pulp->part_edge_size_changes[p];
            pulp->part_size_changes[p] = 0;
            pulp->part_edge_size_changes[p] = 0;

            // if ((double)pulp->part_edge_sizes[p] / pulp->avg_edge_size >
            // pulp->max_e)
            // pulp->max_e = (double)pulp->part_edge_sizes[p] /
            // pulp->avg_edge_size;
          }

          cur_iter += 1.0;
          multiplier = (double)nprocs * ((X - Y) * (cur_iter / tot_iter) + Y);

          // int32_t* tmp = pulp->local_parts;
          // pulp->local_parts = pulp->local_parts_next;
          // pulp->local_parts_next = tmp;
          if (debug)
            printf("Task %d num_swapped_2 %lu \n", procid, num_swapped_2);
          num_swapped_2 = 0;
        }

      }  // end refine iter

      /*if (cur_outer_iter + 1 == outer_iter)
      {
    #pragma omp single
    {
        update_pulp_data(g, pulp);

        if ( (pulp->max_e > edge_balance*BAL_CUTOFF &&
              pulp->max_e*BAL_CHANGE < running_bal) ||
             (double)pulp->cut_size < running_cut*CUT_CHANGE )
        {
          outer_iter += 1;
          tot_iter += (double)(balance_iter + refine_iter);
          running_cut = (double)pulp->cut_size;
          running_bal = pulp->max_e;
          multiplier = (double)nprocs*( (X - Y)*(cur_iter/tot_iter) + Y );
        }
    }
      }*/

    }  // end outer loop

    clear_thread_queue(&tq);
    clear_thread_comm(&tc);
    clear_thread_pulp(&tp);
  }  // end parallel

  // part_eval(g, pulp);

  if (verbose) {
    elt = omp_get_wtime() - elt;
    printf("Task %d, pulp_ve_weighted() time %9.6f (s)\n", procid, elt);
  }
  if (debug) {
    printf("Task %d pulp_ve_weighted() success\n", procid);
  }

  return 0;
}