// Includes.
#include <shrUtils.h>

#define NRAWR 500
#define NRAWC 3000
#define NGPSR 2
#define NGPSC 3000
#define NDIR 201
#define NDIC 2021
#define NIMGR 201
#define NIMGC 2021
#define BLOCK_SIZE 24

typedef struct { 
	int width; 
	int height; 
	int stride; 
	float* elements; 
} Matrix;

typedef struct { 
	int width; 
	int height; 
	int stride; 
	float2* elements; 
} Matrix2;


/* Headers */
void creatematrix (Matrix*,int,int);
void creatematrix2(Matrix2*,int,int);
__global__ void kernel (Matrix, Matrix, Matrix2, Matrix);

Matrix RAW, IMG, GPS;
Matrix2 DI;

int main()
{ 
	int i, j;
	Matrix dr, di, dg;
	Matrix2 dd;

	creatematrix(&IMG,NDIR,NDIC);
					
	creatematrix(&RAW,NRAWR,NRAWC);
	for(j=0;j<NRAWR;j++)
	{
		for(i=0;i<NRAWC;i++)
		{
			RAW.elements[j*RAW.width+i]=j+1;
		}
	}

	i=NGPSC/(BLOCK_SIZE*BLOCK_SIZE);
	i++;
	i=i*(BLOCK_SIZE*BLOCK_SIZE);

	creatematrix(&GPS,NGPSR,i);
	for(j=0;j<GPS.width;j++)
	{
			GPS.elements[j]=0.3*(j+1);
			GPS.elements[GPS.width+j]=0;
	}

	creatematrix2(&DI,NDIR,NDIC);
	for(j=0;j<NDIR;j++)
	{
		for(i=0;i<NDIC;i++)
		{
			DI.elements[j*DI.width+i].x=-10+i*0.5;
			DI.elements[j*DI.width+i].y=200+j;
		}
	}

	//Allocate and copy in device memory.

	cudaFuncSetCacheConfig(kernel, cudaFuncCachePreferL1);

	di.height=IMG.height;
	di.width=di.stride=IMG.width;
	i = di.height*di.width*sizeof(float);
	cudaMalloc((void**) &di.elements, i);

	dr.height=RAW.height;
	dr.width=dr.stride=RAW.width;
	i = dr.height*dr.width*sizeof(float);
	cudaMalloc((void**) &dr.elements, i);
	cudaMemcpy(dr.elements, RAW.elements, i,
                              cudaMemcpyHostToDevice);

	dg.height=GPS.height;
	dg.width=dg.stride=GPS.width;
	i = dg.height*dg.width*sizeof(float);
	cudaMalloc((void**) &dg.elements, i);
	cudaMemcpy(dg.elements, GPS.elements, i,
                              cudaMemcpyHostToDevice) ;

	dd.height=DI.height;
	dd.width=dd.stride=DI.width;
	i = dd.height*dd.width*sizeof(float2);
	cudaMalloc((void**) &dd.elements, i);
	cudaMemcpy(dd.elements, DI.elements, i,
                              cudaMemcpyHostToDevice) ;

	int divisorx=BLOCK_SIZE, divisory=BLOCK_SIZE, 
		Numx = DI.width, Numy =  DI.height;
      
     while(divisorx<Numx)
     {
		 while( Numx % divisorx != 0)
			 divisorx=divisorx+1;
		 if(Numx==divisorx)
		 {
			 printf("t %ld es PRIMO n",Numx);    
			return 0; 
		}
      break;
      }

	 while(divisory<Numy)
     {
		 while( Numy % divisory != 0)
			 divisory=divisory+1;
		 if(Numy==divisory)
		 {
			 printf("t %ld es PRIMO n",Numy); 
			 return 0;
		 }
      break;
      }

	dim3 threads (BLOCK_SIZE, BLOCK_SIZE);

	dim3 grid (dd.width / divisorx, dd.height / divisory);

	kernel<<< grid, threads >>>(dr, dg, dd, di);

	i=di.height*di.width*sizeof(float);

	cudaMemcpy(IMG.elements, di.elements, i,
				cudaMemcpyDeviceToHost);

	cudaFree(dr.elements);
	cudaFree(dd.elements);
	cudaFree(dg.elements);
	cudaFree(di.elements);

	return 0;
		}

void creatematrix(Matrix* mat,int row,int column)
{
	int size;
	mat->height=row;
	mat->width=mat->stride=column;
	size = mat->height*mat->width*sizeof(float);
	mat->elements = (float*)malloc(size);
}

void creatematrix2(Matrix2* mat,int row,int column)
{
	int size;
	mat->height=row;
	mat->width=mat->stride=column;
	size = mat->height*mat->width*sizeof(float2);
	mat->elements = (float2*)malloc(size);
}

__global__ void kernel (Matrix dr, Matrix dg, Matrix2 dd, Matrix di)
{
	float temp, minx, sigma=998, resolution=2;
	short i, j, k, l, m, n, a, basis;

	unsigned short jmpcol =(dd.width/gridDim.x) * (1+blockIdx.x);
	unsigned short indy= blockIdx.x * (dd.width/gridDim.x) + threadIdx.x;

	unsigned short jmprow = (dd.height/gridDim.y) * (1+blockIdx.y);	
	unsigned short indx= blockIdx.y * (dd.height/gridDim.y) + threadIdx.y;

	unsigned short thnum = (BLOCK_SIZE*threadIdx.y)+threadIdx.x;

	float2 var;

	__shared__ float sdg[2][BLOCK_SIZE*BLOCK_SIZE];
	__shared__ int widthg;
	__shared__ int widthd;
	__shared__ int widthr;

	widthg=dg.width;
	widthd=dd.width;
	widthr=dr.width;

	a=BLOCK_SIZE*BLOCK_SIZE;
	n= widthg/a;

	for(l=0;l<n;l++)
	{
		sdg[0][thnum]=dg.elements[(a*l)+thnum];
		sdg[1][thnum]=dg.elements[(a*l)+thnum+widthg];
		__syncthreads();

		m = NGPSC - (a*l);

		if(m>a)
			m=a;

		for(i=indx;i<jmprow;i+=BLOCK_SIZE)
		{
			for(j=indy;j<jmpcol;j+=BLOCK_SIZE)
			{
 				var=dd.elements[i*widthd+j];
				for (k=0;k<m;k++)
				{
					minx = sdg[0][k]-var.x;
					temp = minx*minx;
					minx = sdg[1][k]-var.y;
					temp += minx*minx;
					temp = sqrtf (temp);
					if(temp<sigma)
					{
						basis=temp/resolution;
						di.elements[i*widthd+j]+=(dr.elements[(basis*widthr)+k]+dr.elements[(basis+1)*widthr+k])/2;
					}
				}
			}
		}
		__syncthreads();
	}
}