Sequential numerical solution methods are commonly used for solving the phonon Boltzmann transport equation (BTE) because of simplicity of implementation and low storage requirements. However, they exhibit poor convergence for low Knudsen numbers. This is because sequential solution procedures couple the phonon BTEs in physical space efficiently but the coupling is inefficient in wave vector (**K**) space. As the Knudsen number decreases, coupling in **K** space becomes dominant and convergence rates fall. Since materials like silicon have **K**-resolved Knudsen numbers that span two to five orders of magnitude at room temperature, diffuse-limit solutions are not feasible for all **K** vectors. Consequently, nongray solutions of the BTE experience extremely slow convergence. In this paper, we develop a coupled-ordinates method for numerically solving the phonon BTE in the relaxation time approximation. Here, interequation coupling is treated implicitly through a point-coupled direct solution of the **K**-resolved BTEs at each control volume. This implicit solution is used as a relaxation sweep in a geometric multigrid method which promotes coupling in physical space. The solution procedure is benchmarked against a traditional sequential solution procedure for thermal transport in silicon. Significant acceleration in computational time, between 10 and 300 times, over the sequential procedure is found for heat conduction problems.