A matrix-free parallel solution method for the three-dimensional heterogeneous Helmholtz equation