This paper describes an approach to numerically approximate the time evolution of multibody systems with flexible (compliant) components. Its salient attribute is that at each time step, both the formulation of the system equations of motion and their numerical solution are carried out using parallel computing on graphics processing unit cards. The equations of motion are obtained using the absolute nodal coordinate formulation, yet any other multibody dynamics formalism would fit equally well the overall solution strategy outlined herein. The implicit numerical integration method adopted relies on a Newton-Krylov methodology and a parallel direct sparse solver to precondition the underlying linear system. The proposed approach, implemented in a software infrastructure available under an open-source BSD-3 license, leads to improvements in overall simulation times of up to one order of magnitude when compared with matrix-free parallel solution approaches that do not use preconditioning.