  We introduce innovative algorithms for computing exact or approximate
(minimum-norm) solutions to $Ax=b$ or the {\it normal equation} $A^TAx=A^Tb$,
where $A$ is an $m \times n$ real matrix of arbitrary rank. We present more
efficient algorithms when $A$ is symmetric PSD. First, we introduce the {\it
Triangle Algorithm} (TA), a {\it convex-hull membership} algorithm that given
$b_k=Ax_k$ in the ellipsoid $E_{A,\rho}=\{Ax: \Vert x \Vert \leq \rho\}$, it
either computes an improved approximation $b_{k+1}=Ax_{k+1}$ or proves $b \not
\in E_{A,\rho}$. We then give a dynamic variant of TA, the {\it Centering
Triangle Algorithm} (CTA), generating residual, $r_k=b -Ax_k$ via the iteration
of $F_1(r)=r-(r^THr/r^TH^2r)Hr$, where $H=AA^T$. If $A$ is symmetric PSD, $H$
can be taken as $A$. Next, for each $t=1, \dots, m$, we derive $F_t(r)=r-
\sum_{i=1}^t \alpha_{t,i}(r) H^i r$ whose iterations correspond to a Krylov
subspace method with restart. If $\kappa^+(H)$ is the ratio of the largest to
smallest positive eigenvalues of $H$, when $Ax=b$ is consistent, in
$k=O({\kappa^+(H)}{t^{-1}} \ln \varepsilon^{-1})$ iterations of $F_t$, $\Vert
r_k \Vert \leq \varepsilon$. Each iteration takes $O(tN+t^3)$ operations, $N$
the number of nonzero entries in $A$. By directly applying $F_t$ to the normal
equation, we get $\Vert A^TAx_k - A^Tb \Vert \leq \varepsilon$ in
$O({\kappa^+(AA^T)}{t}^{-1} \ln \varepsilon^{-1})$ iterations. On the other
hand, given any residual $r$, we compute $s$, the degree of its minimal
polynomial with respect to $H$ in $O(sN+s^3)$ operations. Then $F_s(r)$ gives
the minimum-norm solution of $Ax=b$ or an exact solution of $A^TAx=A^Tb$. The
proposed algorithms are simple to implementation and theoretically robust. We
present sample computational results, comparing the performance of CTA with CG
and GMRES methods. The results support CTA as a highly competitive option.
