In this work we study the stellar-dynamical hardening of unequal mass massive black hole (MBH) binaries in the central regions of galactic nuclei. We present a comprehensive set of direct N-body simulations of the problem, varying both the total mass and the mass ratio of the MBH binary. Our initial model starts as an axisymmetric, rotating galactic nucleus, to describe the situation right after the galaxies have merged, but the black holes are still unbound to each other. We confirm that results presented in earlier works (Berczik et al. 2006; Khan et al. 2013; Wang et al. 2014) about the solution of the “last parsec problem” (sufficiently fast black hole coalescence for black hole growth in cosmological context) are robust for both for the case of unequal black hole masses and large particle numbers. The MBH binary hardening rate depends on the reduced mass ratio through a single parameter function, which quantitatively quite well agrees with standard 3 body scattering theory (see e.g., Hills 1983). Based on our results we conclude that MBH binaries at high redshifts are expected to merge with a factor of ~ 2 more efficiently, which is important to determine the possible overall gravitational wave signals. However, we have not yet fully covered all the possible parameter space, in particular with respect to the preceding of the galaxy mergers, which may lead to a wider variety of initial models, such as initially more oblate and / or even significantly triaxial galactic nuclei. Our N-body simulations were carried out on a new special supercomputers using the hardware acceleration with graphic processing units (GPUs).