A Newton-Krylov algorithm is presented for the compressible Navier-Stokes equations in three dimensions on unstructured grids. The algorithm uses a preconditioned matrix-free Krylov method to solve the linear system that arises in the Newton iterations. Incomplete factorization is used as the preconditioner, based on an approximate Jacobian matrix after the reverse Cuthill-McKee reordering of the unknowns. Approximate viscous operators that involve only the nearest neighboring terms are studied to construct an efficient and effective preconditioner. Two incomplete factorization techniques are examined: one based on a level-of-fill approach while the other utilizes a threshold strategy. The performance of the algorithm is demonstrated with numerical studies over the ONERA M6 wing and the DLR-F6 wing-body configuration. A ten-order-of-magnitude residual reduction can be obtained in 20-25 hours on a single processor for a mesh with roughly 500,000 nodes.