   s_red[threadIdx.x] = 0;
 
   if(id < entries){
-    const float diff = u[id] - newu[id];
+    const datafloat diff = u[id] - newu[id];
     s_red[threadIdx.x] = diff*diff;
   }
 


#define BDIM 256

/* 

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


   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;

