This paper presents an efficient approach for stabilizing solution and accelerating convergence of a harmonic balance equation system for an efficient analysis of turbomachinery unsteady flows due to flutter and blade row interaction. The proposed approach combines the Runge–Kutta method with the lower upper symmetric Gauss Seidel (LU-SGS) method and the block Jacobi method. The LU-SGS method, different from its original application as an implicit time marching scheme, is used as an implicit residual smoother with under-relaxation, allowing big Courant–Friedrichs–Lewy (CFL) numbers (in the order of hundreds), leading to significant convergence speedup. The block Jacobi method is introduced to implicitly integrate the time spectral source term of a harmonic balance equation system, in order to reduce the complexity of the direct implicit time integration by the LU-SGS method. The implicit treatment of the time spectral source term thus greatly augments the stability region of a harmonic balance equation system in the case of grid-reduced frequency well above ten. Validation of the harmonic balance flow solver was carried out using linear cascade test data. Flutter analysis of a transonic rotor and blade row interaction analyses for a transonic compressor stage were presented to demonstrate the stabilization and acceleration effect by the combination of the LU-SGS and the block Jacobi methods. The influence of the number of Jacobi iterations on solution stabilization is also investigated, showing that two Jacobi iterations are sufficient for stability purpose, which is much more efficient than existing methods of its kind in the open literature.