       p[threadIdx.x][threadIdx.y] = p[threadIdx.x][threadIdx.y] + sigma*grad_u;
       p[threadIdx.x][threadIdx.y] = p[threadIdx.x][threadIdx.y]/max(1.0f, length(p[threadIdx.x][threadIdx.y]));
 
-      if (x<=0)
-        p[threadIdx.x][threadIdx.y].x = 0.0f;
-      if (y<=0)
-        p[threadIdx.x][threadIdx.y].y = 0.0f;
-
-      if (x >= width-1)
-        p[threadIdx.x][threadIdx.y].x = 0.0f;
-      if (y >= height-1)
-        p[threadIdx.x][threadIdx.y].y = 0.0f;
-
-      // Write back
+      // Write Back
       if ( (threadIdx.x > 0) && (threadIdx.y > 0) )
+      {
         device_p[c2] = p[threadIdx.x][threadIdx.y];
-    }
-
-        syncthreads();
-
-    // texture fetches
-    float f = tex2D(rof_tex_f, x+0.5f, y+0.5f);
-    float u = tex2D(rof_tex_u, x+0.5f, y+0.5f);
-
-    // Primal update
-    if ( (threadIdx.x > 0) && (threadIdx.y > 0) )
-    {
-      u_[threadIdx.x][threadIdx.y] = u;
-
-      float divergence = p[threadIdx.x][threadIdx.y].x - p[threadIdx.x-1][threadIdx.y].x +
-                         p[threadIdx.x][threadIdx.y].y - p[threadIdx.x][threadIdx.y-1].y;
-
-      // update primal variable
-      u = (u + tau*(divergence + lambda*f))/(1.0f+tau*lambda);
-
-      if ( (threadIdx.x < BSX-1) && (threadIdx.y < BSY-1) )
-      {
-        // Write back
         device_u[c1] = u;
-        device_u_[c1] = u + theta*(u-u_[threadIdx.x][threadIdx.y]);
       }
     }
   }

    if (threadIdx.y==BSY-1)
      u_[threadIdx.x][threadIdx.y+1] = tex2D(rof_tex_u_, x+0.5f, y+1.5f);

//    syncthreads();

    float2 grad_u = make_float2(u_[threadIdx.x+1][threadIdx.y] - u_[threadIdx.x][threadIdx.y],
                                u_[threadIdx.x][threadIdx.y+1] - u_[threadIdx.x][threadIdx.y] );

////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
__global__ void update_shared_single(float* device_u, float2* device_p,
                                     float tau, float theta, float sigma, float lambda,
                                     int width, int height, int stride1, int stride2)
{
