  Consider solving large sparse range symmetric singular linear systems $ A
{\bf x}= {\bf b} $ which arise, for instance, in the discretization of
convection diffusion equations with periodic boundary conditions, and partial
differential equations for electromagnetic fields using the edge-based finite
element method.
  In theory, the Generalized Minimal Residual (GMRES) method converges to the
least squares solution for inconsistent systems if the coefficient matrix $A$
is range symmetric, i.e. $ {\rm R}(A)= {\rm R}(A^{ \rm T } )$, where $ {\rm
R}(A)$ is the range space of $A$.
  We derived the necessary and sufficient conditions for GMRES to determine a
least squares solution of inconsistent and consistent range symmetric systems
assuming exact arithmetic except for the computation of the elements of the
Hessenberg matrix.
  In practice, GMRES may not converge due to numerical instability. In order to
improve the convergence, we propose using the pseudoinverse for the solution of
the severely ill-conditioned Hessenberg systems in GMRES. Numerical experiments
on inconsistent systems indicate that the method is effective and robust.
Finally, we further improve the convergence of the method by reorthogonalizing
the Modified Gram-Schmidt procedure.
