<|start_filename|>include/label.hpp<|end_filename|>
#ifndef _label_hpp_
#define _label_hpp_
#include <cuda_runtime.h>
#include "edge.hpp"
#include "io.hpp"
#ifndef CSP
#define CSP 1
#endif

__device__ float *A_dd, *B_dd, *Theta_dd;

__global__ void assign(float *A, float *B, float *Theta) {
  A_dd = A;
  B_dd = B;
  Theta_dd = Theta;
}

class label {
 public:
  int vid;
 public:
  __device__ label() {
    return;
  }

  __device__ virtual int get_vid() {
    return vid;
  }

  __device__ label
  expand(edge e) {
    return label();
  }

  __device__ static bool pass_constraints_check() {
    return false;
  }
  __device__ static bool dominate(label a, label b) {
    return false;
  }

  ~label() {
  }
};

#if CSP == 2
class c_label {
public	:
    int vid;
    float d;
    float delay;
    int hop;
    int father_id;
    c_label *pre;

public	:
    __device__ c_label() {
        return;
    }

    c_label(int vid) : vid(vid) {
        d = 0;
        delay = 0;
        hop = 0;
        father_id = -1;
        pre = NULL;
    }

    __device__ c_label(int vid, float d, float delay, int hop) :
    vid(vid), d(d), delay(delay), hop(hop) {
        pre = NULL;

        return;
    }

    __device__ static void copy(c_label *a, c_label *b) {
        a->vid = b->vid;
        a->d = b->d;
        a->delay = b->delay;
        a->hop = b->hop;
        a->father_id = b->father_id;
        a->pre = b->pre;
    }

    __device__ int get_vid() {
        return vid;
    }

    template<typename Edge_t>
    __device__ void expand(int next_v, Edge_t e, c_label *son) {
        son->vid = next_v;
        son->d = this->d + e.len;
        son->delay = this->delay + e.delay;
        son->hop = this->hop + 1;
        son->father_id = this->vid;
        son->pre = this->pre;
    }

    __device__ static bool dominate(c_label a, c_label b) {
        if (a.d <= b.d && a.delay <= b.delay && a.hop <= b.hop) {
            return true;
        }
        return false;
    }
};
#elif CSP == 3
class c_label {
public	:
    int vid;
    float d;
    float delay;
    int hop;
    float maxcost;
    int father_id;
    c_label *pre;

public	:
    __device__ c_label() {
        return;
    }

    c_label(int vid) : vid(vid) {
        d = 0;
        delay = 0;
        hop = 0;
        maxcost = 0;
        father_id = -1;
        pre = NULL;
    }

    __device__ c_label(int vid, float d, float delay, int hop, float maxcost) :
    vid(vid), d(d), delay(delay), hop(hop), maxcost(maxcost) {
        pre = NULL;

        return;
    }

    __device__ static void copy(c_label *a, c_label *b) {
        a->vid = b->vid;
        a->d = b->d;
        a->delay = b->delay;
        a->hop = b->hop;
        a->maxcost = b->maxcost;
        a->father_id = b->father_id;
        a->pre = b->pre;
    }

    __device__ int get_vid() {
        return vid;
    }

    template<typename Edge_t>
    __device__ void expand(int next_v, Edge_t e, c_label *son) {
        son->vid = next_v;
        son->d = this->d + e.len;
        son->delay = this->delay + e.delay;
        son->hop = this->hop + 1;
        son->maxcost = max(this->maxcost, e.len);
        son->father_id = this->vid;
        son->pre = this->pre;
    }

    __device__ static bool dominate(c_label a, c_label b) {
        if (a.d <= b.d && a.delay <= b.delay && a.hop <= b.hop && a.maxcost <= b.maxcost) {
            return true;
        }
        return false;
    }
};
#elif CSP == 1
class c_label {
 public    :
  int vid;
  float d;
  float delay;
//  int father_id;
//  c_label *pre;

 public    :
  __device__ c_label() {
    return;
  }

  c_label(int vid) : vid(vid) {
    d = 0;
    delay = 0;
//    father_id = -1;
//    pre = NULL;
  }

  __device__ c_label(int vid, float d, float delay) :
      vid(vid), d(d), delay(delay) {
//    pre = NULL;

    return;
  }
  __device__ static bool equal(c_label *a, c_label *b) {
    return a->vid == b->vid && a->d == b->d && a->delay == b->delay;
  }

  __device__ static void copy(c_label *a, c_label *b) {
    a->vid = b->vid;
    a->d = b->d;
    a->delay = b->delay;
//    a->father_id = b->father_id;
//    a->pre = b->pre;
  }

  __device__ int get_vid() {
    return vid;
  }

  template<typename Edge_t>
  __device__ void expand(int next_v, Edge_t e, c_label *son) {
    son->vid = next_v;
    son->d = this->d + e.len;
    son->delay = this->delay + e.delay;
//    son->father_id = this->vid;
//    son->pre = this->pre;
//    printf("%d->%d %lf %lf %lf %lf %lf %lf\n", this->vid, next_v, e.delay, this->delay, son->delay, this->d, e.len, son->d);
//    printf("%d->%d %lf %lf\n", this->vid, next_v, son->d, son->delay);
  }

  template<typename Edge_t>
  void expand_cpu(int next_v, Edge_t e, c_label *son) {
    son->vid = next_v;
    son->d = this->d * e.len;
    son->delay = this->delay * e.delay;
//    son->father_id = this->vid;
//    son->pre = this->pre;
//    printf("%d->%d %lf %lf %lf\n", this->vid + 1, next_v + 1, e.delay, this->delay, son->delay);
  }

  __device__ static bool dominate(c_label a, c_label b) {
    return a.d <= b.d && a.delay <= b.delay;
  }
  __device__ static bool dominate_vid(c_label a, c_label b) {
    if (a.vid != b.vid) {
      return false;
    }
    return a.d <= b.d && a.delay <= b.delay;
  }

  static bool dominate_cpu(c_label a, c_label b) {
    return a.d <= b.d && a.delay <= b.delay;
  }
};
#elif CSP == 4 //2-OSP
class c_label {
public	:
    int vid;
    float d;
    float delay;
    int father_id;
    c_label *pre;

public	:
    __device__ c_label() {
        return;
    }

    c_label(int vid) : vid(vid) {
        d = 0;
        delay = 0;
        father_id = -1;
        pre = NULL;
    }

    __device__ c_label(int vid, float d, float delay) :
    vid(vid), d(d), delay(delay) {
        pre = NULL;

        return;
    }

    __device__ static void copy(c_label *a, c_label *b) {
        a->vid = b->vid;
        a->d = b->d;
        a->delay = b->delay;
        a->father_id = b->father_id;
        a->pre = b->pre;
    }

    __device__ int get_vid() {
        return vid;
    }

    template<typename Edge_t>
    __device__ void expand(int next_v, Edge_t e, c_label *son) {
        son->vid = next_v;
        son->d = this->d + e.len;
        son->delay = this->delay + e.delay;
        son->father_id = this->vid;
        son->pre = this->pre;
    }

    __device__ static bool dominate(c_label a, c_label b) {
        if (a.d <= b.d && a.delay <= b.delay) {
            return true;
        }
        return false;
    }
};
#elif CSP == 5 // social
class c_label {
 public    :
  int vid;
  float d;
  float delay;
//  int father_id;
//  c_label *pre;

 public    :
  __device__ c_label() {
    return;
  }

  c_label(int vid) : vid(vid) {
    d = 1;
    delay = 1;
//    father_id = -1;
//    pre = NULL;
  }

  __device__ c_label(int vid, float d, float delay) :
      vid(vid), d(d), delay(delay) {
//    pre = NULL;

    return;
  }
  __device__ static bool equal(c_label *a, c_label *b) {
    return a->vid == b->vid && a->d == b->d && a->delay == b->delay;
  }

  __device__ static void copy(c_label *a, c_label *b) {
    a->vid = b->vid;
    a->d = b->d;
    a->delay = b->delay;
//    a->father_id = b->father_id;
//    a->pre = b->pre;
  }

  __device__ int get_vid() {
    return vid;
  }

  template<typename Edge_t>
  __device__ void expand(int next_v, Edge_t e, c_label *son) {
    son->vid = next_v;
    son->d = this->d * e.len;
    son->delay = this->delay * e.delay;
//    son->father_id = this->vid;
//    son->pre = this->pre;
//    printf("%d->%d %lf %lf %lf %lf %lf %lf\n", this->vid, next_v, e.delay, this->delay, son->delay, this->d, e.len, son->d);
//    printf("%d->%d %lf %lf\n", this->vid, next_v, son->d, son->delay);
  }

  template<typename Edge_t>
  void expand_cpu(int next_v, Edge_t e, c_label *son) {
    son->vid = next_v;
    son->d = this->d * e.len;
    son->delay = this->delay * e.delay;
//    son->father_id = this->vid;
//    son->pre = this->pre;
//    printf("%d->%d %lf %lf %lf\n", this->vid + 1, next_v + 1, e.delay, this->delay, son->delay);
  }

  __device__ static bool dominate(c_label a, c_label b) {
    return a.d >= b.d && a.delay >= b.delay;
  }
  __device__ static bool dominate_vid(c_label a, c_label b) {
    if (a.vid != b.vid) {
      return false;
    }
    return a.d >= b.d && a.delay >= b.delay;
  }

  static bool dominate_cpu(c_label a, c_label b) {
    return a.d >= b.d && a.delay >= b.delay;
  }
};
#else // network qos
class c_label {
 public    :
  int vid;
  float d;
  float delay;
//  int father_id;
//  c_label *pre;

 public    :
  __device__ c_label() {
    return;
  }

  c_label(int vid) : vid(vid) {
    d = 0;
    delay = 21474836;
//    father_id = -1;
//    pre = NULL;
  }

  __device__ c_label(int vid, float d, float delay) :
      vid(vid), d(d), delay(delay) {
//    pre = NULL;

    return;
  }
  __device__ static bool equal(c_label *a, c_label *b) {
    return a->vid == b->vid && a->d == b->d && a->delay == b->delay;
  }

  __device__ static void copy(c_label *a, c_label *b) {
    a->vid = b->vid;
    a->d = b->d;
    a->delay = b->delay;
//    a->father_id = b->father_id;
//    a->pre = b->pre;
  }
  __device__ int get_vid() {
    return vid;
  }

  template<typename Edge_t>
  __device__ void expand(int next_v, Edge_t e, c_label *son) {
    son->vid = next_v;
    son->d = this->d + e.len;
    son->delay = min(this->delay, e.delay);
//    son->father_id = this->vid;
//    son->pre = this->pre;
//    printf("%d->%d %lf %lf %lf %lf %lf %lf\n", this->vid, next_v, e.delay, this->delay, son->delay, this->d, e.len, son->d);
//    printf("%d->%d %lf %lf\n", this->vid, next_v, son->d, son->delay);
  }

