/* NOTE: ADAPTED and MODIFIED by Hoang-Trong Minh Tuan @2010 George Mason University  */

/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * NVIDIA Corporation and its licensors retain all intellectual property and 
 * proprietary rights in and to this software and related documentation. 
 * Any use, reproduction, disclosure, or distribution of this software 
 * and related documentation without an express license agreement from
 * NVIDIA Corporation is strictly prohibited.
 *
 * Please refer to the applicable NVIDIA end user license agreement (EULA) 
 * associated with this source code for terms and conditions that govern 
 * your use of this NVIDIA software.
 * 
 */ 

/*
    Parallel reduction

    This sample shows how to perform a reduction operation on an array of values
    to produce a single value.

    Reductions are a very common computation in parallel algorithms.  Any time
    an array of values needs to be reduced to a single value using a binary 
    associative operator, a reduction can be used.  Example applications include
    statistics computaions such as mean and standard deviation, and image 
    processing applications such as finding the total luminance of an
    image.

    This code performs sum reductions, but any associative operator such as
    min() or max() could also be used.

    It assumes the input size is a power of 2.

    COMMAND LINE ARGUMENTS

    "--shmoo":         Test performance for 1 to 32M elements with each of the 7 different kernels
    "--n=<N>":         Specify the number of elements to reduce (default 1048576)
    "--threads=<N>":   Specify the number of threads per block (default 128)
    "--kernel=<N>":    Specify which kernel to run (0-6, default 6)
    "--maxblocks=<N>": Specify the maximum number of thread blocks to launch (kernel 6 only, default 64)
    "--cpufinal":      Read back the per-block results and do final sum of block sums on CPU (default false)
    "--cputhresh=<N>": The threshold of number of blocks sums below which to perform a CPU final reduction (default 1)
    "-type=<T>":       The datatype for the reduction, where T is "int", "float", or "double" (default int)
*/

// Utilities and system includes
#include <algorithm>
#include <cuda.h>
#include "cutil_inline.h"

// includes, project
#include "reduction.h"


////////////////////////////////////////////////////////////////////////////////
// declaration, forward
template <class T>
T*
sum_reduce(T* d_idata, int size, int maxThreads=256, int whichKernel=6, int maxBlocks=64);

#define MAX_BLOCK_DIM_SIZE 65535

#ifdef WIN32
#define strcasecmp strcmpi
#endif

extern "C"
bool isPow2(unsigned int x)
{
    return ((x&(x-1))==0);
}

unsigned int nextPow2( unsigned int x ) {
    --x;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;
    return ++x;
}

#define MIN(x,y) ((x < y) ? x : y)

