   dim3 threadsPerBlock(__BLOCK_DIM, __BLOCK_DIM);
 
   if (diff >= 0)
-    subtract_psf_kernel<<<numBlocks, threadsPerBlock>>>(res_p, psf_p + diff, stopx, stopy, diff, linsize, found_place_res.val * gain);
+    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, found_place_res.val * gain);
+    subtract_psf_kernel<<<numBlocks, threadsPerBlock>>>(res_p - diff, psf_p, stopx, stopy, diff, linsize, peak_x_gain);
 }
 
 // All pointers are device pointers

  , 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;
}


#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 peak_x_gain
  ) {
  const lli diff = peak_psf_pos - peak_res_pos;
  lldiv_t
      resxy = div(peak_res_pos, (lli)linsize)
    , psfxy = div(peak_psf_pos, (lli)linsize)
    ;
  const int
      stopx = linsize - abs (psfxy.rem - resxy.rem)
