This paper presents a fast GPU-based solution to compensate positron range effects in heterogeneous media for iterative PET reconstruction. We assume a factorized approach, where projections are decomposed to phases according to the main physical effects. Positron range is the first effect in this chain, which causes a spatially varying blurring according to local material properties. In high-resolution small animal PET systems, the average free path length of positrons may be many times longer than the linear size of voxels. This means that positron range significantly compromises the reconstruction quality if it is not compensated, and also that the material dependent blurring should have a very large support so its voxel space calculation would take prohibitively long. Frequency domain filtering does not have such computational complexity problems, but its direct form is ruled out by the fact that we need a spatially variant filtering in heterogeneous media. To handle heterogeneous media, we execute Fast Fourier Transforms for each material type and for appropriately modulated tracer densities and merge these partial results into a density that describes the composed, heterogeneous medium. Fast Fourier Transform requires the filter kernels on the same resolution as the tracer density is defined, so we also present efficient methods for re-sampling the probability densities of positron range for different resolutions and basis functions. The algorithm is implemented on the GPU, built into the TeraTomo™ system, and requires just a few seconds on high resolution voxel arrays.