In this paper, we consider an efficient iterative approach to the solution of the discrete Helmholtz equation with Dirichlet, Neumann and Sommerfeld-like boundary conditions based on a compact sixth order approximation scheme and preconditioned Krylov subspace methodology. A sixth order compact scheme for the 3D Helmholtz equation with different boundary conditions is developed to reduce approximation and pollution errors, thereby softening the point-per-wavelength constraint. The resulting systems of finite-difference equations are solved by different preconditioned Krylov subspace-based methods. In the majority of test problems, the preconditioned Generalized Minimal Residual (GMRES) method is the superior choice, but in the case of sufficiently fine grids a simple stationary two-level algorithm proposed in this paper in combination with a lower order approximation preconditioner presents an efficient alternative to the GMRES method. In the analysis of the lower order preconditioning developed here, we introduce the term "\(k\)-th order preconditioned matrix" in addition to the commonly used "an optimal preconditioner". The necessity of the new criterion is justified by the fact that the condition number of the preconditioned matrix \( AA^{-1}_p \) in some of our test problems improves with the decrease of the grid step size. In a simple 1D case, we are able to prove this analytically. This new parameter could serve as a guide in the construction of new preconditioners. The lower order direct preconditioner used in our algorithms is based on a combination of the separation of variables technique and Fast Fourier Transform (FFT) type methods. The resulting numerical methods allow efficient implementation on parallel computers. Numerical results confirm the high efficiency of the proposed iterative approach.