A new and efficient computational method for constrained multibody systems is proposed. In the proposed method, local parametrization method is employed to apply the same solution method for position, velocity, and acceleration analyses since the coefficient matrices for each analysis have an identical matrix pattern. The skyline solution method is used to overcome numerical inefficiency when solving large scaled equations. Also, subsystem mpartitioning method is derived systematically to perform parallel processing for real time simulation. To show the numerical accuracy and efficiency of the proposed method, three numerical problems are solved.