The efficiency of three Krylov subspace methods with their ILU0-preconditioned version in solving the systems with the nondiagonal sparse matrix is examined. The systems have arisen from the discretization of Poisson’s equation using the 4th and 6th-order compact schemes. Four matrix-vector multiplication techniques based on four sparse matrix storage schemes are considered in the algorithm of the Krylov subspace methods and their effects are explored. The convergence history, error reduction, iteration-resolution relation and CPU-time are addressed. The efficacy of various methods is evaluated against a benchmark scenario in which the conventional second-order central difference scheme is employed to discretize Poisson’s equation. The Krylov subspace methods, paired with four distinct matrix-vector multiplication strategies across three discretization approaches, are tested and implemented within an incompressible fluid flow solver to solve the elliptic segment of the equations. The resulting solution process CPU-time surface gives a new vision regarding speeding up a CFD code with proper selection of discretization stencil and matrix-vector multiplication technique.