The modified BGK equation adapted to various flow regimes can be presented by the aid of the basic characteristics on molecular movement and collision approaching to equilibrium. The discrete velocity ordinate method is developed and applied to the velocity distribution function to remove its continuous dependency on the velocity space, and then the velocity distribution function equation is cast into hyperbolic conservation law form with nonlinear source terms. Based on the unsteady time-splitting method and the non-oscillatory, containing no free parameters, and dissipative (NND) scheme, the gas kinetic finite difference second-order scheme is constructed for the computation of the discrete velocity distribution functions. The mathematical model on the interaction of molecules with solid surface is studied and used in the numerical method. Four types of numerical quadrature rules, such as the modified Gauss-Hermite formula, the composite Newton-Cotes integration method, the Gauss-Legendre numerical quadrature rule, and the Golden Section number-theoretic integral method, are developed and applied to the discretized velocity space to evaluate the macroscopic flow parameters at each point in the physical space. As a result, a unified simplified gas kinetic algorithm is established for the flows from rarefied transition to continuum regime. Based on analyzing the inner parallel degree of the unified algorithm, the parallel strategy adapted to the gas kinetic numerical algorithm is studied, and then the HPF parallel processing software for the unified algorithm is developed. To test the present method, the one-dimensional shock-tube problems, the flows past two-dimensional circular cylinder, and the flows around three-dimensional sphere and spacecraft shape with various Knudsen numbers are simulated. The computational results are found in high resolution of the flow fields and good agreement with the theoretical, DSMC, N-S, and experimental results.