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

    else
      theta = 1.0f;

    update_shared_single<<<dimGrid, dimBlock>>>(device_u->data(), device_u_->data(), device_p->data(),
                                                tau, theta, sigma, lambda, width, height,
                                                device_u->stride(), device_p->stride());
    sigma /= theta;


  if(x<width && y<height)
  {
    __shared__ float2 p[BSX][BSY];
    p[threadIdx.x][threadIdx.y] = tex2D(rof_tex_p, x+0.5f, y+0.5f);

    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];

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

