Compute the QR factorization of a matrix, with pivoting.

Parameters:
a     - An n-by-m arraylike, m >= n.
dtype - (optional) The data type to use for computations.
        Default is np.float.

Returns:
q    - An m-by-m orthogonal matrix (q q^T = ident)
r    - An n-by-m upper triangular matrix
pmut - An n-element permutation vector

The returned values will satisfy the equation

np.dot(r, q) == a[:,pmut]

The outputs are computed indirectly via the function
_qr_factor_packed. If you need to compute q and r matrices in
production code, there are faster ways to do it. This function is for
testing _qr_factor_packed.

The permutation vector pmut is a vector of the integers 0 through
n-1. It sorts the rows of 'a' by their norms, so that the
pmut[i]'th row of 'a' has the i'th biggest norm.

def _qr_factor_full(a, dtype=np.float):
    """Compute the QR factorization of a matrix, with pivoting.

Parameters:
a     - An n-by-m arraylike, m >= n.
dtype - (optional) The data type to use for computations.
        Default is np.float.

Returns:
q    - An m-by-m orthogonal matrix (q q^T = ident)
r    - An n-by-m upper triangular matrix
pmut - An n-element permutation vector

The returned values will satisfy the equation

np.dot(r, q) == a[:,pmut]

The outputs are computed indirectly via the function
_qr_factor_packed. If you need to compute q and r matrices in
production code, there are faster ways to do it. This function is for
testing _qr_factor_packed.

The permutation vector pmut is a vector of the integers 0 through
n-1. It sorts the rows of 'a' by their norms, so that the
pmut[i]'th row of 'a' has the i'th biggest norm."""

    n, m = a.shape

    # Compute the packed Q and R matrix information.

    packed, pmut, rdiag, acnorm = \
        _manual_qr_factor_packed(a, dtype)

    # Now we unpack. Start with the R matrix, which is easy: we just
    # have to piece it together from the strict lower triangle of 'a'
    # and the diagonal in 'rdiag'.

    r = np.zeros((n, m))

    for i in range(n):
        r[i,:i] = packed[i,:i]
        r[i,i] = rdiag[i]

    # Now the Q matrix. It is the concatenation of n Householder
    # transformations, each of which is defined by a row in the upper
    # trapezoidal portion of 'a'. We extract the appropriate vector,
    # construct the matrix for the Householder transform, and build up
    # the Q matrix.

    q = np.eye(m)
    v = np.empty(m)

    for i in range(n):
        v[:] = packed[i]
        v[:i] = 0

        hhm = np.eye(m) - 2 * np.outer(v, v) / np.dot(v, v)
        q = np.dot(hhm, q)

    return q, r, pmut