  template<typename Edge_t>
  void expand_cpu(int next_v, Edge_t e, c_label *son) {
    son->vid = next_v;
    son->d = this->d + e.len;
    son->delay = min(this->delay, e.delay);
//    son->father_id = this->vid;
//    son->pre = this->pre;
//    printf("%d->%d %lf %lf %lf\n", this->vid + 1, next_v + 1, e.delay, this->delay, son->delay);
  }

  __device__ static bool dominate(c_label a, c_label b) {
    if (a.d <= b.d && a.delay >= b.delay) {
        return true;
    }
    return false;
  }
  __device__ static bool dominate_vid(c_label a, c_label b) {
    if (a.vid != b.vid) {
      return false;
    }
    if (a.d <= b.d && a.delay >= b.delay) {
      return true;
    }
    return false;
  }

  static bool dominate_cpu(c_label a, c_label b) {
    return a.d <= b.d && a.delay >= b.delay;
  }
};
#endif

#endif

<|start_filename|>include/graph.hpp<|end_filename|>
#ifndef _graph_hpp_
#define _graph_hpp_
#include <cstdio>
#include <vector>
#include <cassert>
#include <algorithm>
#include <queue>
#include <set>
#include <list>
#include <limits>
#include <iostream>
#include <string>
#include <map>

#include "gpu_graph.hpp"
#include "edge.hpp"

using namespace std;
template<typename Edge_t>
class Graph;

template<typename Edge_t>
class gpuGraph;

class Edge {
 public:
  int id;
  int src, snk;
  float len, delay;
 public:
  Edge() {
    ;
  }
  Edge(int id) :
      id(id) {
    ;
  }
  Edge(int gg, int v, int I);

  int assign(
      int id,
      int src,
      int snk,
      float len = 0.0,
      float delay = 0.0) {
    this->id = id;
    this->src = src;
    this->snk = snk;
    this->len = len;
    this->delay = delay;

    return 0;
  }

  int getOwner() const {
    return src;
  }
  int operator()() const;
  int getVersion() const;
  bool operator<(const Edge &o) const;
};

template<typename Graph_t>
class Vertex {
 public:
  int id;
  vector<Edge *> edges;
  vector<int> next;
  vector<Edge *> edges_r;
  vector<int> pre;
 public:
  Vertex() {
    ;
  }
  Vertex(int id) :
      id(id) {
    ;
  }
  int assign(Edge *e, Graph_t *g) {
    this->edges.push_back(e);
    this->next.push_back(e->snk);
    Vertex *son = &(g->v_l[e->snk]);
    son->edges_r.push_back(e);
    son->pre.push_back(this->id);

    return 0;
  }
};

template<typename Edge_t>
class Graph {
 public:

  using MyV = Vertex<Graph<Edge_t> >;
  vector<MyV> v_l;
  vector<Edge *> e_l;
  int N, M;

  // CSR part
  // set as public for simpler manipulatign
  int *_count_row; // CSR pointer
  int *_eid_col; // CSR edge id
  int *_vid_col; //CSR tail v id
  int *e_sink; // eid to sink
  int *e_csr;  // eid to csr pos
  float *_len; // CSR edge len
  float *_delay; // CSR edge delay

  Edge_t *el;

  gpuGraph<Edge_t> *gg; // gpu graph on device
  gpuGraph<Edge_t> *gg_h; // gpu graph on host

 public:
  Graph() {
    // data allocation is outside... since public
    ;
  }
  Graph(int N, int M) :
      N(N), M(M), e_l(M) {
    v_l.reserve(N);
    e_l.reserve(2 * M);
    for (int i = 0; i < N; i++) {
      MyV *v = new MyV(i);
      v_l.push_back(*v);
    }
    for (int i = 0; i < M; i++) {
      e_l[i] = new Edge(i);
    }
  }

  int adj_to_csr();
  int make_gpu_graph();
};

template<typename Edge_t>
int Graph<Edge_t>::adj_to_csr() {
  // allocate memory for CSR
  _count_row = (int *) malloc(sizeof(int) * (N + 1));
  _eid_col = (int *) malloc(sizeof(int) * M);
  _vid_col = (int *) malloc(sizeof(int) * M);
  e_sink = (int *) malloc(sizeof(int) * M);
  e_csr = (int *) malloc(sizeof(int) * M);
  _len = (float *) malloc(sizeof(float) * M);
  _delay = (float *) malloc(sizeof(float) * M);

  el = (Edge_t *) malloc(sizeof(Edge_t) * M);
  _count_row[0] = 0;
  for (int i = 0; i < N; i++) {
    _count_row[i + 1] = _count_row[i] + v_l[i].edges.size();
    if (v_l[i].edges.empty()) {
      continue;
    }
    Edge *e;
    for (int j = _count_row[i], t = 0; j < _count_row[i + 1]; j++, t++) {
      e = v_l[i].edges[t];
      e_csr[e->id] = j;
      e_sink[e->id] = e->src;
      _eid_col[j] = e->id;
      _vid_col[j] = e->snk;
      _len[j] = e->len;
      _delay[j] = e->delay;
    }
  }

  return 0;
}

template<typename Edge_t>
int Graph<Edge_t>::make_gpu_graph() {
  gg_h = new gpuGraph<Edge_t>();
  gg_h->N = N;
  cudaMalloc(&(gg_h->_count_row_d), sizeof(int) * (N + 1));
  cudaMalloc(&(gg_h->_eid_col_d), sizeof(int) * M);
  cudaMalloc(&(gg_h->_vid_col_d), sizeof(int) * M);
  cudaMalloc(&(gg_h->_sink_d), sizeof(int) * M);
  cudaMalloc(&(gg_h->_csr_d), sizeof(int) * M);
  cudaMalloc(&(gg_h->_len_data_d), sizeof(float) * M);
  cudaMalloc(&(gg_h->_delay_data_d), sizeof(float) * M);
  cudaMalloc(&(gg_h->el_d), sizeof(Edge_t) * M);

  cudaMemcpy(gg_h->_count_row_d, _count_row, sizeof(int) * (N + 1), cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->_eid_col_d, _eid_col, sizeof(int) * M, cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->_vid_col_d, _vid_col, sizeof(int) * M, cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->_sink_d, e_sink, sizeof(int) * M, cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->_csr_d, e_csr, sizeof(int) * M, cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->_len_data_d, _len, sizeof(float) * M, cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->_delay_data_d, _delay, sizeof(float) * M, cudaMemcpyHostToDevice);
  cudaMemcpy(gg_h->el_d, el, sizeof(Edge_t) * M, cudaMemcpyHostToDevice);

  cudaMalloc(&(this->gg), sizeof(gpuGraph<Edge_t>));
  cudaMemcpy(this->gg, this->gg_h, sizeof(gpuGraph<Edge_t>), cudaMemcpyHostToDevice);

  return 0;
}

#endif


<|start_filename|>include/edge.hpp<|end_filename|>
#ifndef _edge_hpp_
#define _edge_hpp_

template<typename Edge_t>
class Graph;

class edge {
 public:
  __device__ edge() {
  }

  ~edge() {
  }
};

class c_edge {
 public:
  float len;
  float delay;
 public:
  c_edge() {
  }
  c_edge(float l, float d) :
      len(l), delay(d) {
  }
};

#endif


<|start_filename|>algorithm/dominance.cuh<|end_filename|>
#ifndef _dominance_cuh_
#define _dominance_cuh_

#include <cuda_runtime.h>
#include "../include/label.hpp"
#include "../include/constraint.hpp"
#include "../include/edge.hpp"
#include "../include/graph.hpp"
#include "../include/util.hpp"

#define SearchSpace 128
#define GlobalSpaceFactor 2
#define _K 4

template<typename Label_t>
__device__ void second_level_pruning(
    Label_t *cur_v,
    Label_t *v_dominance_list,
    int &dominance_list_len,
    Label_t *v_dynamic_list,
    int &dynamic_list_len,
    Label_t *local_buffer,
    int *expand_buffer_size,
    int *mutex) {
  __shared__ int flag;
  // -1 is dominated
  // 1 dominate others
  // 0 non, append

  flag = 0;
  int vid = cur_v->vid;
  // each thread now work on 1 comparison
  int lid = threadIdx.y * blockDim.x + threadIdx.x;
  if (dominance_list_len <= SearchSpace) {
    while (lid < dominance_list_len && lid < SearchSpace && flag != -1) {
      if (Label_t::dominate(v_dominance_list[lid], *cur_v)) {
        flag = -1;
      }
      lid += blockDim.x * blockDim.y;
    }
  } else {
    while (lid < SearchSpace && flag != -1) {
      if (Label_t::dominate(v_dominance_list[lid], *cur_v)) {
        flag = -1;
      }
      lid += blockDim.x * blockDim.y;
    }
    lid = threadIdx.y * blockDim.x + threadIdx.x;
    while (lid < dynamic_list_len && flag != -1) {
      if (Label_t::dominate_vid(v_dynamic_list[lid], *cur_v)) {
        flag = -1;
      }
      lid += blockDim.x * blockDim.y;
    }
  }
  __syncthreads();
  if (flag == -1) {
    cur_v->vid = -1;
    return;
  }
  lid = threadIdx.y * blockDim.x + threadIdx.x;
  if (dominance_list_len <= SearchSpace) {
    while (lid < dominance_list_len && lid <= SearchSpace && flag != 1) {
      // analysis dominance situation
      if (Label_t::dominate(*cur_v, v_dominance_list[lid])) {
        bool leave = true;
        while (leave) {
          if (atomicCAS(mutex + vid, 0, 1) == 0) {
            if (Label_t::dominate(*cur_v, v_dominance_list[lid])) {
              Label_t::copy(&v_dominance_list[lid], cur_v);
//            cur_v->pre = &v_dominance_list[lid];
              flag = 1;
            }
            __threadfence_block();
            atomicExch(mutex + vid, 0);
            leave = false;
          }
        }
      }
      lid += blockDim.x * blockDim.y;
    }
  } else {
    while (lid < SearchSpace && flag != 1) {
      // analysis dominance situation
      if (Label_t::dominate(*cur_v, v_dominance_list[lid])) {
        bool leave = true;
        while (leave) {
          if (atomicCAS(mutex + vid, 0, 1) == 0) {
            if (Label_t::dominate(*cur_v, v_dominance_list[lid])) {
              Label_t::copy(&v_dominance_list[lid], cur_v);
              flag = 1;
            }
            __threadfence_block();
            atomicExch(mutex + vid, 0);
            leave = false;
          }
        }
      }
      lid += blockDim.x * blockDim.y;
    }
    __syncthreads();
    lid = threadIdx.y * blockDim.x + threadIdx.x;
    while (lid < dynamic_list_len && flag != 1) {
      // analysis dominance situation
      if (Label_t::dominate_vid(*cur_v, v_dynamic_list[lid])) {
        bool leave = true;
        while (leave) {
          if (atomicCAS(mutex + (vid / _K), 0, 1) == 0) {
            if (Label_t::dominate_vid(*cur_v, v_dynamic_list[lid])) {
              Label_t::copy(&v_dynamic_list[lid], cur_v);
              flag = 1;
            }
            __threadfence_block();
            atomicExch(mutex + (vid / _K), 0);
            leave = false;
          }
        }
      }
      lid += blockDim.x * blockDim.y;
    }
  }
  __syncthreads();
  if (threadIdx.x == 0 && threadIdx.y == 0) {
    if (flag == 0) {
      bool leave = true;
      while (leave) {
        if (dominance_list_len < SearchSpace) {
          if (atomicCAS(mutex + vid, 0, 1) == 0) {
            Label_t::copy(&v_dominance_list[dominance_list_len], cur_v);
            (dominance_list_len)++;
            __threadfence_block();
            atomicExch(mutex + vid, 0);
            leave = false;
          }
        } else {
          if (atomicCAS(mutex + (vid / _K), 0, 1) == 0) {
            Label_t::copy(&v_dynamic_list[dynamic_list_len], cur_v);
            (dynamic_list_len)++;
            (dominance_list_len)++;
            __threadfence_block();
            atomicExch(mutex + (vid / _K), 0);
            leave = false;
          }
        }
      }
    }
    // fill into bucket
    Label_t::copy(&local_buffer[*expand_buffer_size], cur_v);
    (*expand_buffer_size)++;
    cur_v->vid = -1;
  }
  __syncthreads();
}

