I finally implemented Baraff's algorithm from the paper, "Linear-Time Dynamics using Lagrange Multipliers". I have compared the results of the Baraff method with those of the original method (Smith's ODE) using the piston assembly simulation. The holy grail seems to be sparse LCP, but now I am content with the current implementation. According to the paper, the idea is to first solve for the primary constraints using the sparse LDLT decomposition, and then solve for the auxiliary constraints using the LCP method. The followings are the essential steps.
A = J2 M - 1 K 
Given the external force F and the velocity vt at the current time step, the goal is to find the velocity at the next time step, vt + h, where h is the time step size. The system equation can be formulated as
where M is the system inertia matrix, J1 is the primary constraint Jacobian, J2 is the auxiliary constraint Jacobian, λ is the primary constraint impulse, and μ is the auxiliary constraint impulse.
First, solve
using sparse LDLT decomposition to yield
Next, formulate A as
where K is defined as
where B is solved from the following equation
Note that the LDLT factorization has been done for the matrix
in the previous step, so solving for B is rather efficient.
Finally, formulate and solve the following equation for μ using the LCP method.
The velocity vt + h is given by
Comments