  We present the first detailed comparison between million-body globular
cluster simulations computed with a H\'enon-type Monte Carlo code, CMC, and a
direct $N$-body code, NBODY6++GPU. Both simulations start from an identical
cluster model with $10^6$ particles, and include all of the relevant physics
needed to treat the system in a highly realistic way. With the two codes
"frozen" (no fine-tuning of any free parameters or internal algorithms of the
codes) we find excellent agreement in the overall evolution of the two models.
Furthermore, we find that in both models, large numbers of stellar-mass black
holes (> 1000) are retained for 12 Gyr. Thus, the very accurate direct $N$-body
approach confirms recent predictions that black holes can be retained in
present-day, old globular clusters. We find only minor disagreements between
the two models and attribute these to the small-$N$ dynamics driving the
evolution of the cluster core for which the Monte Carlo assumptions are less
ideal. Based on the overwhelming general agreement between the two models
computed using these vastly different techniques, we conclude that our Monte
Carlo approach, which is more approximate, but dramatically faster compared to
the direct $N$-body, is capable of producing a very accurate description of the
long-term evolution of massive globular clusters even when the clusters contain
large populations of stellar-mass black holes.
