Compute the packed pivoting Q-R factorization of a matrix.

Parameters:
a     - An n-by-m matrix, m >= n. This will be *overwritten*
        by this function as described below!
enorm - A Euclidian-norm-computing function.
finfo - A Numpy finfo object.

Returns:
pmut   - An n-element permutation vector
rdiag  - An n-element vector of the diagonal of R
acnorm - An n-element vector of the norms of the rows
         of the input matrix 'a'.

Computes the transposed Q-R factorization of the matrix 'a', with
pivoting, in a packed form, in-place. The packed information can be
used to construct matrices Q and R such that

  A P = R Q or, in Python,
  np.dot(r, q) = a[pmut]

where q is m-by-m and q q^T = ident and r is n-by-m and is lower
triangular. The function _qr_factor_full can compute these
matrices. The packed form of output is all that is used by the main LM
fitting algorithm.

"Pivoting" refers to permuting the rows of 'a' to have their norms in
nonincreasing order. The return value 'pmut' maps the unpermuted rows
of 'a' to permuted rows. That is, the norms of the rows of a[pmut] are
in nonincreasing order.

The parameter 'a' is overwritten by this function. Its new value
should still be interpreted as an n-by-m array. It comes in two
parts. Its strict lower triangular part contains the struct lower
triangular part of R. (The diagonal of R is returned in 'rdiag' and
the strict upper trapezoidal part of R is zero.) The upper trapezoidal
part of 'a' contains Q as factorized into a series of Householder
transformation vectors. Q can be reconstructed as the matrix product
of n Householder matrices, where the i'th Householder matrix is
defined by

H_i = I - 2 (v^T v) / (v v^T)

where 'v' is the pmut[i]'th row of 'a' with its strict lower
triangular part set to zero. See _qr_factor_full for more information.

'rdiag' contains the diagonal part of the R matrix, taking into
account the permutation of 'a'. The strict lower triangular part of R
is stored in 'a' *with permutation*, so that the i'th row of R has
rdiag[i] as its diagonal and a[pmut[i],:i] as its upper part. See
_qr_factor_full for more information.

'acnorm' contains the norms of the rows of the original input
matrix 'a' without permutation.

The form of this transformation and the method of pivoting first
appeared in Linpack.

def _qr_factor_packed(a, enorm, finfo):
    """Compute the packed pivoting Q-R factorization of a matrix.

Parameters:
a     - An n-by-m matrix, m >= n. This will be *overwritten*
        by this function as described below!
enorm - A Euclidian-norm-computing function.
finfo - A Numpy finfo object.

Returns:
pmut   - An n-element permutation vector
rdiag  - An n-element vector of the diagonal of R
acnorm - An n-element vector of the norms of the rows
         of the input matrix 'a'.

Computes the transposed Q-R factorization of the matrix 'a', with
pivoting, in a packed form, in-place. The packed information can be
used to construct matrices Q and R such that

  A P = R Q or, in Python,
  np.dot(r, q) = a[pmut]

where q is m-by-m and q q^T = ident and r is n-by-m and is lower
triangular. The function _qr_factor_full can compute these
matrices. The packed form of output is all that is used by the main LM
fitting algorithm.

"Pivoting" refers to permuting the rows of 'a' to have their norms in
nonincreasing order. The return value 'pmut' maps the unpermuted rows
of 'a' to permuted rows. That is, the norms of the rows of a[pmut] are
in nonincreasing order.

The parameter 'a' is overwritten by this function. Its new value
should still be interpreted as an n-by-m array. It comes in two
parts. Its strict lower triangular part contains the struct lower
triangular part of R. (The diagonal of R is returned in 'rdiag' and
the strict upper trapezoidal part of R is zero.) The upper trapezoidal
part of 'a' contains Q as factorized into a series of Householder
transformation vectors. Q can be reconstructed as the matrix product
of n Householder matrices, where the i'th Householder matrix is
defined by

H_i = I - 2 (v^T v) / (v v^T)

where 'v' is the pmut[i]'th row of 'a' with its strict lower
triangular part set to zero. See _qr_factor_full for more information.

'rdiag' contains the diagonal part of the R matrix, taking into
account the permutation of 'a'. The strict lower triangular part of R
is stored in 'a' *with permutation*, so that the i'th row of R has
rdiag[i] as its diagonal and a[pmut[i],:i] as its upper part. See
_qr_factor_full for more information.

'acnorm' contains the norms of the rows of the original input
matrix 'a' without permutation.

The form of this transformation and the method of pivoting first
appeared in Linpack."""

    machep = finfo.eps
    n, m = a.shape

    if m < n:
        raise ValueError('"a" must be at least as tall as it is wide')

    acnorm = np.empty(n, finfo.dtype)
    for j in range(n):
        acnorm[j] = enorm(a[j], finfo)

    rdiag = acnorm.copy()
    wa = acnorm.copy()
    pmut = np.arange(n)

    for i in range(n):
        # Find the row of a with the i'th largest norm, and note it in
        # the pivot vector.

        kmax = rdiag[i:].argmax() + i

        if kmax != i:
            temp = pmut[i]
            pmut[i] = pmut[kmax]
            pmut[kmax] = temp

            rdiag[kmax] = rdiag[i]
            wa[kmax] = wa[i]

            temp = a[i].copy()
            a[i] = a[kmax]
            a[kmax] = temp

        # Compute the Householder transformation to reduce the i'th
        # row of A to a multiple of the i'th unit vector.

        ainorm = enorm(a[i,i:], finfo)

        if ainorm == 0:
            rdiag[i] = 0
            continue

        if a[i,i] < 0:
            # Doing this apparently improves FP precision somehow.
            ainorm = -ainorm

        a[i,i:] /= ainorm
        a[i,i] += 1

        # Apply the transformation to the remaining rows and update
        # the norms.

        for j in range(i + 1, n):
            a[j,i:] -= a[i,i:] * np.dot(a[i,i:], a[j,i:]) / a[i,i]

            if rdiag[j] != 0:
                rdiag[j] *= np.sqrt(max(1 - (a[j,i] / rdiag[j])**2, 0))

                if 0.05 * (rdiag[j] / wa[j])**2 <= machep:
                    # What does this do???
                    wa[j] = rdiag[j] = enorm(a[j,i+1:], finfo)

        rdiag[i] = -ainorm

    return pmut, rdiag, acnorm