text: 
  #include <CL/sycl.hpp>
  #include <dpct/dpct.hpp>
  #include <math.h>
  #include <float.h>
  #include <fstream>
  #include <string>
  #include <cstring>
  #include <sstream>
  #include <iostream>
  #include <chrono>
  
  using namespace std::chrono;
  typedef high_resolution_clock myclock;
  typedef duration<float> myduration;
  
  #define MAX_WG_SIZE 256
  
  template <typename T>
  T* mem_alloc (const int align, const size_t size) {
    return (T*) aligned_alloc(align, size * sizeof(T));
  }
  
  template <typename T>
  void mem_free (T* p) {
    free(p);
  }
  
  
  float gammafunction(unsigned int n)
  {   
    if(n == 0)
      return 0.0f;
    float x = ((float)n + 0.5f) * sycl::log((float)n) -
              ((float)n - 1.0f) * sycl::log(sycl::exp(1.0f));
    return x;
  }
  
  void epi(const unsigned int* dev_data_zeros, 
                      const unsigned int* dev_data_ones, 
                      float* dev_scores, 
                      const int num_snp, 
                      const int PP_zeros, 
                      const int PP_ones,
                      const int mask_zeros, 
                      const int mask_ones,
                      sycl::nd_item<3> item_ct1) 
  {
  
    int i, j, tid, p, k;
    float score = FLT_MAX;
  
    j = item_ct1.get_local_range().get(2) * item_ct1.get_group(2) +
        item_ct1.get_local_id(2);
    i = item_ct1.get_local_range().get(1) * item_ct1.get_group(1) +
        item_ct1.get_local_id(1);
    tid = i * num_snp + j;
  
    if (j > i && i < num_snp && j < num_snp) {
      unsigned int ft[2 * 9];
      for(k = 0; k < 2 * 9; k++) ft[k] = 0;
  
      unsigned int t00, t01, t02, t10, t11, t12, t20, t21, t22;
      unsigned int di2, dj2;
      unsigned int* SNPi;
      unsigned int* SNPj;
  
      // Phenotype 0
      SNPi = (unsigned int*) &dev_data_zeros[i * 2];
      SNPj = (unsigned int*) &dev_data_zeros[j * 2];
      for (p = 0; p < 2 * PP_zeros * num_snp - 2 * num_snp; p += 2 * num_snp) {
        di2 = ~(SNPi[p] | SNPi[p + 1]);
        dj2 = ~(SNPj[p] | SNPj[p + 1]);
  
        t00 = SNPi[p] & SNPj[p];
        t01 = SNPi[p] & SNPj[p + 1];
        t02 = SNPi[p] & dj2;
        t10 = SNPi[p + 1] & SNPj[p];
        t11 = SNPi[p + 1] & SNPj[p + 1];
        t12 = SNPi[p + 1] & dj2;
        t20 = di2 & SNPj[p];
        t21 = di2 & SNPj[p + 1];
        t22 = di2 & dj2;
  
        ft[0] += sycl::popcount(t00);
        ft[1] += sycl::popcount(t01);
        ft[2] += sycl::popcount(t02);
        ft[3] += sycl::popcount(t10);
        ft[4] += sycl::popcount(t11);
        ft[5] += sycl::popcount(t12);
        ft[6] += sycl::popcount(t20);
        ft[7] += sycl::popcount(t21);
        ft[8] += sycl::popcount(t22);
      }
  
      // remainder
      p = 2 * PP_zeros * num_snp - 2 * num_snp;
      di2 = ~(SNPi[p] | SNPi[p + 1]);
      dj2 = ~(SNPj[p] | SNPj[p + 1]);
      di2 = di2 & mask_zeros;
      dj2 = dj2 & mask_zeros;
  
      t00 = SNPi[p] & SNPj[p];
      t01 = SNPi[p] & SNPj[p + 1];
      t02 = SNPi[p] & dj2;
      t10 = SNPi[p + 1] & SNPj[p];
      t11 = SNPi[p + 1] & SNPj[p + 1];
      t12 = SNPi[p + 1] & dj2;
      t20 = di2 & SNPj[p];
      t21 = di2 & SNPj[p + 1];
      t22 = di2 & dj2;
  
      ft[0] += sycl::popcount(t00);
      ft[1] += sycl::popcount(t01);
      ft[2] += sycl::popcount(t02);
      ft[3] += sycl::popcount(t10);
      ft[4] += sycl::popcount(t11);
      ft[5] += sycl::popcount(t12);
      ft[6] += sycl::popcount(t20);
      ft[7] += sycl::popcount(t21);
      ft[8] += sycl::popcount(t22);
  
      // Phenotype 1
      SNPi = (unsigned int*) &dev_data_ones[i * 2];
      SNPj = (unsigned int*) &dev_data_ones[j * 2];
      for(p = 0; p < 2 * PP_ones * num_snp - 2 * num_snp; p += 2 * num_snp)
      {
        di2 = ~(SNPi[p] | SNPi[p + 1]);
        dj2 = ~(SNPj[p] | SNPj[p + 1]);
  
        t00 = SNPi[p] & SNPj[p];
        t01 = SNPi[p] & SNPj[p + 1];
        t02 = SNPi[p] & dj2;
        t10 = SNPi[p + 1] & SNPj[p];
        t11 = SNPi[p + 1] & SNPj[p + 1];
        t12 = SNPi[p + 1] & dj2;
        t20 = di2 & SNPj[p];
        t21 = di2 & SNPj[p + 1];
        t22 = di2 & dj2;
  
        ft[9] += sycl::popcount(t00);
        ft[10] += sycl::popcount(t01);
        ft[11] += sycl::popcount(t02);
        ft[12] += sycl::popcount(t10);
        ft[13] += sycl::popcount(t11);
        ft[14] += sycl::popcount(t12);
        ft[15] += sycl::popcount(t20);
        ft[16] += sycl::popcount(t21);
        ft[17] += sycl::popcount(t22);
      }
      p = 2 * PP_ones * num_snp - 2 * num_snp;
      di2 = ~(SNPi[p] | SNPi[p + 1]);
      dj2 = ~(SNPj[p] | SNPj[p + 1]);
      di2 = di2 & mask_ones;
      dj2 = dj2 & mask_ones;
  
      t00 = SNPi[p] & SNPj[p];
      t01 = SNPi[p] & SNPj[p + 1];
      t02 = SNPi[p] & dj2;
      t10 = SNPi[p + 1] & SNPj[p];
      t11 = SNPi[p + 1] & SNPj[p + 1];
      t12 = SNPi[p + 1] & dj2;
      t20 = di2 & SNPj[p];
      t21 = di2 & SNPj[p + 1];
      t22 = di2 & dj2;
  
      ft[9] += sycl::popcount(t00);
      ft[10] += sycl::popcount(t01);
      ft[11] += sycl::popcount(t02);
      ft[12] += sycl::popcount(t10);
      ft[13] += sycl::popcount(t11);
      ft[14] += sycl::popcount(t12);
      ft[15] += sycl::popcount(t20);
      ft[16] += sycl::popcount(t21);
      ft[17] += sycl::popcount(t22);
  
      // compute score
      score = 0.0f;
      for(k = 0; k < 9; k++)
        score += gammafunction(ft[k] + ft[9 + k] + 1) - gammafunction(ft[k]) - gammafunction(ft[9 + k]);
      score = sycl::fabs(score);
      if(score == 0.0f)
        score = FLT_MAX;
      dev_scores[tid] = score;
    }
  }
  
  
  int main(int argc, char **argv)
  {
    dpct::device_ext &dev_ct1 = dpct::get_current_device();
    sycl::queue &q_ct1 = dev_ct1.default_queue();
    int i, j, x;
    int num_pac = atoi(argv[1]);  // #samples
    int num_snp = atoi(argv[2]);  // #SNPs
    int iteration = atoi(argv[3]);// #kernel run
    int block_snp = 64;
  
    srand(100);
    unsigned char *SNP_Data = mem_alloc<unsigned char>(64, num_pac * num_snp);
    unsigned char *Ph_Data = mem_alloc<unsigned char>(64, num_pac);
  
    // generate SNPs between 0 and 2
    for (i = 0; i < num_pac; i++)
      for(j = 0; j < num_snp; j++)
        SNP_Data[i * num_snp + j] = rand() % 3;
  
    // generate phenotype between 0 and 1
    for(int i = 0; i < num_pac; i++) Ph_Data[i] = rand() % 2;
  
    // transpose the SNP data
    unsigned char *SNP_Data_trans = mem_alloc<unsigned char>(64, num_pac * num_snp);
  
    for (i = 0; i < num_pac; i++) 
      for(j = 0; j < num_snp; j++) 
        SNP_Data_trans[j * num_pac + i] = SNP_Data[i * num_snp + j];
  
    int phen_ones = 0;
    for(i = 0; i < num_pac; i++)
      if(Ph_Data[i] == 1)
        phen_ones++;
  
    // transform SNP data to a binary format
  
    int PP_zeros = ceil((1.0*(num_pac - phen_ones))/32.0);
    int PP_ones = ceil((1.0*phen_ones)/32.0);
  
    unsigned int *bin_data_zeros = mem_alloc<unsigned int>(64, num_snp * PP_zeros * 2);
    unsigned int *bin_data_ones = mem_alloc<unsigned int>(64, num_snp * PP_ones * 2);
    memset(bin_data_zeros, 0, num_snp*PP_zeros*2*sizeof(unsigned int));
    memset(bin_data_ones, 0, num_snp*PP_ones*2*sizeof(unsigned int));
  
    for(i = 0; i < num_snp; i++)
    {
      int x_zeros = -1;
      int x_ones = -1;
      int n_zeros = 0;
      int n_ones = 0;
  
      for(j = 0; j < num_pac; j++){
        unsigned int temp = (unsigned int) SNP_Data_trans[i * num_pac + j];
  
        if(Ph_Data[j] == 1){
          if(n_ones%32 == 0){
            x_ones ++;
          }
          // apply 1 shift left to 2 components
          bin_data_ones[i * PP_ones * 2 + x_ones*2 + 0] <<= 1;
          bin_data_ones[i * PP_ones * 2 + x_ones*2 + 1] <<= 1;
          // insert '1' in correct component
          if(temp == 0 || temp == 1){
            bin_data_ones[i * PP_ones * 2 + x_ones*2 + temp ] |= 1;
          }
          n_ones ++;
        } else {
          if(n_zeros%32 == 0){
            x_zeros ++;
          }
          // apply 1 shift left to 2 components
          bin_data_zeros[i * PP_zeros * 2 + x_zeros*2 + 0] <<= 1;
          bin_data_zeros[i * PP_zeros * 2 + x_zeros*2 + 1] <<= 1;
          // insert '1' in correct component
          if(temp == 0 || temp == 1){
            bin_data_zeros[i * PP_zeros * 2 + x_zeros*2 + temp] |= 1;
          }
          n_zeros ++;
        }
      }
    }
  
    unsigned int mask_zeros = 0xFFFFFFFF;
    for(int x = num_pac - phen_ones; x < PP_zeros * 32; x++)
      mask_zeros = mask_zeros >> 1;
  
    unsigned int mask_ones = 0xFFFFFFFF;
    for(x = phen_ones; x < PP_ones * 32; x++)
      mask_ones = mask_ones >> 1;
  
    // transpose the binary data structures
    unsigned int* bin_data_ones_trans = mem_alloc<unsigned int>(64, num_snp * PP_ones * 2);
  
    for(i = 0; i < num_snp; i++)
      for(j = 0; j < PP_ones; j++)
      {
        bin_data_ones_trans[(j * num_snp + i) * 2 + 0] = bin_data_ones[(i * PP_ones + j) * 2 + 0];
        bin_data_ones_trans[(j * num_snp + i) * 2 + 1] = bin_data_ones[(i * PP_ones + j) * 2 + 1];
      }
  
    unsigned int* bin_data_zeros_trans = mem_alloc<unsigned int>(64, num_snp * PP_zeros * 2);
  
    for(i = 0; i < num_snp; i++)
      for(j = 0; j < PP_zeros; j++)
      {
        bin_data_zeros_trans[(j * num_snp + i) * 2 + 0] = bin_data_zeros[(i * PP_zeros + j) * 2 + 0];
        bin_data_zeros_trans[(j * num_snp + i) * 2 + 1] = bin_data_zeros[(i * PP_zeros + j) * 2 + 1];
      }
  
    float *scores = mem_alloc<float>(64, num_snp * num_snp);
    for(x = 0; x < num_snp * num_snp; x++) scores[x] = FLT_MAX;
  
    auto start = myclock::now();
  
    unsigned int* d_data_zeros;
    d_data_zeros =
        sycl::malloc_device<unsigned int>(num_snp * PP_zeros * 2, q_ct1);
    q_ct1
        .memcpy(d_data_zeros, bin_data_zeros_trans,
                num_snp * PP_zeros * 2 * sizeof(unsigned int))
        .wait();
  
    unsigned int* d_data_ones;
    d_data_ones = sycl::malloc_device<unsigned int>(num_snp * PP_ones * 2, q_ct1);
    q_ct1
        .memcpy(d_data_ones, bin_data_ones_trans,
                num_snp * PP_ones * 2 * sizeof(unsigned int))
        .wait();
  
    float* d_scores;
    d_scores = sycl::malloc_device<float>(num_snp * num_snp, q_ct1);
  
    // setup kernel ND-range
    int num_snp_m = num_snp;
    while(num_snp_m % block_snp != 0) num_snp_m++;
  
    sycl::range<3> grid(1, num_snp_m, num_snp_m / block_snp);
    sycl::range<3> block(1, 1, block_snp);
  
    // epistasis detection kernel
    for (int i = 0; i < iteration; i++) {
      q_ct1.memcpy(d_scores, scores, sizeof(float) * num_snp * num_snp).wait();
      /*
      DPCT1049:0: The workgroup size passed to the SYCL kernel may exceed
       * the limit. To get the device limit, query
       * info::device::max_work_group_size. Adjust the workgroup size if needed.
  
       */
      q_ct1.submit([&](sycl::handler &cgh) {
        cgh.parallel_for(sycl::nd_range<3>(grid * block, block),
                         [=](sycl::nd_item<3> item_ct1) {
                           epi(d_data_zeros, d_data_ones, d_scores, num_snp,
                               PP_zeros, PP_ones, mask_zeros, mask_ones,
                               item_ct1);
                         });
      });
    }
    q_ct1.memcpy(scores, d_scores, sizeof(float) * num_snp * num_snp).wait();
  
    auto end = myclock::now();
    myduration elapsed = end - start;
    std::cout << "Total offloading time: " << elapsed.count() << " sec" << std::endl;
  
    // compute the minimum score on a host
    float score = scores[0];
    int solution = 0;
    for (int i = 1; i < num_snp * num_snp; i++) {
      if (score > scores[i]) {
        score = scores[i];
        solution = i;
      }
    }
  
    std::cout << "Score: " << score << std::endl;
    std::cout << "Solution: " << solution / num_snp << ", " << solution % num_snp << std::endl;
    if ( (fabsf(score - 83.844f) > 1e-3f) || (solution / num_snp != 1253) || 
        (solution % num_snp != 25752) )
      std::cout << "FAIL\n";
    else
      std::cout << "PASS\n";
  
    sycl::free(d_data_zeros, q_ct1);
    sycl::free(d_data_ones, q_ct1);
    sycl::free(d_scores, q_ct1);
  
    mem_free(bin_data_zeros);
    mem_free(bin_data_ones);
    mem_free(bin_data_zeros_trans);
    mem_free(bin_data_ones_trans);
    mem_free(scores);
    mem_free(SNP_Data);
    mem_free(SNP_Data_trans);
    mem_free(Ph_Data);
    return 0;
  }
  