This work discusses linear systems numerical solution arising from the discretized Navier-Stokes equations in fluid mechanics applications. These systems exhibit a saddle point structure and solving them is essential for accurately simulating fluid flows. Traditional methods employ a Navier-Stokes-based root-finding algorithm which involves solving a saddle point system at each Newton iteration. Due to the large scale and ill-conditioning of the linear systems involved, efficient linear solvers are necessary. In this work, a new approach based on the restarted regularized preconditioned generalized minimal residual solver is developed. The MUltifrontal Massively Parallel Sparse direct Solver (MUMPS) is proposed to accelerate and enhance the performance of the regularized preconditioner. The effectiveness of various preconditioners is assessed through numerical experiments on the Lid-driven cavity and Flow over backward facing step systems. The performance of the restarted solvers is evaluated based on computational time and the number of restarted GMRES iterations, providing insights into their efficiency and effectiveness.