Modern computers are generally equipped with multi-core central processing units (CPUs). There has been great interest to introduce a parallelized harmonic balance method to make full use of computing resources and improve the efficiency in dealing with nonlinear dynamic analysis of high-dimensional spatially discretized models of continuous systems. In this work, an optimized efficient Galerkin averaging-incremental harmonic balance method for solving high-dimensional models of continuous systems based on parallel computing is introduced. Optimized parallel implementation based on tensor contraction is introduced in time-domain series calculations and the quasi-Newton method is used in the iteration procedure, which greatly accelerate computational speeds of both serial and parallel implementations. Especially, the parallel implementation achieves high parallel efficiency when multiple CPU cores are used. Due to its high computational efficiency and good robustness, the proposed method has the potential to be used as a powerful universal solver and analyzer for general types of continuous systems.