Steady-state Water Distribution Network models compute pipe flows and nodal heads for assumed nodal demands, pipe hydraulic resistances, etc. The nonlinear mathematical problem is based on energy and mass conservation laws which is solved by using global linearization techniques, such as global gradient algorithm (GGA). The matrix of coefficients of the linear system inside GGA belongs to the class of sparse, symmetric and positive definite. Therefore a fast solver for the linear system is important in order to achieve the computational efficiency, especially when multiple runs are required. This work aims at testing three main strategies for the solution of linear systems inside GGA. The tests are performed on eight real networks by sampling nodal demands, considering the pressure-driven and demand-driven modelling to evaluate the robustness of solvers. The results show that there exists a robust specialized direct method which is superior to all the other alternatives. Furthermore, it is found that the number of times the linear system is solved inside the GGA does not depend on the specific solver, if a small regularization to the linear problem is applied, and that pressure-driven modelling requires a greater number which depends on the size and topology of the network and not only on the level of pressure deficiency.