Large scale optimization of systems governed by partial differential equations (PDEs) is a frontier problem in scientific computation. The state-of-the-art for solving such problems is reduced-space quasi-Newton sequential quadratic programming (SQP) methods. These take full advantage of existing PDE solver technology and parallelize well. However, their algorithmic scalability is questionable; for certain problem classes they can be very slow to converge. In this paper we propose a full-space Newton-Krylov SQP method that uses the reduced-space quasi-Newton method as a preconditioner. The new method is fully parallelizable; exploits the structure of and available parallel algorithms for the PDE forward problem; and is quadratically convergent close to a local minimum. We restrict our attention to boundary value problems and we solve a model optimal flow control problem, with both Stokes and Navier-Stokes equations as constraints. Algorithmic comparisons, scalability results, and parallel performance on a Cray T3E-900 are presented. On the model problems solved, the new method is a factor of 5-10 faster than reduced space quasi-Newton SQP, and is scalable provided a good forward preconditioner is available.