A primitive variable spectral method for simulating incompressible viscous flows inside a finite cylinder is presented. One element of originality of the proposed method is that the radial discretization of the Fourier coefficients depends on the Fourier mode, its dimension decreasing with the increase of the azimuthal modal number. This principle was introduced independently by Matsushima and Marcus and by Verkley in polar coordinates and is adopted here for the first time to formulate a 3D cylindrical Galerkin projection method. A second element of originality is the use of a special basis of Jacobi polynomials introduced recently for the radial dependence in the solution of Dirichlet problems. In this basis the radial operators are represented by matrices of minimal sparsity — diagonal stiffness and tridiagonal mass — provided here in closed form for the first time, and lead to a Helmholtz operator characterized by a favorable condition number. Finally, a new method is presented for eliminating the singular behaviour of the solution originated by the rotation of the lid with respect to the cylindrical wall. Thanks to these elements, the resulting Navier-Stokes spectral solver guarantees the differentiability to any order of the solution in the entire computational domain and does not suffer from the time-step stability restriction occurring in spectral methods with a point clustering close to the axis. Several test examples are offered that demonstrate the spectral accuracy of the solution method under different representative conditions.