A fast numerical algorithm for nonlinear finite sound beams is studied.Multiple couplings in nonlinear finite sound beams are found,including full coupling between the harmonics,recursive coupling between field points along the axial direction,and local coupling along the radial direction.Calculation efficiency is improved by a spatial-domain decomposition of acoustic field along radial direction.The method is achieved through a multi-threading parallel computation.Calculation results for nonlinear focused Gaussian beams show that as problem size and the number of threads are reasonablly matched,this algorithm can significantly improve calculation speed while ensuring accuracy.Calculation results are compared with theoretical analysis.