   }
   fclose(fp);
 
-  const datafloat avgTimePerIteration = timeTaken/((datafloat) iterationsTaken);
-
-  printf("Top right current          : %7.9e\n"     , h_u[N*(N+2) + N]);
-  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);
   printf("Maximum absolute error     : %7.9e\n"     ,  maxError);
 
   // Free all the mess
-  cudaFree(c_u);
-  cudaFree(c_u2);
-  cudaFree(c_res);
-
   free(h_u);
   free(h_u2);
   free(h_res);
-}
+  free(h_rhs);
+}

  }
}

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

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


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

void solve(int N, datafloat tol, datafloat *h_rhs, datafloat *h_res, datafloat *h_u, datafloat *h_u2 ){

  const int iterationChunk = 100; // needs to be multiple of 2
  int iterationsTaken = 0;
  
  // Setup jacobi kernel block-grid sizes
  dim3 jBlock(BX,BY);
  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);

  // Device Arrays
  datafloat *c_u, *c_u2, *c_rhs, *c_res;

