/*
MPLG-ECL-GC: This code compresses the input graph and decompresses it on-the-fly to compute Maximal Independent Set.
Copyright 2022 Texas State University

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials
provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written
permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN
IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Authors: Noushin Azami and Martin Burtscher

URL: The latest version of this code is available at
https://cs.txstate.edu/~burtscher/research/MPLG/.

Publication: This work is described in detail in the following paper.
Noushin Azami and Martin Burtscher. Compressed In-memory Graphs for Accelerating GPU-based Analytics. Proceedings of the 12th SC Workshop on Irregular Applications: Architectures and Algorithms. November 2022.
*/


#include <stdlib.h>
#include <stdio.h>
#include <cuda.h>
#include <set>
#include <fstream>
#include <algorithm>
#include "ECLgraph.h"

static const int Device = 0;
static const int ThreadsPerBlock = 256;

typedef unsigned char stattype;
static const stattype in = 0xfe;
static const stattype out = 0;

/* main computation kernel */

template <int logn>
static __global__ void findmins(const int nodes, const int* const __restrict__ nidx, volatile stattype* const __restrict__ nstat, const int* const __restrict__ encoded)
{
  const int from = threadIdx.x + blockIdx.x * ThreadsPerBlock;
  const int incr = gridDim.x * ThreadsPerBlock;

  int missing;
  do {
    missing = 0;
    for (int v = from; v < nodes; v += incr) {
      const stattype nv = nstat[v];
      if (nv & 1) {
        const int end = nidx[v + 1];
        const int beg = nidx[v];
        long cur = (long)beg * logn;
        int shift = cur % 32;
        int i = nidx[v];
        int pos = cur / 32;
        int res = __funnelshift_rc(encoded[pos], encoded[pos + 1], shift);
        int nl = res & ((1 << logn) - 1);
        while ((i < end) && ((nv > nstat[nl]) || ((nv == nstat[nl]) && (v > nl)))) {
          i++;
          cur += logn;
          shift = (shift + logn) % 32;
          pos = cur / 32;
          res = __funnelshift_rc(encoded[pos], encoded[pos + 1], shift);
          nl = res & ((1 << logn) - 1);
        }
        if (i < end) {
          missing = 1;
        } else {
          cur = (long)beg * logn;
          shift = cur % 32;
          for (int i = beg; i < end; i++) {
            const int pos = cur / 32;
            int res = __funnelshift_rc(encoded[pos], encoded[pos + 1], shift);
            const int nli = res & ((1 << logn) - 1);
            nstat[nli] = out;
            cur += logn;
            shift = (shift + logn) % 32;
          }
          nstat[v] = in;
        }
      }
    }
  } while (missing != 0);
}

/* hash function to generate random values */

// source of hash function: https://stackoverflow.com/questions/664014/what-integer-hash-function-are-good-that-accepts-an-integer-hash-key
static __device__ unsigned int hash(unsigned int val)
{
  val = ((val >> 16) ^ val) * 0x45d9f3b;
  val = ((val >> 16) ^ val) * 0x45d9f3b;
  return (val >> 16) ^ val;
}

/* prioritized-selection initialization kernel */

static __global__ void init(const int nodes, const int edges, const int* const __restrict__ nidx, stattype* const __restrict__ nstat)
{
  const int from = threadIdx.x + blockIdx.x * ThreadsPerBlock;
  const int incr = gridDim.x * ThreadsPerBlock;

  const float avg = (float)edges / nodes;
  const float scaledavg = ((in / 2) - 1) * avg;

  for (int i = from; i < nodes; i += incr) {
    stattype val = in;
    const int degree = nidx[i + 1] - nidx[i];
    if (degree > 0) {
      float x = degree - (hash(i) * 0.00000000023283064365386962890625f);
      int res = __float2int_rn(scaledavg / (avg + x));
      val = (res + res) | 1;
    }
    nstat[i] = val;
  }
}

struct GPUTimer
{
  cudaEvent_t beg, end;
  GPUTimer() {cudaEventCreate(&beg);  cudaEventCreate(&end);}
  ~GPUTimer() {cudaEventDestroy(beg);  cudaEventDestroy(end);}
  void start() {cudaEventRecord(beg, 0);}
  float stop() {cudaEventRecord(end, 0);  cudaEventSynchronize(end);  float ms;  cudaEventElapsedTime(&ms, beg, end);  return 0.001f * ms;}
};