////////////////////////////////////////////////////////////////////////////////
// Compute the number of threads and blocks to use for the given reduction kernel
// For the kernels >= 3, we set threads / block to the minimum of maxThreads and
// n/2. For kernels < 3, we set to the minimum of maxThreads and n.  For kernel 
// 6, we observe the maximum specified number of blocks, because each thread in 
// that kernel can process a variable number of elements.
////////////////////////////////////////////////////////////////////////////////
void getNumBlocksAndThreads(int whichKernel, int n, int maxBlocks, int maxThreads, int &blocks, int &threads)
{
    
    
    
    if (whichKernel < 3)
    {
        threads = (n < maxThreads) ? nextPow2(n) : maxThreads;
        blocks = (n + threads - 1) / threads;
    }
    else
    {
        threads = (n < maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
        blocks = (n + (threads * 2 - 1)) / (threads * 2);
    }
        

    if (whichKernel == 6)
        blocks = MIN(maxBlocks, blocks);
}
void getNumBlocksAndThreads_6(int n, int maxBlocks, int maxThreads, int &blocks, int &threads)
{
    threads = (n < maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
    blocks = (n + (threads * 2 - 1)) / (threads * 2);
    blocks = MIN(maxBlocks, blocks);
}


////////////////////////////////////////////////////////////////////////////////////
//TUAN: this section of code hereafter is developed by Tuan, Hoang-Trong
////////////////////////////////////////////////////////////////////////////////////
#define NUMBLOCK_TESLA 64
#define NUMBLOCK_FERMI 64
#define MAXBLOCK_TESLA 64
#define MAXBLOCK_FERMI 64
template <class T>
T*  init_fast_tesla() {
    T * d_odata = NULL;
    cudaMalloc((void**) &d_odata, NUMBLOCK_TESLA*sizeof(T));
    return d_odata;
}
template <class T>
T* init_fast_fermi() {
    T* d_odata = NULL;
    cudaMalloc((void**) &d_odata, NUMBLOCK_FERMI*sizeof(T));
    return d_odata;
}
template <class T>
void finalize_fast_tesla(T * d_odata){
    cudaFree(d_odata);
}
template <class T>
void finalize_fast_fermi(T * d_odata){
    cudaFree(d_odata);
}

//forward
template <class T>
T
sum_reduce6_fast_tesla(T* d_idata, int size, T* d_odata);

template <class T>
T* 
sum_reduce6_fast_tesla_test(T* d_idata, int size, T* d_odata);

template <class T>
T
sum_reduce6_fast_fermi(T* d_idata, int size);

/*
  Call this function only once */
extern "C" double* init_double_tesla(){
    return init_fast_tesla<double>();
}
extern "C" float* init_float_tesla(){
    return init_fast_tesla<float>();
}
extern "C" int* init_int_tesla(){
    return init_fast_tesla<int>();
}
/*
  Call this function only once */
extern "C" double* init_double_fermi(){
    return init_fast_fermi<double>();
}
extern "C" float*  init_float_fermi(){
    return init_fast_fermi<float>();
}
extern "C" int* init_int_fermi(){
    return init_fast_fermi<int>();
}

/*
  Call this function only once */
extern "C"  void finalize_double_tesla(double * d_idata){
    finalize_fast_tesla<double>(d_idata);
}
extern "C"  void finalize_float_tesla(float * d_idata){
    finalize_fast_tesla<float>(d_idata);
}
extern "C"  void finalize_int_tesla(int * d_idata){
    finalize_fast_tesla<int>(d_idata);
}
/*
  Call this function only once */
extern "C"  void finalize_double_fermi(double * d_idata){
    finalize_fast_fermi<double>(d_idata);
}
extern "C"  void finalize_float_fermi(float * d_idata){
    finalize_fast_fermi<float>(d_idata);
}
extern "C"  void finalize_int_fermi(int * d_idata){
    finalize_fast_fermi<int>(d_idata);
}

// NON-working version (BUG)
extern "C" double sum_reduce_double_2host_tesla(double* d_idata, int n, double *d_odata) {
    int numBlocks= NUMBLOCK_TESLA;
    cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(double), cudaMemcpyDeviceToDevice);
    double sumval;
    sumval = sum_reduce6_fast_tesla<double>(d_idata, n, d_odata);
    return sumval;
}
// WORKING version
extern "C" double sum_reduce_double_2host_tesla_test(double* d_idata, int n, double *d_odata) {
    int numBlocks= NUMBLOCK_TESLA;
    cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(double), cudaMemcpyDeviceToDevice);
    double *res = sum_reduce6_fast_tesla_test<double>(d_idata, n, d_odata);
    double sumval;
    cudaMemcpy(&sumval, d_odata, sizeof(double), cudaMemcpyDeviceToHost);
    return sumval;
}

extern "C" double sum_reduce_double_2host_fermi(double* p_arr, int n) {
    double sumval;
    sumval = sum_reduce6_fast_fermi<double>(p_arr, n);
    return sumval;
}

double* sum_reduce_double(double* p_arr, int n, int maxThreads=256, int whichKernel=6, int maxBlocks=64) {
    return sum_reduce<double>(p_arr, n, maxThreads, whichKernel, maxBlocks);
}

int* sum_reduce_int(int* p_arr, int n, int maxThreads=256, int whichKernel=6, int maxBlocks=64) {
    return sum_reduce<int>(p_arr, n, maxThreads, whichKernel, maxBlocks);
}

float* sum_reduce_float(float* p_arr, int n, int maxThreads=256, int whichKernel=6, int maxBlocks=64) {
    return sum_reduce<float>(p_arr, n, maxThreads, whichKernel, maxBlocks);
}


