Optimal flow control problems are important for applications in science and engineering. Solving such problems usually requires the solution of a large linear generalized saddle-point system. This linear system is sparse and highly indefinite. In order to solve such systems using Krylov subspace methods, efficient preconditioners are necessary to enhance their robustness and accelerate the convergence. Standard block preconditioning techniques for the generalized saddle-point systems require an efficient approximation of the Schur complement. This is a big challenge since the Schur complement is large and dense, and therefore computationally expensive to approximate. For some problems, it is even impossible to approximate the Schur complement efficiently. In this dissertation, we propose a new class of preconditioners for optimal flow control problems using multilevel sequentially semiseparable (MSSS) matrix computations. In contrast to standard block preconditioners, MSSS preconditioners do not approximate the Schur complement of the generalized saddle-point system but compute an approximate factorization of the global block system in linear computational complexity. This is a big advantage over block preconditioners. The key to this global factorization is that the Schur complements in this factorization usually have low off-diagonal rank. Therefore, these Schur complements can be approximated by matrices with a low rank off-diagonal structure. For this, MSSS matrix computations are very well suited. Theoretical analysis shows that MSSS preconditioners yield a spectrum of the preconditioned system matrix that is contained in a circle centered at $(1, 0)$. This radius can be controlled arbitrarily by properly choosing a parameter in the MSSS preconditioner computations. This in turn implies that the convergence of MSSS preconditioned systems can be independent of the mesh size and regularization parameter for PDE-constrained optimization problems while for computational fluid dynamics problems, the convergence is independent of the mesh size and Reynolds number. Mesh size independent and wave number independent convergence can be also obtained when applying the MSSS preconditioning technique to the Helmholtz problem. Numerical results verify the convergence property. In this dissertation, we also studied the problem of optimal in-domain control of the Navier-Stokes equation. We use a simplified wind farm control example to formulate such a problem. Compared with standard PDE-constrained optimization problems where the controls are distributed throughout the whole domain, this in-domain control problem is even more difficult to solve since the control only acts on a few parts of the domain. This in turn gives a linear system of the generalized saddle-point type. Block preconditioners cannot give satisfactory performance for such problem because the Schur complement for such system is very difficult or even impossible to approximate efficiently. Applying MSSS preconditioners to this problem gives superior performance compared to block preconditioning techniques.