#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// Thread block size
#define BLOCK_SIZE 16
#define TILE_SIZE  16

#define WA 512   // Matrix A width
#define HA 512   // Matrix A height

__global__ void
poisson( float* A, int wA)
{

   // 1. 2D Thread ID
   int iid=0;
   //float value = 0.000000001;
   //float delta_u=-1.0;
   do
   {
   int tx = blockIdx.x * TILE_SIZE + threadIdx.x;
   int ty = blockIdx.y * TILE_SIZE + threadIdx.y;
   //int tx = threadIdx.x;
   //int ty = threadIdx.y;

    if (ty>1 && ty<WA-1 && tx>1 && tx<HA-1)
        {
        //float tmp_delta_u=A[ty * wA + tx];

        float elementA = A[(ty*wA) + (tx+1)];
        elementA += A[(ty*wA) + (tx-1)];
        elementA += A[((ty+1)*wA) + tx];
        elementA += A[((ty-1)*wA) + tx];
        elementA *= 0.25;

        //tmp_delta_u-=elementA ;
        A[ty * wA + tx] = elementA ;

        //if (fabs(tmp_delta_u)>delta_u)
        //    delta_u=fabs(tmp_delta_u);
        };
    }
    //while (value<delta_u);
    while (iid++<20000);

}

void randomInit(float* data, int size)
{
    for (int i = 0; i < size; ++i)
        {
        data[i] = i%13;
        };
}

/////////////////////////////////////////////////////////
// Program main
/////////////////////////////////////////////////////////

int
main(int argc, char** argv)
{

    srand(10);

    unsigned int size_A = WA * HA;
    unsigned int mem_size_A = sizeof(float) * size_A;
    float* h_A = (float*) malloc(mem_size_A);

    randomInit(h_A, size_A);

    float* d_A;
    cudaMalloc((void**) &d_A, mem_size_A);

    cudaMemcpy(d_A, h_A, mem_size_A,
               cudaMemcpyHostToDevice);


    dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
    dim3 dimGrid( WA/BLOCK_SIZE,HA/BLOCK_SIZE);

    poisson<<<dimGrid, dimBlock>>>(d_A,WA);
    cudaMemcpy(h_A, d_A, mem_size_A,
               cudaMemcpyDeviceToHost);

    printf("\n\nMatrix C (Results)\n");
    for(int i = 0; i < size_A; i++)
    {
       printf("%.1f ", h_A[i]);
       if(((i + 1) % WA) == 0)
          printf("\n");
    }
    printf("\n");

    // 7. clean up memory
    free(h_A);
    cudaFree(d_A);

}