In contact mechanics, enforcing non-penetration constraints through Lagrange multipliers leads to an unconstrained optimization problem, whose mixed finite elements discretization naturally yields saddle-point linear systems, coupling displacement unknowns and active contact constraints. These 2 × 2 block linear systems are difficult to solve being both indefinite and ill-conditioned, hence the paramount importance of their preconditioning. In our work, we focus on developing efficient numerical methods to solve saddle-point systems arising from large-scale 3D contact problems on large parallel computers. The proposed method is implemented using the PETSc library and is intended to be included in MANTA, a new HPC simulation tool developed at CEA. The approach proposed here is based on a projection method described in [Ainsworth, 2001], which transforms the indefinite saddle-point system into a positive definite reduced system. This modified system involves the dense inverse of a submatrix related to the constraint block, leading to significant memory demands. This requires both a matrix-free implementation of the coefficient matrix and a sparse approximation of this inverse in the preconditioning matrix. Our first contribution is hence to propose several preconditioning matrices, with the goal to have them both cheap to construct and light on memory, yet accurate enough to be efficient in conjunction with the algorithm chosen for the preconditioning. Our second contribution, based on [Golub and Greif, 2003], is to introduce a scaling factor in front of one term of the system, that greatly improves the conditioning of the preconditioned system. Several issues related to the parallel distribution of the coefficient and preconditioning matrices are also addressed. Using the Conjugate Gradient method as the iterative linear solver with an Algebraic Multigrid Method for the preconditioning, the current numerical tests show great scalability potential of the method.