template<typename Label_t>
__device__
void second_level_pruning_only(
    Label_t *cur_v,
    Label_t *v_dominance_list,
    int &dominance_list_len,
    Label_t *v_dynamic_list,
    int &dynamic_list_len,
    Label_t *local_buffer,
    int *expand_buffer_size) {
  __shared__ int flag;
  // -1 is dominated
  // 1 dominate others
  // 0 non, append
  flag = 0;
  int lid = threadIdx.y * blockDim.x + threadIdx.x;
  while (lid < dominance_list_len && lid < SearchSpace && flag != -1) {
    if (Label_t::dominate(v_dominance_list[lid], *cur_v) && !Label_t::equal(&v_dominance_list[lid], cur_v)) {
      flag = -1;
    }
    lid += blockDim.x * blockDim.y;
  }
  __syncthreads();
  if (flag == -1) {
    cur_v->vid = -1;
    return;
  }
  if (threadIdx.x == 0) {
    // fill into bucket
    Label_t::copy(&local_buffer[*expand_buffer_size], cur_v);
    (*expand_buffer_size)++;
    cur_v->vid = -1;
  }
  __syncthreads();
}

#endif


<|start_filename|>include/gpu_graph.hpp<|end_filename|>
#ifndef _gpu_graph_hpp_
#define _gpu_graph_hpp_
#include <cuda_runtime.h>

template<typename Edge_t>
class gpuGraph {
 public:
  int N;
  bool direct;

  int *_count_row_d; // CSR pointer
  int *_eid_col_d; // CSR edge id
  int *_vid_col_d; //CSR tail v id
  int *_sink_d; // eid to sink
  int *_csr_d; // eid to csr pos
  float *_len_data_d; // CSR edge len
  float *_delay_data_d; // CSR edge delay

  Edge_t *el_d;
 public:
  gpuGraph() {
    // cudaMalloc(gg_d);
    // cudaMemcpy(gg_d, g->gg);
  }
  ~gpuGraph() {
    //	cudaFree(csr_count_row_d);
    //	cudaFree(csr_eid_col_d);
    //	cudaFree(csr_vid_col_d);
    //	cudaFree(eid_cap_d);
    //	cudaFree(eid_sink_d);
    //	cudaFree(eid_csr_d);
    //	cudaFree(csr_len_data_d);
    //	cudaFree(csr_delay_data_d);
  }
};

#endif /* GPUGRAPH_H_ */


<|start_filename|>algorithm/cspp.cuh<|end_filename|>
#ifndef _cspp_cuh_
#define _cspp_cuh_

#include <cuda_runtime.h>
#include "../include/label.hpp"
#include "../include/constraint.hpp"
#include "../include/edge.hpp"
#include "../include/graph.hpp"
#include "../include/util.hpp"
#include "thrust/scan.h"
#include "thrust/execution_policy.h"
#include "thrust/device_ptr.h"
#include "bucket.cuh"
#include "../include/TimeMeasurer.hpp"
#include "../include/AR.hpp"
#include <math.h>

#include <cuda_profiler_api.h>

#define USE_AR

#define LRG_FLOAT 1000000
#define NUM_SM 56
#define ITER 1000
#ifdef USE_AR
#define AR_DEGREE 3
#define THETA 0.95
#else
#define MAX_LABELS_IN_FRONTIER 1000000
#endif
int v_per_block = (BLOCK_SIZE / 32) * VSIZE;

template<typename Label_t>
double calculate_p(double ar_n, int n, int max_labels) {
  // M / size = N - P - theta*P + AR(n)P/n
  double est = max_labels * 1.0 / sizeof(Label_t) / (ar_n / n - THETA);
  if (est < 0) {
    return (double) max_labels;
  }
  return min(est, (double) max_labels);
}

int max_frontier(int totalLabel) {
  size_t free, total;
  cudaMemGetInfo(&free, &total);
  return free / sizeof(c_label) - totalLabel;
}

// 2 Cons massive with check function
template<typename Label_t, typename Edge_t>
int RUN_CSPP(
    Graph<Edge_t> *g,
    int ntask,
    vector<Label_t> &init_state_l,
    vector<int> &sink_l,
    vector<int *> &path_l,
    vector<int> &size_of_result_l,
    vector<c_constraint<Label_t>> &C_l);

template<typename Label_t>
__global__ void ck_CSPP_path(int source, Label_t **dominance_list, int *list_len, int *path, int *mutex);

template<typename Label_t>
__global__ void ck_init(int *dominance_list_len, Label_t **dominance_list, Label_t source) {
  dominance_list_len[source.vid] = 1;
  Label_t::copy(dominance_list[source.vid], &source);
}

