Solving linear systems by orthogonal tridiagonalization (GMINRES and/or GLSQR)

A general matrix A can be reduced to tridiagonal form by orthogonal
transformations on the left and right: UTAV = T. We can arrange that the
rst columns of U and V are proportional to given vectors b and c. An iterative
form of this process was given by Saunders, Simon, and Yip (SINUM 1988) and
used to solve square systems Ax = b and ATy = c simultaneously. (One of the
resulting solvers becomes MINRES when A is symmetric and b = c.)

The approach was rediscovered by Reichel and Ye (NLAA 2008) with emphasis
on rectangular A and least-squares problems Ax ~ b. The resulting solver was
regarded as a generalization of LSQR (although it doesn't become LSQR in
any special case). Careful choice of c was shown to improve convergence.

In his last year of life, Gene Golub became interested in \GLSQR" for
estimating cTx = bTy without computing x or y. Golub, Stoll, and Wathen
(ETNA 2008) revealed that the orthogonal tridiagonalization is equivalent to a
certain block Lanczos process. This reminds us of Golub, Luk, and Overton
(TOMS 1981): a block Lanczos approach to computing singular vectors.