  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 or transneptunian objects in the 5:2
and 3:1 mean-motion resonances, demonstrating that long-term perturbations
remove the initially stable Neptune-crossing members.
