  We review the recent optimizations of gravitational $N$-body kernels for
running them on graphics processing units (GPUs), on single hosts and massive
parallel platforms. For each of the two main $N$-body techniques, direct
summation and tree-codes, we discuss the optimization strategy, which is
different for each algorithm. Because both the accuracy as well as the
performance characteristics differ, hybridizing the two algorithms is essential
when simulating a large $N$-body system with high-density structures containing
few particles, and with low-density structures containing many particles. We
demonstrate how this can be realized by splitting the underlying Hamiltonian,
and we subsequently demonstrate the efficiency and accuracy of the hybrid code
by simulating a group of 11 merging galaxies with massive black holes in the
nuclei.
