A finite-element matrix method of the dynamical equations for constrained flexible multibody systems undergoing large rotations is presented. The algorithmic procedure is based on recursive formulation where all the kinematical expressions as well as the final governing equations of motion are in a matrix form suited for real time computation. The advantages of the method stem from the partitioning of the hybrid coordinates used in the analysis and the identification of the matrices associated with the partial velocities. A parallel processing algorithm based on the procedures outlined is also developed. A discussion on the parallel processors implementation and their utility in simulation of complex constrained multibody systems is presented. A spatial robotic manipulator is simulated to illustrate the performance of the algorithm.