  Brown and Walker (1997) showed that GMRES determines a least squares solution
of $ A x = b $ where $ A \in {\bf R}^{n \times n} $ without breakdown for
arbitrary $ b, x_0 \in {\bf R}^n $ if and only if $A$ is range-symmetric, i.e.
$ {\cal R} (A^{\rm T}) = {\cal R} (A) $, where $ A $ may be singular and $ b $
may not be in the range space ${\cal R} A)$ of $A$. In this paper, we propose
applying GMRES to $ A C A^{\rm T} z = b $, where $ C \in {\bf R}^{n \times n} $
is symmetric positive definite. This determines a least squares solution $ x =
CA^{\rm T} z $ of $ A x = b $ without breakdown for arbitrary (singular) matrix
$A \in {\bf R}^{n \times n}$ and $ b \in {\bf R}^n $. To make the method
numerically stable, we propose using the pseudoinverse with an appropriate
threshold parameter to suppress the influence of tiny singular values when
solving the severely ill-conditioned Hessenberg systems which arise in the
Arnoldi process of GMRES when solving inconsistent range-symmetric systems.
Numerical experiments show that the method taking $C$ to be the identity matrix
and the inverse matrix of the diagonal matrix whose diagonal elements are the
diagonal of $A A^{\rm T}$ gives a least squares solution even when $A$ is not
range-symmetric, including the case when $ {\rm index}(A) >1$.
