A numerical method for computing three-dimensional Stokes flow driven by a doubly-periodic array of regularized forces is presented. In the non-periodic direction either a free boundary or a homogeneous Dirichlet condition is enforced. The method consists of finding a regularized Green's function in Fourier space analytically. Then only an inverse fast Fourier transform (inverse FFT) has to be computed. Accuracy is verified by comparing numerical results to a solution that is independent of the method. In an Ewald splitting, the FFT method can be used to compute the smooth component of the flow, which allows for a splitting parameter as small as a few grid cells. This selection makes the sum in physical space converge extremely fast. Numerical examples demonstrate that fact. Since the forces are regularized, in some cases splitting is not even needed, depending on the relative sizes of the numerical parameters. The method is applied to model the flow created by carpets of nodal cilia based on cilium shape.