template<typename Label_t, typename Edge_t>
int RUN_CSPP(
    Graph<Edge_t> *g,
    int ntask,
    vector<Label_t> &init_state_l,
    vector<int> &sink_l,
    vector<int *> &path_l,
    vector<int> &size_of_result_l,
    vector<c_constraint<Label_t>> &C_l) {
  int N = g->N;
  TimeMeasurer t;
  cuda_err_chk(cudaGetLastError());

  Label_t *buffer;
  cudaMallocManaged(&buffer, sizeof(Label_t) * N * SearchSpace);
  Label_t *bucket;
  cudaMallocManaged(&bucket, sizeof(Label_t) * N * SearchSpace);
  cuda_err_chk(cudaGetLastError());

  int *scan_out, *status;
  cudaMallocManaged(&scan_out, sizeof(int) * N * SearchSpace);
  cudaMallocManaged(&status, sizeof(int) * N * SearchSpace);

  int *mutex;
  cudaMallocManaged(&mutex, sizeof(int) * N * SearchSpace);
  cuda_err_chk(cudaGetLastError());

  cudaMemset(mutex, 0, sizeof(int) * N * SearchSpace);
  Label_t **dominance_list_d;
  cudaMallocManaged(&dominance_list_d, sizeof(Label_t *) * N);
  Label_t **dynamic_list_d;
  cudaMallocManaged(&dynamic_list_d, sizeof(Label_t *) * N / _K);
  cudaDeviceSynchronize();
  cuda_err_chk(cudaGetLastError());

  int *dominance_list_len;
  cudaMallocManaged(&dominance_list_len, sizeof(int) * N);
  int *dynamic_list_len;
  cudaMallocManaged(&dynamic_list_len, sizeof(int) * N / _K);

  Label_t *dominance_list;
  Label_t *dynamic_list;
  vector<Label_t *> dominance_list_l;
  vector<Label_t *> dynamic_list_l;
  cudaMallocManaged(&dominance_list, sizeof(Label_t) * N * SearchSpace);
  cudaMallocManaged(&dynamic_list, sizeof(Label_t) * N * SearchSpace * GlobalSpaceFactor);

  for (int T = 0; T < N; T++) {
    dominance_list_l.push_back(&dominance_list[T * SearchSpace]);
    if (T % _K == 0)
      dynamic_list_l.push_back(&dynamic_list[T * SearchSpace * GlobalSpaceFactor]);
  }
  cudaMemcpy(dominance_list_d, &dominance_list_l[0], sizeof(Label_t *) * N, cudaMemcpyHostToDevice);
  cudaMemcpy(dynamic_list_d, &dynamic_list_l[0], sizeof(Label_t *) * N / _K, cudaMemcpyHostToDevice);
  cudaDeviceSynchronize();
  int expanded = 0;
  cuda_err_chk(cudaGetLastError());

  int DOP;
  std::vector<double> frontier_size;
  DOP = 100000;
  frontier_size.push_back(0);
  double *coefficients = (double *) malloc(AR_DEGREE * sizeof(double));
  int frontier_sum = 0;
  int max_bucket_size = 0;
  for (int i = 0; i < ntask; i++) {
    int total_label = 0;
    int MAX_LABELS_IN_FRONTIER = max_frontier(frontier_sum);

    DOP = MAX_LABELS_IN_FRONTIER / ITER;
    frontier_size.clear();
    int size_of_bucket = 1;
    int new_size_of_bucket = 0;
    cudaMemset(status, 0, sizeof(int) * N * SearchSpace);
    cudaMemset(bucket, 0, sizeof(Label_t) * N * SearchSpace);
    Label_t source = init_state_l[i];

    cudaMemcpy(bucket, &source, sizeof(Label_t), cudaMemcpyHostToDevice);
    cudaMemset(dominance_list, 0, sizeof(Label_t) * N * SearchSpace);
    cuda_err_chk(cudaGetLastError());

    cudaMemset(dynamic_list, 0, sizeof(Label_t) * N * SearchSpace * GlobalSpaceFactor);
    cudaMemset(dominance_list_len, 0, sizeof(int) * N);
    cudaDeviceSynchronize();
    cuda_err_chk(cudaGetLastError());
    cudaMemset(dynamic_list_len, 0, sizeof(int) * N / _K);
    ck_init<Label_t> <<< 1, 1 >>> (dominance_list_len, dominance_list_d, source);
    cudaDeviceSynchronize();
    cuda_err_chk(cudaGetLastError());
    thrust::device_ptr<int> status_ptr = thrust::device_pointer_cast(status);
    thrust::device_ptr<int> scan_out_ptr = thrust::device_pointer_cast(scan_out);
    frontier_sum = 0;
    double est = 0;
    bool f = false;
    int round = 0;
    cudaDeviceSynchronize();
    TimeMeasurer otherT;
    TimeMeasurer expandT;
    TimeMeasurer pruneT;
    TimeMeasurer t_ex;
    TimeMeasurer t_measure;

    t.tic();
    t_measure.resume();
    while (size_of_bucket != 0) {
      round++;
      cudaDeviceSynchronize();
      cuda_err_chk(cudaGetLastError());
      dim3 block(32, BLOCK_SIZE / 32);
      int num_of_block, SS;
      if (size_of_bucket > DOP) {
        expandT.resume();
//        printf(">>>>>>>>>> %d %d %d %lf\n", round, parallelism, size_of_bucket, parallelism * 1.0 / size_of_bucket);
        num_of_block = DOP / v_per_block;
        int left_frontier = size_of_bucket - DOP;
        if (DOP % v_per_block) {
          num_of_block++;
        }
        SS = N * SearchSpace / num_of_block;
        ck_expand_bucket_with_filter<Label_t, Edge_t> <<< num_of_block, block >>>(SS,
                                                                                  &bucket[left_frontier],
                                                                                  buffer,
                                                                                  status,
                                                                                  g->gg,
                                                                                  mutex,
                                                                                  dominance_list_d,
                                                                                  dominance_list_len,
                                                                                  C_l[0],
                                                                                  DOP,
                                                                                  dynamic_list_d,
                                                                                  dynamic_list_len);
        cudaDeviceSynchronize();
        cuda_err_chk(cudaGetLastError());
        expandT.pause();

        pruneT.resume();
        ck_recheck_bucket<Label_t, Edge_t><<<num_of_block, 128>>>(SS,                                                                           bucket,
                                                                           buffer,
                                                                           status,
                                                                           dominance_list_d,
                                                                           dominance_list_len,
                                                                           left_frontier,
                                                                           dynamic_list_d,
                                                                           dynamic_list_len);
        cudaDeviceSynchronize();
        frontier_sum += DOP;
        pruneT.pause();

      } else {
        expandT.resume();
        num_of_block = size_of_bucket / v_per_block;
        if (size_of_bucket % v_per_block) {
          num_of_block++;
        }
//        REACH_HERE;
//        printf("NUM blk %d %d %d\n", num_of_block, size_of_bucket, v_per_block);
        SS = N * SearchSpace / num_of_block;
        ck_expand_bucket_with_filter<Label_t, Edge_t> <<< num_of_block, block >>>(SS,
                                                                                  bucket,
                                                                                  buffer,
                                                                                  status,
                                                                                  g->gg,
                                                                                  mutex,
                                                                                  dominance_list_d,
                                                                                  dominance_list_len,
                                                                                  C_l[0],
                                                                                  size_of_bucket,
                                                                                  dynamic_list_d,
                                                                                  dynamic_list_len);
        cudaDeviceSynchronize();
        frontier_sum += size_of_bucket;
        expandT.pause();
      }
      cuda_err_chk(cudaGetLastError());
      otherT.resume();
      ////////// bucket Time
      thrust::inclusive_scan(status_ptr, status_ptr + num_of_block, scan_out_ptr);
      cudaDeviceSynchronize();
      cudaMemcpy(&new_size_of_bucket, &scan_out[num_of_block - 1], sizeof(int), cudaMemcpyDeviceToHost);
      cudaDeviceSynchronize();
      max_bucket_size = max(max_bucket_size, new_size_of_bucket);
      ////////// bucket Time
      otherT.pause();
//      printf("%d %d\n", round, new_size_of_bucket);

      // change K size
      double log_bucket = log2(new_size_of_bucket);
      if (log_bucket < frontier_size.back()
          && frontier_size.back() < frontier_size[frontier_size.size() - 2]
          && frontier_size[frontier_size.size() - 2] < frontier_size[frontier_size.size() - 3]) {
        // ending part
        DOP = min(MAX_LABELS_IN_FRONTIER, DOP + DOP / 2);
        frontier_size.push_back(log_bucket);
      } else {
        frontier_size.push_back(log_bucket);
        if (frontier_size.size() > AR_DEGREE) {
          if (!AR::AutoRegression(&frontier_size[0], frontier_size.size(), AR_DEGREE, coefficients)) {
            printf("AR error!\n");
          } else {
            est = 0;
            for (int j = 0; j < AR_DEGREE; j++) {
              est += coefficients[j] * frontier_size[frontier_size.size() - j - 1];
            }
            est = exp2(est);
//            printf("EST is %lf %d\n", est, new_size_of_bucket);
            if (est > new_size_of_bucket) {
              double p = calculate_p<Label_t>(est, new_size_of_bucket, MAX_LABELS_IN_FRONTIER);
              DOP = (int) std::round(p);
//              parallelism = max((int) (std::round(p / NUM_SM / v_per_block)), 1) * NUM_SM * v_per_block;
//              printf("EST is %lf %d\n", est, new_size_of_bucket);
//              printf("<><><>< %d %lf\n", parallelism, p);
              f = true;
            }
          }
        }
      }

      if (new_size_of_bucket == 0) {
        break;
      }
      otherT.resume();
      ////////// bucket Time
      ck_fill_buckets<Label_t> << < num_of_block, BLOCK_SIZE >> > (SS, buffer, status, scan_out, bucket);
      cudaDeviceSynchronize();
      size_of_bucket = new_size_of_bucket;
      ////////// bucket Time
      otherT.pause();
    }
    t.toc();
    t.print_ms("Execution Time");
    cudaDeviceSynchronize();

    // Get result back
//    int *dm_size = (int *) malloc(sizeof(int) * N);
//    cudaMemcpy(dm_size, dynamic_list_len, sizeof(int) * N / _K, cudaMemcpyDeviceToHost);

//    int pos = SearchSpace * sink_l[0];
//    auto t1 = dominance_list[pos].d;
//    auto t2 = dominance_list[pos].delay;
//    for (int i = 0; i < dominance_list_len[sink_l[0]]; i++) {
//      if (dominance_list[pos+i].d < t1) {
//        t1 = dominance_list[pos+i].d;
//        t2 = dominance_list[pos+i].delay;
//      }
//    }
//    std::cout << std::fixed << "<><><><><><><>" << t2 << " " << t1 << std::endl;


//    t.print_ns("Execution Time");
//    std::cout << "Bucket management: " << double (otherT.total_time_ns()) / 1000000 << " ms" << std::endl;
//    std::cout << "Expand Time: " << double (t_ex.total_time_ns()) / 1000000 << " ms" << std::endl;
//    std::cout << "Prune Time: " << double (pruneT.total_time_ns() + (expandT.total_time_ns() - t_ex.total_time_ns())) / 1000000 << " ms" << std::endl;
//    std::cout << "K == " << parallelism << std::endl;
    cudaDeviceSynchronize();
  }
  return 0;
}

template<typename Label_t>
__global__ void ck_CSPP_path(int source, Label_t **dominance_list, int *list_len, int *path, int *mutex) {
  __shared__ int len;
  len = *list_len;
  __shared__
  Label_t final_label;
  final_label.d = LRG_FLOAT;

  // find the best result in dominance_list
  __syncthreads();
  if (len == 0) {
//    printf("%d\n", source);
    return;
  }
  int x = threadIdx.x;
  __syncthreads();
  while (x < len) {
//    printf("%d -- %.2lf %d\n", x, (*dominance_list)[x].d, (*dominance_list)[x].delay);
    if (final_label.d > (*dominance_list)[x].d) {
      bool leave = true;
      while (leave) {
        if (atomicCAS(mutex, 0, 1) == 0) {
          if (final_label.d > (*dominance_list)[x].d) {
            Label_t::copy(&final_label, &(*dominance_list)[x]);
            __threadfence();
          }
          atomicExch(mutex, 0);
          leave = false;
        }
      }
    }
    x += blockDim.x;
  }
  __syncthreads();
  // generate path
  if (threadIdx.x == 0) {
    Label_t *ptr_label = &final_label;
    int h = 0;
    while (ptr_label != NULL) {
//      printf("%d ", ptr_label->father_id);
      path[h++] = ptr_label->father_id;
      ptr_label = ptr_label->pre;
    }
//    printf("%d\n", source);
  }
  __syncthreads();
}

#endif


<|start_filename|>algorithm/cspp_cpu.hpp<|end_filename|>
#ifndef _cspp_cpu_hpp_
#define _cspp_cpu_hpp_
#include <vector>
#include <map>
#include "../include/label.hpp"
#include "../include/constraint.hpp"
#include "../include/edge.hpp"
#include "../include/graph.hpp"
#include "../include/util.hpp"
#include "../include/TimeMeasurer.hpp"
#include "../include/graph.hpp"
#include <omp.h>

#define LRG_FLOAT 1000000

template<typename Label_t>
void add_new_frontier(std::map<int, std::vector<Label_t>> &next_round, Label_t &new_label) {
  if (next_round.find(new_label.vid) != next_round.end()) {
    bool dominate = true;
    for (auto it = next_round[new_label.vid].begin(); it != next_round[new_label.vid].end();) {
      if (Label_t::dominate_cpu(*it, new_label)) {
        // new label is dominated
        dominate = false;
        break;
      }
      if (Label_t::dominate_cpu(new_label, *it)) {
        // new label dominates some other labels
        it = next_round[new_label.vid].erase(it);
      } else {
        it++;
      }
    }
    if (!dominate) {
      return; // go to next edge
    }
    next_round[new_label.vid].push_back(new_label);
  } else {
    std::vector<Label_t> new_frontier;
    new_frontier.push_back(new_label);
    next_round[new_label.vid] = new_frontier;
  }
}

