Summary: Complementarity solvers are continually being challenged by modelers demanding improved reliability and scalability. Building upon a strong theoretical background, the semismooth algorithm has the potential to meet both of these requirements. We discuss relevant theory associated with the algorithm and then describe a sophisticated implementation in detail. Particular emphasis is given to the use of preconditioned iterative methods to solve the (nonsymmetric) systems of linear equations generated at each iteration and robust methods for dealing with singularity. Results on the MCPLIB test suite indicate that the code is reliable and efficient and scales well to very large problems.