In this thesis we will show that unsteady, incompressible fluid-flow problems, as described by the three-dimensional Navier-Stokes equations, can efficiently be simulated by a parallel computer code based on a spectral element discretization. This efficiency relies a great deal on the numerous improvements to spectral methods on the theoretical, algorithmic and computational levels which have been proposed over the last decade. Consequently, the application of spectral methods is no longer limited to regular problems in simple geometries, but is extended to complex situations. The spectral element method in particular has received much attention, because it combines the practical advantages of the well-established finite element method with the "spectral" ability to reduce the number of degrees of freedom to obtain a prescribed level of accuracy. Furthermore, it embodies the concept of domain decomposition which is closely related to parallelism. The latter aspect is important for an effective implementation on the leading ("affordable") supercomputers of the nineties, which are almost exclusively based on parallel processors with distributed memory. The continuous Navier-Stokes equations describe the velocity and pressure of a fluid in a domain and the spectral element method reduces them to a set of discrete equations. This discretization technique has become rather mature and furnishes the framework for this thesis. The main challenge, however, lies in the development of fast algorithms for the solution of the discrete systems. To this end, we will analyze existing and new methods of decoupling the velocity components from the pressure. Iterative methods have shown their merits for the solution of the resulting, decoupled set of equations, provided that they are properly preconditioned. This condition has led to the introduction of a large number of preconditioning methods for velocity and pressure operators. The common factor of these techniques is that they are based on fast, local solves, combined with a strategy to deal with the element interfaces. Another important issue is the time discretization of the unsteady equations. It is common practice to choose an implicit time-integration scheme for the linear terms, and an explicit one for the nonlinear terms. The way in which these two schemes are combined is not trivial, especially when two important conditions, high-order accuracy in time and stability, have to be satisfied simultaneously. In fact, the concepts of precision and stability are related: the more stable a time scheme, the larger the maximum time step that can be used. From a computational viewpoint, it is desirable to advance in time with the maximum dowed time step. An acceptable level of accuracy in time is then only guaranteed by high-order time schemes. The construction of a high-performance code requires not only readily parallelizable numerical solution methods, but also fast algorithms to manage the (little) communication involved. The use of iterative methods for the solution of the discrete systems enhances the parallel implementation, especially when the preconditioners are block diagonal (communication free). The communication algorithms are analyzed and optimized, and the use of very fast, architecture-specific routines has to be balanced against the portability of the code. All the algorithmic, numerical, computational, and parallel improvements are first tested individually on small problems. The parallel spectral element code is based on the outcome of these tests and is used for the simulation of the three-dimensional flow over a backward-facing step.