Efficient parallel computation of complex flows is essential to bring modern computer power to bear on fluid calculations where complicated physical descriptions are required. We have developed an efficient, parallel computational method for solving generalized Stokes flow problems that arise when operator splitting of the velocity and pressure fields is used in Newton's method for solution of nonlinear, steady-state flow problems. The key to the parallelization is the incorporation of a preconditioned iterative matrix solution. The linear system that results from finite element discretization of the generalized Stokes problem is asymmetric, indefinite, and block singular. At each Newton step, this system is solved using an algorithm that combines a parallel preconditioner with a Krylov subspace method. The parallel preconditioner, called the block complement and additive levels method (BCALM) preconditioner, is based on treating pressure unknowns separately from the velocities and gradients. A pressure preconditioner is constructed from factorization of the Schur complement of the pressures using a Jacobi-type iteration. The viscous operator is preconditioned using the additive Schwarz method. The resulting iterative method is demonstrated to have high parallel efficiency, subject to effective domain decomposition. The iterative solver is developed in the context of simulation of natural convection modeled by the Boussnesq approximation. For natural convection in a rectangular cavity heated from the sides, in the limit of high Grashof number, the linear system that arises during solution for steady state using Newton's method is stiff, asymmetric and indefinite. For this model problem, the preconditioner is shown to be robust and the overall iterative solution is highly efficient relative to other solution methods.