  This manuscript presents an efficient solver for the linear system that
arises from the Hierarchical Poincar\'e-Steklov (HPS) discretization of three
dimensional variable coefficient Helmholtz problems. Previous work on the HPS
method has tied it with a direct solver. This work is the first efficient
iterative solver for the linear system that results from the HPS
discretization. The solution technique utilizes GMRES coupled with a locally
homogenized block-Jacobi preconditioner. The local nature of the discretization
and preconditioner naturally yield the matrix-free application of the linear
system. Numerical results illustrate the performance of the solution technique.
This includes an experiment where a problem approximately 100 wavelengths in
each direction that requires more than a billion unknowns to achieve
approximately 4 digits of accuracy takes less than 20 minutes to solve.