static void computeMIS(const int nodes, const int edges, const int* const __restrict__ nidx, stattype* const __restrict__ nstat, const int* const __restrict__ encoded, const int logn, const int elems)
{
  cudaSetDevice(Device);
  cudaDeviceProp deviceProp;
  cudaGetDeviceProperties(&deviceProp, Device);
  if ((deviceProp.major == 9999) && (deviceProp.minor == 9999)) {fprintf(stderr, "ERROR: there is no CUDA capable device\n\n");  exit(-1);}
  const int SMs = deviceProp.multiProcessorCount;
  const int mTpSM = deviceProp.maxThreadsPerMultiProcessor;
  printf("gpu: %s with %d SMs and %d mTpSM (%.1f MHz core and %.1f MHz mem)\n", deviceProp.name, SMs, mTpSM, deviceProp.clockRate * 0.001, deviceProp.memoryClockRate * 0.001);

  int* nidx_d;
  stattype* nstat_d;
  int* encoded_d;

  if (cudaSuccess != cudaMalloc((void **)&encoded_d, elems * sizeof(int))) {fprintf(stderr, "ERROR: could not allocate nlist_d\n\n");  exit(-1);}
  if (cudaSuccess != cudaMalloc((void **)&nidx_d, (nodes + 1) * sizeof(int))) {fprintf(stderr, "ERROR: could not allocate nidx_d\n\n");  exit(-1);}
  if (cudaSuccess != cudaMalloc((void **)&nstat_d, nodes * sizeof(stattype))) {fprintf(stderr, "ERROR: could not allocate nstat_d\n\n");  exit(-1);}

  if (cudaSuccess != cudaMemcpy(nidx_d, nidx, (nodes + 1) * sizeof(int), cudaMemcpyHostToDevice)) {fprintf(stderr, "ERROR: copying to device failed\n\n");  exit(-1);}
  if (cudaSuccess != cudaMemcpy(encoded_d, encoded, elems * sizeof(int), cudaMemcpyHostToDevice)) {fprintf(stderr, "ERROR: copying to device failed\n\n");  exit(-1);}

  const int blocks = SMs * mTpSM / ThreadsPerBlock;
  GPUTimer timer;
  timer.start();
  init<<<blocks, ThreadsPerBlock>>>(nodes, edges, nidx_d, nstat_d);
  switch (logn) {
    case 16: findmins<16><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 17: findmins<17><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 18: findmins<18><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 19: findmins<19><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 20: findmins<20><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 21: findmins<21><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 22: findmins<22><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 23: findmins<23><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 24: findmins<24><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 25: findmins<25><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 26: findmins<26><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 27: findmins<27><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
    case 28: findmins<28><<<blocks, ThreadsPerBlock>>>(nodes, nidx_d, nstat_d, encoded_d);
  }
  float runtime = timer.stop();

  printf("compute time: %.6f s\n", runtime);
  printf("throughput: %.6f Mnodes/s\n", nodes * 0.000001 / runtime);
  printf("throughput: %.6f Medges/s\n", edges * 0.000001 / runtime);

  if (cudaSuccess != cudaMemcpy(nstat, nstat_d, nodes * sizeof(stattype), cudaMemcpyDeviceToHost)) {fprintf(stderr, "ERROR: copying from device failed\n\n");  exit(-1);}

  cudaFree(nstat_d);
  cudaFree(encoded_d);
  cudaFree(nidx_d);
}

int main(int argc, char* argv[])
{
  printf("MPLG-ECL-MIS v1 (%s)\n", __FILE__);
  printf("Copyright 2022 Texas State University\n");

  if (argc != 2) {fprintf(stderr, "USAGE: %s input_file_name\n\n", argv[0]);  exit(-1);}

  ECLgraph g = readECLgraph(argv[1]);
  printf("configuration: %d nodes and %d edges (%s)\n", g.nodes, g.edges, argv[1]);
  printf("average degree: %.2f edges per node\n", 1.0 * g.edges / g.nodes);

  stattype* nstatus = (stattype*)malloc(g.nodes * sizeof(nstatus[0]));
  if (nstatus == NULL) {fprintf(stderr, "ERROR: could not allocate nstatus\n\n");  exit(-1);}

  //encode
  int* const encoded = new int [g.edges];
  for (int i = 0; i < g.edges; i++) encoded[i] = 0;
  const int logn = 32 - __builtin_clz(g.nodes);
  long loc = 0;
  for (int j = 0; j < g.edges; j++) {
    const int val = g.nlist[j];
    const int pos = loc / 32;
    const int shift = loc % 32;
    encoded[pos] |= val << shift;
    if (32 - logn < shift) {
      encoded[pos + 1] = (unsigned int)val >> (32 - shift);
    }
    loc += logn;
  }
  const int elems = (loc + 31) / 32;

  computeMIS(g.nodes, g.edges, g.nindex, nstatus, encoded, logn, elems);

  int count = 0;
  for (int v = 0; v < g.nodes; v++) {
    if (nstatus[v] == in) {
      count++;
    }
  }
  printf("elements in set: %d (%.1f%%)\n", count, 100.0 * count / g.nodes);

  /* result verification code */

  for (int v = 0; v < g.nodes; v++) {
    if ((nstatus[v] != in) && (nstatus[v] != out)) {fprintf(stderr, "ERROR: found unprocessed node in graph\n\n");  exit(-1);}
    if (nstatus[v] == in) {
      for (int i = g.nindex[v]; i < g.nindex[v + 1]; i++) {
        if (nstatus[g.nlist[i]] == in) {fprintf(stderr, "ERROR: found adjacent nodes in MIS\n\n");  exit(-1);}
      }
    } else {
      int flag = 0;
      for (int i = g.nindex[v]; i < g.nindex[v + 1]; i++) {
        if (nstatus[g.nlist[i]] == in) {
          flag = 1;
        }
      }
      if (flag == 0) {fprintf(stderr, "ERROR: set is not maximal\n\n");  exit(-1);}
    }
  }

  freeECLgraph(g);
  free(nstatus);
  return 0;
}