 
   if(x<width && y<height)
   {
-    // Load p
     __shared__ float2 p[BSX][BSY];
     p[threadIdx.x][threadIdx.y] = tex2D(rof_tex_p, x+0.5f, y+0.5f);
 
-    // Load u_
+    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;
+
     __shared__ float u_[BSX][BSY];
-    u_[threadIdx.x][threadIdx.y] = tex2D(rof_tex_u_, x+0.5f, y+0.5f);
+
+    float f = tex2D(rof_tex_f, x+0.5f, y+0.5f);
+    float u = tex2D(rof_tex_u, x+0.5f, y+0.5f);
+    u_[threadIdx.x][threadIdx.y] = u;
+
+    syncthreads();
+
+    // Primal update
+    if ( (threadIdx.x > 0) && (threadIdx.y > 0) )
+    {
+      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);
+      u_[threadIdx.x][threadIdx.y] = u + theta*(u-u_[threadIdx.x][threadIdx.y]);
+    }
 
     syncthreads();
 


  bindTexture(rof_tex_f, device_f);
  bindTexture(rof_tex_u, device_u);
  bindTexture(rof_tex_u_, device_u_);
  bindTexture(rof_tex_p, device_p);

  float L = sqrt(8.0f);

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

////////////////////////////////////////////////////////////////////////////////
__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)
{
