Abstract This article puts forward a computationally efficient block matrix based precise integration algorithm for solving vibration response subjected to time-variable excitation and nonlinearity, especially for nonhomogeneous dynamic response solution with large-scale degrees-of-freedom. In detail, the nonlinear parts and time-varying inputs of a dynamic system are separated from the original dynamic equations and then simulated within a computing time-step by employing a piecewise interpolation function. A novel closed-form iteration formula is presented in conjunction with the block matrix strategy and modified increment-dimensional precise integration technique. Interestingly, the presented approach is essentially a high-accuracy and parallel algorithm, which exhibits a high prediction accuracy without the limitation of matrix inversion, higher-order derivative, periodicity requirement nor cycle oscillation and instability of high-order interpolation. The feasibility and advantage of the proposed method are verified with two numerical examples.