template<typename Label_t, typename Edge_t>
int RUN_CSPP_CPU(
    Graph<Edge_t> *g,
    int ntask,
    vector<Label_t> &init_state_l,
    vector<int> &sink_l,
    vector<int *> &path_l,
    vector<int> &size_of_result_l,
    vector<c_constraint<Label_t>> &C_l) {

//====================
// timer starts
//====================

  TimeMeasurer total_time;
  total_time.tic();
  for (int i = 0; i < ntask; i++) {
//    std::cout << C_l[i].delay_L << std::endl;
    std::cout << "==================" << std::endl;
    // prepare data
    auto src = init_state_l[i];
    std::vector<Label_t> frontiers;
    std::map<int, std::vector<Label_t>> dominace_list;
    std::map<int, std::vector<Label_t>> next_round;
    frontiers.push_back(src);
    TimeMeasurer t;
    t.tic();
    {
      while (!frontiers.empty()) {
//      std::cout << "==================" << std::endl;
        for (int ff = 0; ff < frontiers.size(); ff++) {
          auto frontier = frontiers[ff];
          // expand each of it
          for (auto e : g->v_l[frontier.vid].edges) {
            Label_t new_label;
            frontier.expand_cpu(e->snk, *e, &new_label);
            if (!C_l[i].pass_constraints_check_cpu(new_label)) {
              continue;
            }
            if (dominace_list.find(new_label.vid) != dominace_list.end()) {
              bool dominate = true;
              for (auto it = dominace_list[new_label.vid].begin(); it != dominace_list[new_label.vid].end();) {
                if (Label_t::dominate_cpu(*it, new_label)) {
                  // new label is dominated
                  dominate = false;
                  break;
                }
                if (Label_t::dominate_cpu(new_label, *it)) {
                  // new label dominates some other labels
                  it = dominace_list[new_label.vid].erase(it);
                } else {
                  it++;
                }
              }
              if (!dominate) {
                continue; // go to next edge
              }
              dominace_list[new_label.vid].push_back(new_label);
              add_new_frontier(next_round, new_label);
            } else {
              std::vector<Label_t> new_dominace_vec;
              new_dominace_vec.push_back(new_label);
              dominace_list[new_label.vid] = new_dominace_vec;
              add_new_frontier(next_round, new_label);
            }
          }
        }
        frontiers.clear();
        for (auto i : next_round) {
          frontiers.insert(frontiers.end(), i.second.begin(), i.second.end());
        }
        next_round.clear();
      }
    }
    t.toc();
    t.print_ms("==");
  }
  total_time.toc();
  total_time.print_ms("Total Time is");

  return 0;
}

template<typename Label_t, typename Edge_t>
int RUN_CSPP_OMP(
    Graph<Edge_t> *g,
    int ntask,
    vector<Label_t> &init_state_l,
    vector<int> &sink_l,
    vector<int *> &path_l,
    vector<int> &size_of_result_l,
    vector<c_constraint<Label_t>> &C_l) {

//====================
// timer starts
//====================

  TimeMeasurer total_time;
  total_time.tic();
  for (int i = 0; i < ntask; i++) {
//    std::cout << C_l[i].delay_L << std::endl;
    std::cout << "==================" << std::endl;
    // prepare data
    auto src = init_state_l[i];
    std::vector<Label_t> frontiers;
    std::vector<Label_t> next_round;

    std::vector<std::vector<Label_t>> dominace_list;
    for (int j = 0; j < g->N; j++) {
      dominace_list.emplace_back(std::vector<Label_t>());
    }
    frontiers.push_back(src);
    TimeMeasurer t;
    t.tic();

    while (!frontiers.empty()) {
#pragma omp parallel for
      for (int ff = 0; ff < frontiers.size(); ff++) {
        auto frontier = frontiers[ff];
        // expand each of it
        for (auto e : g->v_l[frontier.vid].edges) {
          Label_t new_label;
          frontier.expand_cpu(e->snk, *e, &new_label);
          if (!C_l[i].pass_constraints_check_cpu(new_label)) {
            continue;
          }
          bool dominate = true;
          for (auto it = dominace_list[new_label.vid].begin(); it != dominace_list[new_label.vid].end();) {
            if (Label_t::dominate_cpu(*it, new_label)) {
              // new label is dominated
              dominate = false;
              break;
            }
            if (Label_t::dominate_cpu(new_label, *it)) {
              // new label dominates some other labels
#pragma omp critical
              it = dominace_list[new_label.vid].erase(it);
            } else {
              it++;
            }
          }
          if (!dominate) {
            continue; // go to next edge
          }
#pragma omp critical
          {
            dominace_list[new_label.vid].push_back(new_label);
            next_round.push_back(new_label);
          }
        }
      }
#pragma omp single
      {
        frontiers.swap(next_round);
        next_round.clear();
      }
#pragma omp barrier
    }

    for (auto i : dominace_list[g->N - 1]) {
      std::cout << i.d << " " << i.delay << std::endl;
    }
    t.toc();
    t.print_ms("==");
  }
  total_time.toc();
  total_time.print_ms("Total Time is");

  return 0;
}
#endif


<|start_filename|>include/comm.hpp<|end_filename|>
#ifndef _comm_hpp_
#define _comm_hpp_

#include "util.hpp"
#include "graph.hpp"

#include "io.hpp"
#include "../algorithm/cspp.hpp"
#endif


<|start_filename|>include/task.hpp<|end_filename|>
#ifndef _task_hpp_
#define _task_hpp_

#include <vector>

using namespace std;

class task;

class task {
 public:
  int ntask;
  vector<int> source_l_;
  vector<int> sink_l_;
  vector<int> delay_l_;
  vector<int> hop_l_;
  vector<int> cost_l_;

 public:
  task() {
    ;
  }
  task(int ntask) :
      ntask(ntask) {
    ;
  }
  int read_query(char *file_name, int graph_type) {
    char temp[255];
    FILE *topo_f;
    int id, src, snk, delay, hop, cost;
    topo_f = fopen(file_name, "r");
    vector<tuple<int, int, int, int>> task_l;
    if (graph_type == 1) {
      fgets(temp, 255, topo_f);
      while (fscanf(topo_f, "%d", &id) != EOF) {
        fscanf(topo_f, ",%d", &src);
        fscanf(topo_f, ",%d", &snk);
        fscanf(topo_f, ",%d", &delay);
        fscanf(topo_f, ",%d", &hop);
        fscanf(topo_f, ",%d", &cost);
        source_l_.push_back(src);
        sink_l_.push_back(snk);
        delay_l_.push_back(delay);
        hop_l_.push_back(hop);
        cost_l_.push_back(cost);
      }
    } else {
      int u, v, cost_lower, cost_upper, cost_limit;
      for (; fscanf(topo_f, "%d %d %d %d %d\n", &u, &v, &cost_lower, &cost_upper, &cost_limit) == 5;) {
        source_l_.push_back(u);
        sink_l_.push_back(v);
        delay_l_.push_back(cost_limit);
      }
    }
    fclose(topo_f);
    ntask = source_l_.size();

    return 0;
  }

  // used for COLA data set
  int load_query(char *file_name) {
    FILE *topo_f;
    topo_f = fopen(file_name, "r");
    vector<tuple<int, int, int, int>> task_l;
    int u, v, cost_lower, cost_upper, cost_limit;
    for (; fscanf(topo_f, "%d %d %d %d %d\n", &u, &v, &cost_lower, &cost_upper, &cost_limit) == 5;) {
      source_l_.push_back(u);
      sink_l_.push_back(v);
      delay_l_.push_back(cost_limit);
    }
    fclose(topo_f);
    ntask = source_l_.size();

    return 0;
  }
};

#endif


<|start_filename|>algorithm/bucket.cuh<|end_filename|>
#ifndef _bucket_cuh_
#define _bucket_cuh_

#include <cuda_runtime.h>
#include "dominance.cuh"
#include "../include/gpu_graph.hpp"

#define LRG_FLOAT 1000000
#define BLOCK_SIZE 32
#define TABLE_SIZE 96
#define VSIZE 2

#define HASH

template<typename Label_t, typename Edge_t>
__global__ void ck_expand_bucket_with_filter(
    int local_size,
    Label_t *bucket,
    Label_t *buffer,
    int *st,
    gpuGraph<Edge_t> *gg,
    int *mutex,
    Label_t **dominance_list,
    int *dominance_list_len,
    c_constraint<Label_t> C,
    int size_of_bucket,
    Label_t **dynamic_list,
    int *dynamic_list_len) {
  __shared__ Label_t hash_table[TABLE_SIZE];
  __shared__ int maxR;
#ifdef HASH
  __shared__ int mtx;
  __shared__ int hs_cout;
  mtx = 0;
  hs_cout = 0;
#endif
  int tid = threadIdx.y * blockDim.x + threadIdx.x;
  while (tid < TABLE_SIZE) {
    hash_table[tid].vid = -1;
    tid += BLOCK_SIZE;
  }
  st[blockIdx.x] = 0;
  maxR = 0;
  int l, r;
  Label_t *local_buffer = buffer + blockIdx.x * local_size;
  Label_t *cur_v;
  tid = (blockIdx.x * blockDim.y + threadIdx.y) * VSIZE;
  for (int T = 0; T < VSIZE; T++, tid++) {
    if (tid < size_of_bucket) {
      cur_v = bucket + tid;
      l = gg->_count_row_d[cur_v->vid];
      r = gg->_count_row_d[cur_v->vid + 1];
    } else {
      cur_v = nullptr;
      l = r = -1; // give up processing the last empty frontier
    }
    // get the max rounds
    __syncthreads();
    if (threadIdx.x == 0) {
      for (int i = 0; i < BLOCK_SIZE / 32; i++) {
        if (r - l > maxR) {
          maxR = r - l;
        }
      }
    }
    __syncthreads();
//    if (threadIdx.x % 32 == 0) {
//      printf("id %d %d %d %d %d %d\n", blockIdx.x, threadIdx.y * 32 + threadIdx.x, blockIdx.x * blockDim.y + threadIdx.y, tid, r, maxR);
//    }
    l = threadIdx.x + l;
    for (; maxR > 0; maxR -= blockDim.x) {
      if (l < r) {
        Edge_t *next_e = gg->el_d + gg->_eid_col_d[l];
        Label_t new_label;
        cur_v->expand(gg->_vid_col_d[l], *next_e, &new_label);
        if (!C.pass_constraints_check(new_label)) {
          break;
        }
#ifdef HASH
        int pos = new_label.vid % TABLE_SIZE;

//          printf("id %d %d %d ***\n", tid, new_label.vid, pos);
        bool inserted = false;
        while (!inserted) {
          if (hash_table[pos].vid == -1) {
            bool leave = true;
            while (leave) {
              if (atomicCAS(&mtx, 0, 1) == 0) {
                if (hash_table[pos].vid == -1) {
                  Label_t::copy(&hash_table[pos], &new_label);
                  hs_cout++;
                  inserted = true;
                }
                __threadfence();
                atomicExch(&mtx, 0);
                leave = false;
              }
            }
          }
          if (inserted || Label_t::dominate_vid(hash_table[pos], new_label)) {
            break;
          }
          if (Label_t::dominate_vid(new_label, hash_table[pos])) {
            bool leave = true;
            while (leave) {
              if (atomicCAS(&mtx, 0, 1) == 0) {
                if (Label_t::dominate_vid(new_label, hash_table[pos])) {
                  Label_t::copy(&hash_table[pos], &new_label);
                  hs_cout++;
                  inserted = true;
                }
                __threadfence();
                atomicExch(&mtx, 0);
                leave = false;
              }
            }
          }
          pos = (pos + 1) % TABLE_SIZE;
        }
#else
        hash_table[threadIdx.y * blockDim.x + threadIdx.x] = new_label;
        for (int j = 0; j < blockDim.x * blockDim.y; j++) {
            if (hash_table[j].vid > 0) {
                int vid = hash_table[j].vid;
                second_level_pruning<Label_t>(&hash_table[j],
                                              dominance_list[vid],
                                              dominance_list_len[vid],
                                              dynamic_list[vid / _K],
                                              dynamic_list_len[vid / _K],
                                              local_buffer,
                                              &st[blockIdx.x],
                                              mutex);
            }
        }
#endif
      }
      __syncthreads();
#ifdef HASH
      if (hs_cout + BLOCK_SIZE > TABLE_SIZE) {
        hs_cout = 0;
#ifdef DUPLICATE_EDGE
        // works only for duplicated edge
       for (int j = 0; j < BLOCK_SIZE; j++) {
         if (new_label[j].vid < 0) {
           break;
         }
         if (threadIdx.x != j) {
           if (new_label[threadIdx.x].vid == new_label[j].vid && Label_t::dominate(new_label[threadIdx.x], new_label[j])) {
             new_label[j].delay = LRG_FLOAT;
           }
         }
         __syncthreads();
       }
#endif
        for (int j = 0; j < TABLE_SIZE; j++) {
          int vid = hash_table[j].vid;
          if (vid >= 0) {
            second_level_pruning<Label_t>(&hash_table[j],
                                          dominance_list[vid],
                                          dominance_list_len[vid],
                                          dynamic_list[vid / _K],
                                          dynamic_list_len[vid / _K],
                                          local_buffer,
                                          &st[blockIdx.x],
                                          mutex);

          }
        }
      }
#else
      ;;
#endif
      l += blockDim.x;
      __syncthreads();
    }
  }

#ifdef HASH

  for (int j = 0; j < TABLE_SIZE; j++) {
    int vid = hash_table[j].vid;
    if (vid >= 0) {
      second_level_pruning<Label_t>(&hash_table[j],
                                    dominance_list[vid],
                                    dominance_list_len[vid],
                                    dynamic_list[vid / _K],
                                    dynamic_list_len[vid / _K],
                                    local_buffer,
                                    &st[blockIdx.x],
                                    mutex);

    }
  }
#endif
}

