**LA_GELSY** computes the minimum-norm least squares solution to one or
more real or complex linear systems using a complete orthogonal
factorization of . Matrix is rectangular and may be rank-deficient.
The vectors and corresponding solution vectors are
the columns of matrices denoted and , respectively.

The routine computes a factorization of with column pivoting:

where is the largest leading submatrix whose estimated condition number is less than . The order of , , is the effective rank of . is considered to be negligible, and is annihilated by orthogonal (unitary) transformations from the right, yielding the complete orthogonal (unitary) factorization

The minimum-norm least squares solution is then

where consists of the first columns of .