  The numerical simulations of massive collisional stellar systems, such as
globular clusters (GCs), are very time-consuming. Until now, only a few
realistic million-body simulations of GCs with a small fraction of binaries
(5%) have been performed by using the NBODY6++GPU code. Such models took half a
year computational time on a GPU based super-computer. In this work, we develop
a new N-body code, PeTar, by combining the methods of Barnes-Hut tree, Hermite
integrator and slow-down algorithmic regularization (SDAR). The code can
accurately handle an arbitrary fraction of multiple systems (e.g. binaries,
triples) while keeping a high performance by using the hybrid parallelization
methods with MPI, OpenMP, SIMD instructions and GPU. A few benchmarks indicate
that PeTar and NBODY6++GPU have a very good agreement on the long-term
evolution of the global structure, binary orbits and escapers. On a highly
configured GPU desktop computer, the performance of a million-body simulation
with all stars in binaries by using PeTar is 11 times faster than that of
NBODY6++GPU. Moreover, on the Cray XC50 supercomputer, PeTar well scales when
number of cores increase. The ten million-body problem, which covers the region
of ultra compact dwarfs and nuclearstar clusters, becomes possible to be
solved.
