If your discretization of the reaction part of the problem is explicit, you can still use Prof. Bangerth's suggestion. You can quickly compute $K*u$ using the fast Fourier transform on both $K$ and $u$, multiplying the transforms and then transforming back. Then you can easily compute $u\cdot K*u$ at any grid point. However, you may run into stability issues with an explicit time discretization and a reasonably large time-step; reaction-diffusion equations are often stiff. We can probably give you more suggestions if you write out what the kernel is.
Since the kernel is rapidly decreasing far from the origin, you could also try to approximate the convolution $K * u$ using something like the fast multipole method. This idea was originally used to turn the $O(N^2)$ computation of gravitational $N$-body potentials into an $O(N\log N)$ computation of an approximate potential, but it's been extended to more general problems. FMM often shows up on lists of the top 10 greatest algorithms, so it's also worth knowing about on principle.
The fast multipole method and more general hierarchical matrix methods are often used for the numerical solution of integral equations.