   // Copy final solution from device array to host
   cudaMemcpy(h_u, c_u, (N+2)*(N+2)*sizeof(datafloat), cudaMemcpyDeviceToHost);
 
-  // Compute maximum error in finite difference solution
+  const datafloat avgTimePerIteration = timeTaken/((datafloat) iterationsTaken);
+
+  printf("Residual                   : %7.9e\n"     , res);
+  printf("Iterations                 : %d\n"        , iterationsTaken);
+  printf("Average time per iteration : %3.5e ms\n"  , avgTimePerIteration);
+  printf("Bandwidth                  : %3.5e GB/s\n", (1.0e-6)*(6*N*N*sizeof(datafloat))/avgTimePerIteration);
+
+  // free device arrays
+  cudaFree(c_u);
+  cudaFree(c_u2);
+  cudaFree(c_res);
+  cudaFree(c_rhs);
+
+}
+
+int main(int argc, char** argv){
+
+  // parse command line arguements
+  if(argc != 3){
+    printf("Usage: ./main N tol \n");
+    return 0;
+  }
+
+  // Number of internal domain nodes in each direction
+  const int N     = atoi(argv[1]);
+
+  // Termination criterion || unew - u ||_2 < tol 
+  const datafloat tol = atof(argv[2]);
+
+  // 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((N+2)*(N+2), sizeof(datafloat));
+
+  // FD node spacing
+  datafloat delta = 2./(N+1);
+
+  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;
+
+      // solution is u = sin(pi*x)*sin(pi*y) so the rhs is: 
+      h_rhs[i + (N+2)*j] = delta*delta*(2.*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y));
+    }
+  }
+
+  // Solve discrete Laplacian
+  solve(N, tol, h_rhs, h_res, h_u, h_u2);
+
+  // Compute maximum error in finite difference solution and output solution
   FILE *fp = fopen("result.dat", "w");
   datafloat maxError = 0;
   for(int j = 0; j < N+2; ++j){

  cudaEventCreate(&endEvent);

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


  }
}

/* CUDA kernel using shared memory */
__global__ void jacobiShared(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;

  __shared__ float s_u[BY+2][BX+2];

  // load block of data in to shared memory
  for(int n=threadIdx.y; n<BY+2; n+=BY){
    for(int m=threadIdx.x; m<BX+2; m+=BX){

      const int i = blockIdx.x*blockDim.x + m;
      const int j = blockIdx.y*blockDim.y + n;
      
      datafloat val = 0.;

      if(i < N+2 &&  j < N+2)
        val = u[j*(N+2) + i];
      
      s_u[n][m] = val;
    }
  }

  // barrier until all the shared memory entries are loaded
  __syncthreads();

  if((i < N) && (j < N)){

    // Get padded grid ID
    const int pid = (j + 1)*(N + 2) + (i + 1);

    datafloat invD = 0.25;

    newu[pid] = invD*(rhs[pid]
		      + s_u[threadIdx.y+0][threadIdx.x+1]
		      + s_u[threadIdx.y+2][threadIdx.x+1]
		      + s_u[threadIdx.y+1][threadIdx.x+0]
		      + s_u[threadIdx.y+1][threadIdx.x+2]
		      );
  }
}


__global__ void partialReduceResidual(const int entries,
                               datafloat *u,
                               datafloat *newu,
                               datafloat *red){
  __shared__ datafloat s_red[BDIM];

  const int id = blockIdx.x*blockDim.x + threadIdx.x;