template<typename Label_t>
__global__ void ck_fill_buckets(int local_size, Label_t *buffer, int *st, int *histogram, Label_t *bucket) {
  int base;
  int x, index;
  if (blockIdx.x == 0) {
    base = 0;
  } else {
    base = histogram[blockIdx.x - 1];
  }
  __syncthreads();
  index = blockIdx.x * local_size;
  x = threadIdx.x;
  __syncthreads();
  while (x < st[blockIdx.x]) {
    bucket[base + x] = buffer[index + x];
//    printf("%d %.3lf %.3lf\n", bucket[base + x].vid, bucket[base + x].d, bucket[base + x].delay);
    x += blockDim.x;
  }
}

template<typename Label_t, typename Edge_t>
__global__ void ck_recheck_bucket(
    int local_size,
    Label_t *bucket,
    Label_t *buffer,
    int *st,
    Label_t **dominance_list,
    int *dominance_list_len,
    int size_of_bucket,
    Label_t **dynamic_list,
    int *dynamic_list_len) {
  int tid = blockIdx.x;
  Label_t *local_buffer = buffer + blockIdx.x * local_size;
  while (tid < size_of_bucket) {
    int vid = bucket[tid].vid;
//        if (threadIdx.x == 0) {
//            printf("id %d %d %d\n", blockIdx.x, tid, vid);
//        }

    second_level_pruning_only<Label_t>(&bucket[tid],
                                       dominance_list[vid],
                                       dominance_list_len[vid],
                                       dynamic_list[vid / _K],
                                       dynamic_list_len[vid / _K],
                                       local_buffer,
                                       &st[blockIdx.x]);
    tid += gridDim.x;
  }
}

template<typename Label_t, typename Edge_t>
__global__ void ck_expand_bucket_with_nothing(
    int local_size,
    Label_t *bucket,
    Label_t *buffer,
    int *st,
    gpuGraph<Edge_t> *gg,
    int *mutex,
    Label_t **dominance_list,
    int *dominance_list_len,
    c_constraint<Label_t> C,
    int size_of_bucket,
    Label_t **dynamic_list,
    int *dynamic_list_len) {
  __shared__ Label_t hash_table[TABLE_SIZE];
  __shared__ int maxR;
#ifdef HASH
  __shared__ int mtx;
  __shared__ int hs_cout;
  mtx = 0;
  hs_cout = 0;
#endif
  int tid = threadIdx.y * blockDim.x + threadIdx.x;
  while (tid < TABLE_SIZE) {
    hash_table[tid].vid = -1;
    tid += BLOCK_SIZE;
  }
//  st[blockIdx.x] = 0;
  maxR = 0;
  int l, r;
//  Label_t *local_buffer = buffer + blockIdx.x * local_size;
  Label_t *cur_v;
  tid = (blockIdx.x * blockDim.y + threadIdx.y) * VSIZE;
  for (int T = 0; T < VSIZE; T++, tid++) {
    if (tid < size_of_bucket) {
      cur_v = bucket + tid;
      l = gg->_count_row_d[cur_v->vid];
      r = gg->_count_row_d[cur_v->vid + 1];
    } else {
      cur_v = nullptr;
      l = r = -1; // give up processing the last empty frontier
    }
    // get the max rounds
    __syncthreads();
    if (threadIdx.x == 0) {
      for (int i = 0; i < BLOCK_SIZE / 32; i++) {
        if (r - l > maxR) {
          maxR = r - l;
        }
      }
    }
    __syncthreads();
//        if (threadIdx.x % 32 == 0) {
//            printf("id %d %d %d %d %d %d\n", blockIdx.x, threadIdx.y * 32 + threadIdx.x,
//                   blockIdx.x * blockDim.y + threadIdx.y, tid, r, maxR);
//        }
    l = threadIdx.x + l;
    for (; maxR > 0; maxR -= blockDim.x) {
      if (l < r) {
        Edge_t *next_e = gg->el_d + gg->_eid_col_d[l];
        Label_t new_label;
        cur_v->expand(gg->_vid_col_d[l], *next_e, &new_label);
//        if (!C.pass_constraints_check(new_label)) {
//          break;
//        }
#ifdef HASH
        int pos = new_label.vid % TABLE_SIZE;

//          printf("id %d %d %d ***\n", tid, new_label.vid, pos);
        bool inserted = false;
        while (!inserted) {
          if (hash_table[pos].vid == -1) {
            bool leave = true;
            while (leave) {
              if (atomicCAS(&mtx, 0, 1) == 0) {
                if (hash_table[pos].vid == -1) {
//                  Label_t::copy(&hash_table[pos], &new_label);
                  hs_cout++;
                  inserted = true;
                }
                __threadfence();
                atomicExch(&mtx, 0);
                leave = false;
              }
            }
          }
          if (inserted || Label_t::dominate_vid(hash_table[pos], new_label)) {
            break;
          }
          if (Label_t::dominate_vid(new_label, hash_table[pos])) {
            bool leave = true;
            while (leave) {
              if (atomicCAS(&mtx, 0, 1) == 0) {
                if (Label_t::dominate_vid(new_label, hash_table[pos])) {
//                  Label_t::copy(&hash_table[pos], &new_label);
                  hs_cout++;
                  inserted = true;
                }
                __threadfence();
                atomicExch(&mtx, 0);
                leave = false;
              }
            }
          }
          pos = (pos + 1) % TABLE_SIZE;
        }
#else
        hash_table[threadIdx.y * blockDim.x + threadIdx.x] = new_label;
        for (int j = 0; j < blockDim.x * blockDim.y; j++) {
            if (hash_table[j].vid > 0) {
                int vid = hash_table[j].vid;
                //            second_level_pruning<Label_t>(&hash_table[j],
                //                                          dominance_list[vid],
                //                                          dominance_list_len[vid],
                //                                          dynamic_list[vid / _K],
                //                                          dynamic_list_len[vid / _K],
                //                                          local_buffer,
                //                                          &st[blockIdx.x],
                //                                          mutex);
            }
        }
#endif
      }
      __syncthreads();
#ifdef HASH
      if (hs_cout + BLOCK_SIZE > TABLE_SIZE) {
        hs_cout = 0;
#ifdef DUPLICATE_EDGE
        // works only for duplicated edge
       for (int j = 0; j < BLOCK_SIZE; j++) {
         if (new_label[j].vid < 0) {
           break;
         }
         if (threadIdx.x != j) {
           if (new_label[threadIdx.x].vid == new_label[j].vid && Label_t::dominate(new_label[threadIdx.x], new_label[j])) {
             new_label[j].delay = LRG_FLOAT;
           }
         }
         __syncthreads();
       }
#endif
        for (int j = 0; j < TABLE_SIZE; j++) {
          int vid = hash_table[j].vid;
          if (vid >= 0) {
//            second_level_pruning<Label_t>(&hash_table[j],
//                                          dominance_list[vid],
//                                          dominance_list_len[vid],
//                                          dynamic_list[vid / _K],
//                                          dynamic_list_len[vid / _K],
//                                          local_buffer,
//                                          &st[blockIdx.x],
//                                          mutex);

          }
        }
      }
#else
      ;;
#endif
      l += blockDim.x;
      __syncthreads();
    }
  }

#ifdef HASH

  for (int j = 0; j < TABLE_SIZE; j++) {
    int vid = hash_table[j].vid;
    if (vid >= 0) {
//      second_level_pruning<Label_t>(&hash_table[j],
//                                    dominance_list[vid],
//                                    dominance_list_len[vid],
//                                    dynamic_list[vid / _K],
//                                    dynamic_list_len[vid / _K],
//                                    local_buffer,
//                                    &st[blockIdx.x],
//                                    mutex);

    }
  }
#endif
}

#endif


<|start_filename|>include/AR.hpp<|end_filename|>
#ifndef CSP_AR_HPP
#define CSP_AR_HPP
#define MAXENTROPY      0
#define LEASTSQUARES    1
#define TRUE  1
#define FALSE 0
#define MAXCOEFF 100

class AR {
 public:
  static int AutoRegression(double *, int, int, double *, int method = 0);
  static int ARMaxEntropy(double *, int, int, double **, double *, double *, double *, double *);
  static int ARLeastSquare(double *, int, int, double *);
  static int SolveLE(double **, double *, unsigned int);
};

#endif //CSP_AR_HPP


<|start_filename|>include/io.hpp<|end_filename|>
#ifndef _io_hpp_
#define _io_hpp_
#include <stdio.h>
#include <iostream>
#include <vector>
#include <string>
#include "util.hpp"
#include "graph.hpp"
#include "task.hpp"
using namespace std;

extern float *A, *B, *Theta;

