   cudaEventCreate(&endEvent);
 
   cudaEventRecord(startEvent, 0);
+
+  datafloat res;
+
   do {
+
+
     // Call jacobi [iterationChunk] times before calculating residual
     for(int i = 0; i < iterationChunk/2; ++i){
 

    red[blockIdx.x] = s_red[threadIdx.x];
}

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

  if(argc != 3){
    printf("Usage: ./main N tol \n");
    return 0;
  }

  const int N     = atoi(argv[1]);
  const datafloat tol = atof(argv[2]);
  datafloat res = 0;

  const int iterationChunk = 100; // needs to be multiple of 2
  int iterationsTaken = 0;

  // Setup jacobi kernel block-grid sizes
  dim3 jBlock(16,16);
  dim3 jGrid((N + jBlock.x - 1)/jBlock.x, (N + jBlock.y - 1)/jBlock.y);

  // Setup reduceResidual kernel block-grid sizes
  dim3 rBlock(BDIM);
  dim3 rGrid((N*N + rBlock.x - 1)/rBlock.x);

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

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

  datafloat delta = 2./(N+1);
  datafloat normRHS = 0;
  for(int j = 0; j < N+2; ++j){
    for(int i = 0; i < N+2; ++i){
      datafloat x = -1 + delta*i;
      datafloat 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);
    }
  }
  normRHS = sqrt(normRHS);
  
  // Device Arrays
  datafloat *c_u, *c_u2, *c_rhs, *c_res;



#define datafloat double

#define BX 16
#define BY 16

/* 

   Poisson problem:               diff(u, x, 2) + diff(u, y, 2) = f

   Coordinate transform:          x -> -1 + delta*i, 
                                  y -> -1 + delta*j

   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))

   To run with a 402x402 grid until the solution changes less than 1e-7 per iteration (in l2): ./main 400 1e-7  

*/

