Purpose

LA_GELSY computes the minimum-norm least squares solution to one or more real or complex linear systems $A x = b$ using a complete orthogonal factorization of $A$. Matrix $A$ is rectangular and may be rank-deficient. The vectors $b$ and corresponding solution vectors $x$ are the columns of matrices denoted $B$ and $X$, respectively.
The routine computes a $QR$ factorization of $A$ with column pivoting:

$$A \; P = Q \; \left( \begin{array}{cc} R_{11} & R_{12} \\ 0 & R_{22} \end{array} \right)$$

where $R_{11}$ is the largest leading submatrix whose estimated condition number is less than $1/{\bf RCOND}$. The order of $R_{11}$, ${\bf RANK}$, is the effective rank of $A$. $R_{22}$ is considered to be negligible, and $R_{12}$ is annihilated by orthogonal (unitary) transformations from the right, yielding the complete orthogonal (unitary) factorization

$$A \; P = Q \; \left( \begin{array}{cc} T_{11} & 0 \\ 0 & 0 \end{array} \right) Z$$

The minimum-norm least squares solution is then

$$x = P \; Z^H \left( \begin{array}{cc} T_{11}^{-1} Q_1^H b \\ 0 \end{array} \right)$$

where $Q_1$ consists of the first ${\bf RANK}$ columns of $Q$.