 
 #define __BLOCK_DIM 16
 
+typedef long long int lli;
+
 inline
 __host__
 static void subtractPSF(
           double * res_p
   , const double * psf_p
+  , const lli peak_res_pos
+  , const lli peak_psf_pos
   , const int linsize
-  , const double gain
+  , const double peak_x_gain
   ) {
-  typedef long long int lli;
-  const lli diff = (lli)found_place_psf.pos - (lli)found_place_res.pos;
+  const lli diff = peak_psf_pos - peak_res_pos;
   lldiv_t
-      resxy = div(found_place_res.pos, (lli)linsize)
-    , psfxy = div(found_place_psf.pos, (lli)linsize)
+      resxy = div(peak_res_pos, (lli)linsize)
+    , psfxy = div(peak_psf_pos, (lli)linsize)
     ;
   const int
       stopx = linsize - abs (psfxy.rem - resxy.rem)

  , const double threshold
  ) {
  cudaError_t err;

  for (unsigned int i = 0; i < niters; ++i) {
    err = findPeaks(res_p, psf_p, linsize * linsize);
    if (err != cudaSuccess) return err;

    if (abs(found_place_res.val) < threshold) break;

    subtractPSF(res_p, psf_p, linsize, gain);
    mod_p[found_place_res.pos] += found_place_res.val * gain;
    // Wait for the PSF subtraction to finish
    err = cudaDeviceSynchronize();
  }

  return err;
}

  dim3 threadsPerBlock(__BLOCK_DIM, __BLOCK_DIM);

  if (diff >= 0)
    subtract_psf_kernel<<<numBlocks, threadsPerBlock>>>(res_p, psf_p + diff, stopx, stopy, diff, linsize, peak_x_gain);
  else
    subtract_psf_kernel<<<numBlocks, threadsPerBlock>>>(res_p - diff, psf_p, stopx, stopy, diff, linsize, peak_x_gain);
}

// All pointers are device pointers
