We discuss the numerical solution of partial differential equations in a particular class of three-dimensional geometries; the two-dimensional cross section (in the xy-plane) can have a general shape, but is assumed to be invariant with respect to the third direction. Earlier work has exploited such geometries by approximating the solution as a truncated Fourier series in the z-direction. In this paper we propose a new solution algorithm which also exploits the tensor-product feature between the xy-plane and the z-direction. However, the new algorithm is not limited to periodic boundary conditions, but works for general Dirichlet and Neumann type of boundary conditions. The proposed algorithm also works for problems with variable coefficients as long as these can be expressed as a separable function with respect to the variation in the xy-plane and the variation in the z-direction. For problems where the new method is applicable, the computational cost is very competitive with the best iterative solvers. The new algorithm is easy to implement, and useful, both in a serial and parallel context. Numerical results demonstrating the superiority of the method are presented for three-dimensional Poisson and Helmholtz problems using both low order finite elements and high order spectral element discretizations.