   }
 
   const int N     = atoi(argv[1]);
-  const float tol = atof(argv[2]);
-  float res = 0;
+  const datafloat tol = atof(argv[2]);
+  datafloat res = 0;
 
   const int iterationChunk = 100; // needs to be multiple of 2
   int iterationsTaken = 0;

  dim3 rGrid((N*N + rBlock.x - 1)/rBlock.x);

  // Host Arrays
  float *h_u   = (float*) calloc((N+2)*(N+2), sizeof(float));
  float *h_u2  = (float*) calloc((N+2)*(N+2), sizeof(float));
  float *h_rhs = (float*) calloc((N+2)*(N+2), sizeof(float));

  float *h_res = (float*) calloc(rGrid.x, sizeof(float));

  float delta = 2./(N+1);
  float normRHS = 0;
  for(int j = 0; j < N+2; ++j){
    for(int i = 0; i < N+2; ++i){
      float x = -1 + delta*i;
      float y = -1 + delta*j;

      h_rhs[i + (N+2)*j] = delta*delta*(2.*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y));
      normRHS += pow(h_rhs[i+(N+2)*j],2);


   2nd order finite difference:   4*u(j,i) - u(j-1,i) - u(j+1,i) - u(j,i-1) - u(j,i+1) = -delta*delta*f(j,i) 

   define: rhs(j,i) = -delta*delta*f(j,i)

   Jacobi iteration: newu(j,i) = 0.25*(u(j-1,i) + u(j+1,i) + u(j,i-1) + u(j,i+1) + rhs(j,i))

*/

__global__ void jacobi(const int N,
                       const datafloat *rhs,
                       const datafloat *u,
                       datafloat *newu){

  // Get thread indices
  const int i = blockIdx.x*blockDim.x + threadIdx.x;
  const int j = blockIdx.y*blockDim.y + threadIdx.y;

