We present a GPU-accelerated numerical integrator specifically optimized for stability calculations of small bodies in planetary systems. Specifically, the integrator is designed for cases when large numbers of test particles (tens or hundreds of thousands) need to be followed for long durations (millions of orbits) to assess the orbital stability of their initially “close-encounter free” orbits. The GLISSE (Gpu Long-term Integrator for Solar System Evolution) code implements several optimizations to achieve a roughly factor of 100 speed increase over running the same code on a CPU. We explain how various hardware speed bottlenecks can be avoided by the careful code design, although some of the choices restrict the usage to specific types of application. As a first application, we study the long-term stability of small bodies initially on orbits between Uranus and Neptune. We map out in detail the small portion of the phase space in which small bodies can survive for 4.5 billion years of evolution; the ability to integrate large numbers of particles allow us to identify for the first time how instability-inducing mean-motion resonance overlaps sharply define the stable regions. As a second application, we map the boundaries of 4 Gyr stability for transneptunian objects in the 5:2 and 3:1 mean-motion resonances, demonstrating that long-term perturbations remove the initially stable Neptune-crossing members.