/*
   This function is designed to work fast, all parameters are fixed, so less flexible
   You must call the function init_fast_tesla(...) first and only once
*/   
template <class T>
T 
sum_reduce6_fast_tesla(T* d_idata, int size, T* d_odata)
{
    // //use 256 for FERMI
    // // and 128 for Tesla
    // int maxThreads =128;
    // int numThreads = 128;
    // int whichKernel=6;
    // int maxBlocks= MAXBLOCK_TESLA;
    // int numBlocks= NUMBLOCK_TESLA;

    // /* move to sum_reduce_double_host_tesla */
    // // cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(T), cudaMemcpyDeviceToDevice);


    // // execute the kernel
    // reduce6_sum<T>(size, numThreads, numBlocks, d_idata, d_odata);
    // //sum partial block sums on GPU
    // int s=numBlocks;
    // int kernel = whichKernel;

    // // while(s > cpuFinalThreshold) {
    // // while(s > 1) {
    // //     int threads = 0, blocks = 0;
    // //     getNumBlocksAndThreads(kernel, s, maxBlocks, maxThreads, blocks, threads);
                
    // //     reduce_sum<T>(s, threads, blocks, kernel, d_odata, d_odata);
                
    // //     if (kernel < 3)
    // //         s = (s + threads - 1) / threads;
    // //     else
    // //         s = (s + (threads*2-1)) / (threads*2);
    // // }
    // // sum partial sums from each block on CPU        
    // // copy result from device to host
    // T h_odata[NUMBLOCK_TESLA];
    // cutilSafeCallNoSync( cudaMemcpy( &h_odata, d_odata, NUMBLOCK_TESLA*sizeof(T), cudaMemcpyDeviceToHost) );
    T result=0;
    // // cudaMemcpy(&result, d_odata, sizeof(T), cudaMemcpyDeviceToHost);
    // for(int i=0; i<numBlocks; i++) {
    //     result += h_odata[i];
    // }

    return result; 
}

/// Use kernel type 6
template <class T>
T* 
sum_reduce6_fast_tesla_test(T* d_idata, int size, T* d_odata)
{
    //use 256 for FERMI
    // and 128 for Tesla
    int maxThreads =128;
    int numThreads = 128;
    int maxBlocks= MAXBLOCK_TESLA;
    int numBlocks= NUMBLOCK_TESLA;

    /* moved to sum_reduce_double_host_tesla_test */
    // cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(T), cudaMemcpyDeviceToDevice);

    // execute the kernel
    reduce6_sum<T>(size, numThreads, numBlocks, d_idata, d_odata);

    //sum partial block sums on GPU
    int s=numBlocks;

    // while(s > cpuFinalThreshold) {
    while(s > 1) {
        int threads = 0, blocks = 0;
        getNumBlocksAndThreads_6(s, maxBlocks, maxThreads, blocks, threads);
                
        reduce6_sum<T>(s, threads, blocks, d_odata, d_odata);
                
        s = (s + (threads*2-1)) / (threads*2);
    }

    return d_odata; 
}

/// User kernel type 6
template <class T>
T
sum_reduce6_fast_fermi(T* d_idata, int size)
{
    //use 256 for FERMI
    // and 192 for Tesla
    // int maxThreads =256;
    // int numThreads = 256;
    // int maxBlocks=64 ;
    // int numBlocks                         = 64;
    // T* d_odata = NULL;

    // cudaMalloc((void**) &d_odata, numBlocks*sizeof(T));
    // cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(T), cudaMemcpyDeviceToDevice);

    // // execute the kernel
    // reduce_sum<T>(size, numThreads, numBlocks, d_idata, d_odata);
    // //sum partial block sums on GPU
    // int s=numBlocks;

    // // while(s > cpuFinalThreshold) {
    // while(s > 1) {
    //     int threads = 0, blocks = 0;
    //     getNumBlocksAndThreads(kernel, s, maxBlocks, maxThreads, blocks, threads);
                
    //     reduce6_sum<T>(s, threads, blocks, d_odata, d_odata);
                
    //     s = (s + (threads*2-1)) / (threads*2);
    // }

    T result;
    // cudaMemcpy(&result, d_odata, sizeof(T), cudaMemcpyDeviceToHost);
    // cudaFree(d_odata);
    return result; 
}


