   }
 }
 
+/* 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__ double s_red[BDIM];
+  __shared__ datafloat s_red[BDIM];
 
   const int id = blockIdx.x*blockDim.x + threadIdx.x;
 

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

  cudaEventCreate(&endEvent);

  cudaEventRecord(startEvent, 0);

  datafloat res;

  do {


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

