/*
 * 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.
 * 
 */


// nvcc -arch=sm_13 -o norm_const norm_const.cu

#include "../../cuda_by_example/common/book.h"

#define imin(a,b) (a<b?a:b)

const int B_size = 1000000;
const int A_size = 15131;
//const int B_size = 100;
//const int A_size = 1000;

const int threadsPerBlock = 1024;
const int blocksPerGrid =
            imin( 2048, (B_size+threadsPerBlock-1) / threadsPerBlock );


__global__ void normc( unsigned long int *a, unsigned long int *b, int *c ) {

    __shared__ float cache[threadsPerBlock];
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    int cacheIndex = threadIdx.x;

    int	tmp_num_in_A = 0;
    while (tid < B_size) {
        
	   for(int j=0; j<A_size; j++)
	   {
           if( (b[tid*2 + 0] == a[j*2 + 0]) && (b[tid*2 + 1] == a[j*2 + 1]) ) 
		 {
              tmp_num_in_A++;
		    break;
		 }
        }

        tid += blockDim.x * gridDim.x;
    }
    
    // set the cache values
    cache[cacheIndex] = tmp_num_in_A;
    
    // synchronize threads in this block
    __syncthreads();

    // for reductions, threadsPerBlock must be a power of 2
    // because of the following code
    int i = blockDim.x/2;
    while (i != 0) {
        if (cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    }

    if (cacheIndex == 0)
        c[blockIdx.x] = cache[0];
}


int main( void ) {
    unsigned long int   *a, *b; 
    int	c, *partial_c;
    unsigned long int   *dev_a, *dev_b; 
    int	*dev_partial_c;

    cudaEvent_t     start, stop;

    // allocate memory on the cpu side
    a = (unsigned long int*)malloc( A_size*2*sizeof(unsigned long int) );
    b = (unsigned long int*)malloc( B_size*2*sizeof(unsigned long int) );
    partial_c = (int*)malloc( blocksPerGrid*sizeof(int) );


    FILE	*Ain, *Bin;
    Ain = fopen("A.txt", "r");
    Bin = fopen("B.txt", "r");

    unsigned long int	tmp1, tmp2;

    // fill in the host memory with data
    for (int i=0; i<B_size; i++) {
        if(i < A_size) {
           fscanf(Ain, "%lu %lu\n", &tmp1, &tmp2);
	      a[i * 2 + 0] = tmp1;
	      a[i * 2 + 1] = tmp2;
	   }

	   fscanf(Bin, "%lu %lu\n", &tmp1, &tmp2);
	   b[i * 2 + 0] = tmp1;
	   b[i * 2 + 1] = tmp2;
        //b[i * 2 + 0] = 412316861472;
	   //b[i * 2 + 1] = 2199040041218;
    }

    fclose(Ain); fclose(Bin);

    HANDLE_ERROR( cudaEventCreate( &start ) );
    HANDLE_ERROR( cudaEventCreate( &stop ) );
    HANDLE_ERROR( cudaEventRecord( start, 0 ) );

    // allocate the memory on the GPU
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a,
                              A_size*2*sizeof(unsigned long int) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_b,
                              B_size*2*sizeof(unsigned long int) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c,
                              blocksPerGrid*sizeof(int) ) );


    // copy the arrays 'a' and 'b' to the GPU
    HANDLE_ERROR( cudaMemcpy( dev_a, a, A_size*2*sizeof(unsigned long int),
                              cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy( dev_b, b, B_size*2*sizeof(unsigned long int),
                              cudaMemcpyHostToDevice ) ); 

    normc<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b,
                                              dev_partial_c );

    // copy the array 'c' back from the GPU to the CPU
    HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c,
                              blocksPerGrid*sizeof(int),
                              cudaMemcpyDeviceToHost ) );

    // finish up on the CPU side
    c = 0;
    for (int i=0; i<blocksPerGrid; i++) {
        c += partial_c[i];
    }

    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    float   elapsedTime;
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
    								start, stop ) );
    printf( "Time elapsed:  %3.1f ms\n", elapsedTime );

    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );

    FILE *fout;
    fout = fopen("frequency.txt", "w");
    fprintf(fout, "%d\t%lf\n", c, (double)c / (double)B_size);
    fclose(fout);

    printf("Frequency = %d / %d = %f\n\n", c, B_size, (double)c / B_size);


    // free memory on the gpu side
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_b ) );
    HANDLE_ERROR( cudaFree( dev_partial_c ) );

    // free memory on the cpu side
    free( a );
    free( b );
    free( partial_c );
}