/**
 * @brief initialize graph
 * @tparam Graph_t
 * @param file_name
 * @param g
 * @param graph_type 1 for our graph, 2 for xxk graph
 * @return
 */
template<class Graph_t>
int read_graph(char *file_name, Graph_t *g, int graph_type) {
  char temp[255];
  FILE *topo_f;
  int id, src, snk;
  int weight, delay;
  topo_f = fopen(file_name, "r");
  int num_of_vertex, num_of_edge;
  if (graph_type == 1) {
    fscanf(topo_f, "%d", &num_of_vertex);
    fscanf(topo_f, "%c", temp);
    fscanf(topo_f, "%d", &num_of_edge);
    fscanf(topo_f, "%c", temp);
    num_of_edge *= 2;
    g->v_l.clear();
    g->v_l.reserve(num_of_vertex);
    g->N = num_of_vertex;
    g->e_l.clear();
    g->e_l.resize(num_of_edge);
    g->M = num_of_edge;
    using myVertex = Vertex<Graph_t>;
    for (int i = 0; i < num_of_vertex; i++) {
      myVertex *v = new myVertex(i);
      g->v_l.push_back(*v);
    }
    for (int i = 0; i < num_of_edge; i++) {
      g->e_l[i] = new Edge(i);
    }
    while (fscanf(topo_f, "%d", &src) != EOF) {
      fscanf(topo_f, "%d %d %d\n", &snk, &weight, &delay);
//      src--;
//      snk--;
      (g->e_l[id])->assign(id, src, snk, weight, delay);
      g->v_l[src].assign(g->e_l[id], g);
      id += 1;
      (g->e_l[id])->assign(id, snk, src, weight, delay);
      g->v_l[snk].assign(g->e_l[id], g);
      id += 1;
    }
  } else {
    fscanf(topo_f, "%d", &num_of_vertex);
    g->v_l.clear();
    g->v_l.reserve(num_of_vertex);
    g->N = num_of_vertex;
    g->e_l.clear();
    using myVertex = Vertex<Graph_t>;
    for (int i = 0; i < num_of_vertex; i++) {
      myVertex *v = new myVertex(i);
      g->v_l.push_back(*v);
    }
    int w, c, u, v, M = 0;
    for (; fscanf(topo_f, "%d:", &u) == 1;) {
      for (; fscanf(topo_f, "%d", &v) == 1;) {
        if (v == -1) break;
        int ret_val = fscanf(topo_f, "%d %d", &w, &c);
        if (ret_val < 0) {
          cout << "input file error" << endl;
          exit(0);
        }
        g->e_l.push_back(new Edge(M));
        g->e_l.back()->assign(M, u, v, w, c);
        g->v_l[u].assign(g->e_l[M++], g);
//        g->e_l.push_back(new Edge_Index(M));
//        g->e_l.back()->assign(M, v, u, w, c);
//        g->v_l[v].assign(g->e_l[M++], g);
      }
    }
    g->M = g->e_l.size();
    fclose(topo_f);
  }
  return 0;
}

#endif


<|start_filename|>include/util.hpp<|end_filename|>
#ifndef _util_hpp_
#define _util_hpp_

#include <iostream>

#define __DEBUG

#ifdef __DEBUG
#define DEBUG_TEST 1
#else
#define DEBUG_TEST 0
#endif

inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true) {
  if (code != cudaSuccess) {
    fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
  }
}

#define cuda_err_chk(ans) { gpuAssert((ans), __FILE__, __LINE__); }

#endif



<|start_filename|>include/constraint.hpp<|end_filename|>
#ifndef _csp_hpp_
#define _csp_hpp_
#include <cuda_runtime.h>
#include "label.hpp"

template<typename Label_t>
class constraint {
 public:
  constraint() {
    return;
  }
  virtual void set_constraints() {
    return;
  }

  __device__ virtual bool pass_constraints_check(Label_t l) {
    return false;
  }
};

#if CSP == 2
template<typename Label_t>
class c_constraint {
public	:
    float delay_L;
    int hop_L;
public	:
    c_constraint() {
        return;
    }
    c_constraint(float d, int h = 100) :
    delay_L(d), hop_L(h) {
        return;
    }
    void set_constraints(float d, int h) {
        delay_L = d;
        hop_L = h;

        return;
    }

    __device__ bool pass_constraints_check(const Label_t &l) {
        if (l.delay > delay_L || l.hop > hop_L) {
            return false;
        }

        return true;
    }
};
#elif CSP == 3
template<typename Label_t>
class c_constraint {
public	:
    float delay_L;
    int hop_L;
    float maxcost_L;
public	:
    c_constraint() {
        return;
    }
    c_constraint(float d, int h = 100, float m = 10000) :
    delay_L(d), hop_L(h), maxcost_L(m) {
        return;
    }
    void set_constraints(float d, int h, int m) {
        delay_L = d;
        hop_L = h;
        maxcost_L = m;

        return;
    }

    __device__ bool pass_constraints_check(const Label_t &l) {
        if (l.delay > delay_L || l.hop > hop_L || l.maxcost > maxcost_L) {
            return false;
        }

        return true;
    }
};
#elif CSP == 1
template<typename Label_t>
class c_constraint {
 public    :
  float delay_L;
 public    :
  c_constraint() {
    return;
  }
  c_constraint(float d) : delay_L(d) {
    return;
  }
  void set_constraints(int d) {
    delay_L = d;

    return;
  }

  __device__ bool pass_constraints_check(const Label_t &l) {
    return !(l.delay > delay_L);
  }
  bool pass_constraints_check_cpu(const Label_t &l) {
    return !(l.delay > delay_L);
  }
};
#elif CSP == 4 // 2-OSP
template<typename Label_t>
class c_constraint {
public	:
    float delay_L;
public	:
    c_constraint() {
        return;
    }
    c_constraint(float d) :
    delay_L(d) {
        return;
    }
    void set_constraints(int d) {
        delay_L = d;

        return;
    }

    __device__ bool pass_constraints_check(const Label_t &l) {
        if (l.delay > delay_L) {
            return false;
        }

        return true;
    }
};

#elif CSP == 5 //Social
template<typename Label_t>
class c_constraint {
public	:
    float delay_L;
public	:
    c_constraint() {
        return;
    }
    c_constraint(float d) :
    delay_L(d) {
        return;
    }
    void set_constraints(int d) {
        delay_L = d;

        return;
    }

    __device__ bool pass_constraints_check(const Label_t &l) {
        if (l.delay > delay_L) {
            return false;
        }

        return true;
    }
};

#else // network qos
template<typename Label_t>
class c_constraint {
 public    :
  float delay_L;
 public    :
  c_constraint() {
    return;
  }
  c_constraint(float d) : delay_L(d) {
    return;
  }
  void set_constraints(int d) {
    delay_L = d;

    return;
  }

  __device__ bool pass_constraints_check(const Label_t &l) {
    return !(l.delay < delay_L);
  }
  bool pass_constraints_check_cpu(const Label_t &l) {
    return !(l.delay < delay_L);
 }
};
#endif

#endif


<|start_filename|>include/TimeMeasurer.hpp<|end_filename|>
#ifndef _TIMEMEASURER_HPP
#define _TIMEMEASURER_HPP
#include <iostream>
#include <chrono>
using std::chrono::high_resolution_clock;
using std::chrono::milliseconds;
using std::chrono::microseconds;
using std::chrono::nanoseconds;

class TimeMeasurer {
 public:
  TimeMeasurer() {
    ns_count_ = 0;
  }
  ~TimeMeasurer() {}

  void tic() {
    start_time_ = high_resolution_clock::now();
  }

  void toc() {
    end_time_ = high_resolution_clock::now();
  }

  void resume() {
    tic();
  }
  void pause() {
    toc();
    ns_count_ += time_ns();
  }

  long long total_time_ns() const {
    return ns_count_;
  }

  long long time_ms() const {
    return std::chrono::duration_cast<milliseconds>(end_time_ - start_time_).count();
  }

  long long time_us() const {
    return std::chrono::duration_cast<microseconds>(end_time_ - start_time_).count();
  }

  long long time_ns() const {
    return std::chrono::duration_cast<nanoseconds>(end_time_ - start_time_).count();
  }

  void print_ms(std::string msg) const {
    std::cout << msg << ": ";
    print_ms();
  }
  void print_ms() const {
    std::cout << std::chrono::duration_cast<milliseconds>(end_time_ - start_time_).count() << " ms" << std::endl;
  }

  void print_us(std::string msg) const {
    std::cout << msg << ": ";
    print_us();
  }
  void print_us() const {
    std::cout << std::chrono::duration_cast<microseconds>(end_time_ - start_time_).count() << " us" << std::endl;
  }

  void print_ns(std::string msg) const {
    std::cout << msg << ": ";
    print_ns();
  }
  void print_ns() const {
    std::cout << std::chrono::duration_cast<nanoseconds>(end_time_ - start_time_).count() << " ns" << std::endl;
  }

 private:
  TimeMeasurer(const TimeMeasurer &);
  TimeMeasurer &operator=(const TimeMeasurer &);

 private:
  high_resolution_clock::time_point start_time_;
  high_resolution_clock::time_point end_time_;
  long long ns_count_;
};
#endif //_TIMEMEASURER_HPP


<|start_filename|>algorithm/cspp.hpp<|end_filename|>
#ifndef _cspp_hpp_
#define _cspp_hpp_
#include <cuda_runtime.h>
#include <vector>
#include "../include/label.hpp"
#include "../include/constraint.hpp"
#include "../include/graph.hpp"
#include "cspp.cuh"
#include "cspp_cpu.hpp"

template<typename Label_t, typename Edge_t>
class cspp {
  using Graph_t = Graph<Edge_t>;
  using Constraint_t = c_constraint<Label_t>;
 public:
  int *path;

  vector<Label_t> init_state;
  vector<int> sink_l;
  int ntask;
  vector<int *> path_l;
  vector<int> size_of_result_l;

  Graph_t *g;
  vector<Constraint_t> &C_l;
 public:
  // one constraints
  cspp(Graph_t *g, Label_t source_state, int sink, Constraint_t C) : g(g) {
    ntask = 1;
    init_state.push_back(source_state);
    C_l.push_back(C);
    path = (int *) malloc(g->N * sizeof(int));
    path_l.push_back(path);
  }

  cspp(Graph_t *g, int ntask, vector<Label_t> &source_state, vector<int> &sink_l, vector<Constraint_t> &C_l) :
      g(g), ntask(ntask), init_state(source_state), sink_l(sink_l), C_l(C_l) {
    // allocate double size for the path
    long N = g->N;
    int *p = (int *) malloc(sizeof(int) * N * (long) ntask);
    path_l.reserve(ntask);
    for (long i = 0; i < ntask; i++) {
      path_l.push_back(&p[N * i]);
    }
    size_of_result_l.reserve(ntask);
  }