////////////////////////////////////////////////////////////////////////////////
// The main function which runs the reduction test.
/* size : number of elements to reduce (should be power-of-2)
   d_idata : pointer to array of (size) elements in the device memory
   maxThreads : maximum number of threads per block (multiple of warp-size)
   whichKernel : choose the kernel version (0..6)
   maxBlocks : maximum number of thread blocks to launch (work only whichKernel=6)

   NOTE: input data must be on device,
   input data must be continuous
   output data (the sum reduction) still on device
         
*/
////////////////////////////////////////////////////////////////////////////////
template <class T>
T* 
sum_reduce(T* d_idata, int size, int maxThreads=256, int whichKernel=6, int maxBlocks=64)
{

    int numBlocks                         = 0;
    int numThreads                        = 0;
    T* result;
    // int cpuFinalThreshold = 1;
    
    // detect the number of threads per block, and number of blocks per grid
    getNumBlocksAndThreads(whichKernel, size, maxBlocks, maxThreads, numBlocks, numThreads);
    // if (numBlocks == 1) cpuFinalThreshold = 1;

    
    // allocate device memory and data
    T* d_odata = NULL;
    cudaMalloc((void**) &d_odata, numBlocks*sizeof(T));
    cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(T), cudaMemcpyDeviceToDevice);

    // execute the kernel
    reduce_sum<T>(size, numThreads, numBlocks, whichKernel, d_idata, d_odata);
    //sum partial block sums on GPU
    int s=numBlocks;
    int kernel = whichKernel;

    // while(s > cpuFinalThreshold) {
    while(s > 1) {
        int threads = 0, blocks = 0;
        getNumBlocksAndThreads(kernel, s, maxBlocks, maxThreads, blocks, threads);
                
        reduce_sum<T>(s, threads, blocks, kernel, d_odata, d_odata);
                
        if (kernel < 3)
            s = (s + threads - 1) / threads;
        else
            s = (s + (threads*2-1)) / (threads*2);
    }

    // result = d_odata;
    cudaMalloc((void**) &result, sizeof(T));
    cudaMemcpy(result, d_odata, sizeof(T), cudaMemcpyDeviceToDevice);
    cudaFree(d_odata);
    return result;  //TUAN: result is a pointer, need to free in the calling function
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
// MIN REDUCTION
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
/* Find the minimum value on GPU
   Support Double Precision only
   Use kernel 6
*/   
double
min_double_2host_tesla(double* d_idata, int size, double* d_odata)
{
    //use 256 for FERMI
    // and 128 for Tesla
    int maxThreads = 128;
    int numThreads = 128;
    int maxBlocks = MAXBLOCK_TESLA;
    int numBlocks = NUMBLOCK_TESLA;

    cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(double), cudaMemcpyDeviceToDevice);
    
    // execute the kernel
    reduce6_min<double>(size, numThreads, numBlocks, d_idata, d_odata);
    //sum partial block sums on GPU
    int s=numBlocks;

    // while(s > cpuFinalThreshold) {
    while(s > 1) {
        int threads = 0, blocks = 0;
        getNumBlocksAndThreads_6(s, maxBlocks, maxThreads, blocks, threads);
                
        reduce6_min<double>(s, threads, blocks, d_odata, d_odata);
                
        s = (s + (threads*2-1)) / (threads*2);
    }

    double result;
    cudaMemcpy(&result, d_odata, sizeof(double), cudaMemcpyDeviceToHost);
    return result; 
}

double
max_double_2host_tesla(double* d_idata, int size, double* d_odata)
{
    //use 256 for FERMI
    // and 128 for Tesla
    int maxThreads = 128;
    int numThreads = 128;
    int maxBlocks = MAXBLOCK_TESLA;
    int numBlocks = NUMBLOCK_TESLA;
    // double* d_odata = NULL;
    // cudaMalloc((void**) &d_odata, numBlocks*sizeof(double));
    cudaMemcpy(d_odata, d_idata, numBlocks*sizeof(double), cudaMemcpyDeviceToDevice);
    

    // execute the kernel
    reduce6_max<double>(size, numThreads, numBlocks, d_idata, d_odata);
    //sum partial block sums on GPU
    int s=numBlocks;

    // while(s > cpuFinalThreshold) {
    while(s > 1) {
        int threads = 0, blocks = 0;
        getNumBlocksAndThreads_6(s, maxBlocks, maxThreads, blocks, threads);
                
        reduce6_max<double>(s, threads, blocks, d_odata, d_odata);
                
        s = (s + (threads*2-1)) / (threads*2);
    }

    double result;
    cudaMemcpy(&result, d_odata, sizeof(double), cudaMemcpyDeviceToHost);
    return result; 
}

////////////////////////////////////////////////////////////////////////////////
//! Compute sum reduction on CPU
//! We use Kahan summation for an accurate sum of large arrays.
//! http://en.wikipedia.org/wiki/Kahan_summation_algorithm
//! 
//! @param data       pointer to input data
//! @param size       number of input data elements
////////////////////////////////////////////////////////////////////////////////
template<class T>
T reduceCPU(T *data, int size)
{
    T sum = data[0];
    T c = (T)0.0;              
    for (int i = 1; i < size; i++)
        {
            T y = data[i] - c;  
            T t = sum + y;      
            c = (t - sum) - y;  
            sum = t;            
        }
    return sum;
}

extern "C" double reduceCPU_double(double *data, int size)
{
    return reduceCPU<double>(data, size);
}


// double sum_reduce_double_2host(double* p_arr, int n, int maxThreads=256, int whichKernel=6, int maxBlocks=64) {
// double * sumval;
// sumval = sum_reduce<double>(p_arr, n, maxThreads, whichKernel, maxBlocks);
// double res;
// cudaMemcpy(&res, sumval, sizeof(T), cudaMemcpyDeviceToHost);
// return res;
// }
