We study a parallel Newton-Krylov-Schwarz (NKS) based algorithm for solving large sparse systems resulting from a fully implicit discretization of partial differential equations arising from petroleum reservoir simulations. Our NKS algorithm is designed by combining an inexact Newton method with a rank-2 updated quasi-Newton method. In order to improve the computational efficiency, both DDM and SPMD parallelism strategies are adopted. The effectiveness of the overall algorithm depends heavily on the performance of the linear preconditioner, which is made of a combination of several preconditioning components including AMC, relaxed ILU, up scaling, additive Schwarz, CRP-like (constraint residual preconditioning), Watts correction, Shur complement, among others. In the construction of the CRP-like preconditioner, a restarted CMRES is used to solve the projected linear systems. We have tested this algorithm and related parallel software using data from some real applications, and presented numerical results that show that this solver is robust and scalable for large scale calculations in petroleum reservoir simulations.