We study the solution of neutral delay differential equations (NDDEs) by using boundary value methods (BVMs). The BVMs require the solution of nonsymmetric, large and sparse linear systems. The GMRES method with the Strang-type block-circulant preconditioner is proposed to solve these linear systems. We show that if an Ak1,k2 -stable BVM is used for solving an m-by-m system of NDDEs, then our preconditioner is invertible and the spectrum of the preconditioned system is clustered. It follows that when the GMRES method is applied to the preconditioned systems, the method could converge fast. Numerical results are given to show that our method is effective.