This paper focuses on the parallel implementation of a direct N -body method (particleparticle algorithm) and the application of multiple GPUs for galactic dynamics simulations. Application of a hybrid OpenMP-CUDA technology is considered for models with a number of particles N ∼ 10 5 ÷ 10 7 . By means of N -body simulations of gravitationally unstable stellar galactic we have investigated the algorithms parallelization efficiency for various Nvidia Tesla graphics processors (K20, K40, K80). Particular attention was paid to the parallel performance of simulations and accuracy of the numerical solution by comparing single and double floating-point precisions (SP and DP). We showed that the doubleprecision simulations are slower by a factor of 1.7 than the single-precision runs performed on Nvidia Tesla K-Series processors. We also claim that application of the single-precision operations leads to incorrect result in the evolution of the non-axisymmetric gravitating N -body systems. In particular, it leads to significant quantitative and even qualitative distortions in the galactic disk evolution. For instance, after 10 4 integration time steps for the single-precision numbers the total energy, momentum, and angular momentum of a system with N = 2 20 conserve with accuracy of 10 −3 , 10 −2 and 10 −3 respectively, in comparison to the double-precision simulations these values are 10 −5 , 10 −15 and 10 −13 , respectively. Our estimations evidence in favour of usage of the second-order accuracy schemes with double-precision numbers since it is more efficient than in the fourth-order schemes with single-precision numbers.