  int solve() {
    std::cout << "number of tasks " << ntask << std::endl;
    RUN_CSPP<Label_t, Edge_t>(g, ntask, init_state, sink_l, path_l, size_of_result_l, C_l);

    return 0;
  }

  int solve_cpu() {
    RUN_CSPP_CPU<Label_t, Edge_t>(g, ntask, init_state, sink_l, path_l, size_of_result_l, C_l);

    return 0;
  }

  int solve_omp() {
    RUN_CSPP_OMP<Label_t, Edge_t>(g, ntask, init_state, sink_l, path_l, size_of_result_l, C_l);
  }

  int self_validate() {

    return 0;
  }
};

#endif


<|start_filename|>main.cu<|end_filename|>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <cuda_runtime.h>
#include "include/comm.hpp"
#include "include/task.hpp"
#include "include/label.hpp"
#include "include/constraint.hpp"
#include "include/label.hpp"

#ifndef CSP
#define CSP 1
#endif

using namespace std;

static Graph<c_edge> *g;

int pre_process_task(task *tasks, int N) { // to make sure source are within range
  for (int i = 0; i < tasks->ntask; i++) {
    tasks->source_l_[i] %= N;
    tasks->sink_l_[i] %= N;
  }

  return 0;
}

c_edge *create_edge_list(Graph<c_edge> *g) {
  c_edge *l = (c_edge *) malloc(sizeof(c_edge) * g->M);
  for (int i = 0; i < g->M; i++) {
    l[i].len = g->e_l[i]->len;
    l[i].delay = g->e_l[i]->delay;
  }

  return l;
}

int main(int argc, char *argv[]) {
  using Constraint = c_constraint<c_label>;
  using Graph = Graph<c_edge>;
  int graph_type, ntask;
  char graph_f[255], query_f[255];
  FILE *conf;
  conf = fopen("../inputs/config", "r");
  fscanf(conf, "%d", &graph_type);
  fscanf(conf, "%s", graph_f);
  fscanf(conf, "%s", query_f);
  fscanf(conf, "%d", &ntask);

  cudaSetDevice(0);
  g = new Graph();
  task *tasks = new task();
  read_graph(graph_f, g, graph_type);
  if (graph_type == 1) {
    tasks->read_query(query_f, 1);
    pre_process_task(tasks, g->N);
  } else {
    // only for COLA data set
    tasks->load_query(query_f);
  }
  g->adj_to_csr();
  g->el = create_edge_list(g);
  g->make_gpu_graph();

  cspp<c_label, c_edge> *solver;

  // create constraints, and source_l
  vector<Constraint> C_l;
  vector<c_label> source_state_l;
  for (int i = 0; i < ntask; i++) {
    c_label l(tasks->source_l_[i]);
    source_state_l.push_back(l);
    Constraint c(tasks->delay_l_[i]);
    C_l.push_back(c);
  }
  cuda_err_chk(cudaGetLastError());
  std::cout << "Number of V " << g->N << std::endl;
  std::cout << "Number of E " << g->M << std::endl;
  solver = new cspp<c_label, c_edge>(g, ntask, source_state_l, tasks->sink_l_, C_l);
  cuda_err_chk(cudaGetLastError());
  solver->solve();
  cuda_err_chk(cudaGetLastError());
  return 0;
}


<|start_filename|>include/AR.cpp<|end_filename|>
#include "AR.hpp"
#include <cstdlib>
#include <cstdio>
#include <cmath>

int AR::AutoRegression(
    double *inputseries,
    int length,
    int degree,
    double *coefficients,
    int method) {
  double mean;
  int i, t;
  double *w = nullptr;      /* Input series - mean                            */
  double *h = nullptr;
  double *g = nullptr;      /* Used by mempar()                              */
  double *per = nullptr;
  double *pef = nullptr;      /* Used by mempar()                              */
  double **ar = nullptr;      /* AR coefficients, all degrees                  */

  /* Allocate space for working variables */
  if ((w = (double *) malloc(length * sizeof(double))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }
  if ((h = (double *) malloc((degree + 1) * sizeof(double))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }
  if ((g = (double *) malloc((degree + 2) * sizeof(double))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }
  if ((per = (double *) malloc((length + 1) * sizeof(double))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }
  if ((pef = (double *) malloc((length + 1) * sizeof(double))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }

  if ((ar = (double **) malloc((degree + 1) * sizeof(double *))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }
  for (i = 0; i < degree + 1; i++) {
    if ((ar[i] = (double *) malloc((degree + 1) * sizeof(double))) == nullptr) {
      fprintf(stderr, "Unable to malloc memory - fatal!\n");
      exit(-1);
    }
  }

  /* Determine and subtract the mean from the input series */
  mean = 0.0;
  for (t = 0; t < length; t++)
    mean += inputseries[t];
  mean /= (double) length;
  for (t = 0; t < length; t++)
    w[t] = inputseries[t] - mean;

  /* Perform the appropriate AR calculation */
  if (method == MAXENTROPY) {

    if (!ARMaxEntropy(w, length, degree, ar, per, pef, h, g)) {
      fprintf(stderr, "Max entropy failed - fatal!\n");
      exit(-1);
    }
    for (i = 1; i <= degree; i++)
      coefficients[i - 1] = -ar[degree][i];

  } else if (method == LEASTSQUARES) {

    if (!ARLeastSquare(w, length, degree, coefficients)) {
      fprintf(stderr, "Least squares failed - fatal!\n");
      exit(-1);
    }

  } else {

    fprintf(stderr, "Unknown method\n");
    exit(-1);

  }

  if (w != nullptr)
    free(w);
  if (h != nullptr)
    free(h);
  if (g != nullptr)
    free(g);
  if (per != nullptr)
    free(per);
  if (pef != nullptr)
    free(pef);
  if (ar != nullptr) {
    for (i = 0; i < degree + 1; i++)
      if (ar[i] != nullptr)
        free(ar[i]);
    free(ar);
  }

  return (TRUE);
}
int AR::ARMaxEntropy(
    double *inputseries, int length, int degree, double **ar,
    double *per, double *pef, double *h, double *g) {
  int j, n, nn, jj;
  double sn, sd;
  double t1, t2;

  for (j = 1; j <= length; j++) {
    pef[j] = 0;
    per[j] = 0;
  }

  for (nn = 2; nn <= degree + 1; nn++) {
    n = nn - 2;
    sn = 0.0;
    sd = 0.0;
    jj = length - n - 1;
    for (j = 1; j <= jj; j++) {
      t1 = inputseries[j + n] + pef[j];
      t2 = inputseries[j - 1] + per[j];
      sn -= 2.0 * t1 * t2;
      sd += (t1 * t1) + (t2 * t2);
    }
    g[nn] = sn / sd;
    t1 = g[nn];
    if (n != 0) {
      for (j = 2; j < nn; j++)
        h[j] = g[j] + (t1 * g[n - j + 3]);
      for (j = 2; j < nn; j++)
        g[j] = h[j];
      jj--;
    }
    for (j = 1; j <= jj; j++) {
      per[j] += (t1 * pef[j]) + (t1 * inputseries[j + nn - 2]);
      pef[j] = pef[j + 1] + (t1 * per[j + 1]) + (t1 * inputseries[j]);
    }

    for (j = 2; j <= nn; j++)
      ar[nn - 1][j - 1] = g[j];
  }

  return (TRUE);
}

/*
   Least squares method
   Original from <NAME>, Last modified: Aug 13th, 1998
   Modified (for personal style and context) by <NAME>
*/
int AR::ARLeastSquare(
    double *inputseries,
    int length,
    int degree,
    double *coefficients) {
  int i, j, k, hj, hi;
  double **mat;

  if ((mat = (double **) malloc(degree * sizeof(double *))) == nullptr) {
    fprintf(stderr, "Unable to malloc memory - fatal!\n");
    exit(-1);
  }
  for (i = 0; i < degree; i++) {
    if ((mat[i] = (double *) malloc(degree * sizeof(double))) == nullptr) {
      fprintf(stderr, "Unable to malloc memory - fatal!\n");
      exit(-1);
    }
  }

  for (i = 0; i < degree; i++) {
    coefficients[i] = 0.0;
    for (j = 0; j < degree; j++)
      mat[i][j] = 0.0;
  }
  for (i = degree - 1; i < length - 1; i++) {
    hi = i + 1;
    for (j = 0; j < degree; j++) {
      hj = i - j;
      coefficients[j] += (inputseries[hi] * inputseries[hj]);
      for (k = j; k < degree; k++)
        mat[j][k] += (inputseries[hj] * inputseries[i - k]);
    }
  }
  for (i = 0; i < degree; i++) {
    coefficients[i] /= (length - degree);
    for (j = i; j < degree; j++) {
      mat[i][j] /= (length - degree);
      mat[j][i] = mat[i][j];
    }
  }

  /* Solve the linear equations */
  if (!SolveLE(mat, coefficients, degree)) {
    fprintf(stderr, "Linear solver failed - fatal!\n");
    exit(-1);
  }

  for (i = 0; i < degree; i++)
    if (mat[i] != nullptr)
      free(mat[i]);
  if (mat != nullptr)
    free(mat);
  return (TRUE);
}

/*
   Gaussian elimination solver
   Author: <NAME> Last modified: Aug 14th, 1998
   Modified (for personal style and context) by <NAME>
*/
int AR::SolveLE(double **mat, double *vec, unsigned int n) {
  int i, j, k, maxi;
  double vswap, *mswap, *hvec, max, h, pivot, q;

  for (i = 0; i < n - 1; i++) {
    max = fabs(mat[i][i]);
    maxi = i;
    for (j = i + 1; j < n; j++) {
      if ((h = fabs(mat[j][i])) > max) {
        max = h;
        maxi = j;
      }
    }
    if (maxi != i) {
      mswap = mat[i];
      mat[i] = mat[maxi];
      mat[maxi] = mswap;
      vswap = vec[i];
      vec[i] = vec[maxi];
      vec[maxi] = vswap;
    }

    hvec = mat[i];
    pivot = hvec[i];
    if (fabs(pivot) == 0.0) {
      fprintf(stderr, "Singular matrix - fatal!\n");
      return (FALSE);
    }
    for (j = i + 1; j < n; j++) {
      q = -mat[j][i] / pivot;
      mat[j][i] = 0.0;
      for (k = i + 1; k < n; k++)
        mat[j][k] += q * hvec[k];
      vec[j] += (q * vec[i]);
    }
  }
  vec[n - 1] /= mat[n - 1][n - 1];
  for (i = n - 2; i >= 0; i--) {
    hvec = mat[i];
    for (j = n - 1; j > i; j--)
      vec[i] -= (hvec[j] * vec[j]);
    vec[i] /= hvec[i];
  }

  return (TRUE);
}