Compute the Levenberg-Marquardt parameter and solution vector.

Parameters:
a     - n-by-m matrix, m >= n (only the n-by-n component is used)
b     - n-by-n matrix
ddiag - n-vector, diagonal elements of D
delta - positive scalar, specifies scale of enorm(Dx)
par0  - positive scalar, initial estimate of the LM parameter

Returns:
par    - positive scalar, final estimate of LM parameter
x      - n-vector, least-squares solution of LM equation
dxnorm - positive scalar, enorm(D x)
relnormdiff - scalar, (dxnorm - delta) / delta, maybe abs-ified

This routine computes the Levenberg-Marquardt parameter 'par' and a LM
solution vector 'x'. Given an n-by-n matrix A, an n-by-n nonsingular
diagonal matrix D, an m-vector B, and a positive number delta, the
problem is to determine values such that 'x' is the least-squares
solution to

 A x = B
 sqrt(par) * D x = 0

and either

 (1) par = 0, dxnorm - delta <= 0.1 delta or
 (2) par > 0 and |dxnorm - delta| <= 0.1 delta

where dxnorm = enorm(D x).

def _lm_solve_full(a, b, ddiag, delta, par0, dtype=np.float):
    """Compute the Levenberg-Marquardt parameter and solution vector.

Parameters:
a     - n-by-m matrix, m >= n (only the n-by-n component is used)
b     - n-by-n matrix
ddiag - n-vector, diagonal elements of D
delta - positive scalar, specifies scale of enorm(Dx)
par0  - positive scalar, initial estimate of the LM parameter

Returns:
par    - positive scalar, final estimate of LM parameter
x      - n-vector, least-squares solution of LM equation
dxnorm - positive scalar, enorm(D x)
relnormdiff - scalar, (dxnorm - delta) / delta, maybe abs-ified

This routine computes the Levenberg-Marquardt parameter 'par' and a LM
solution vector 'x'. Given an n-by-n matrix A, an n-by-n nonsingular
diagonal matrix D, an m-vector B, and a positive number delta, the
problem is to determine values such that 'x' is the least-squares
solution to

 A x = B
 sqrt(par) * D x = 0

and either

 (1) par = 0, dxnorm - delta <= 0.1 delta or
 (2) par > 0 and |dxnorm - delta| <= 0.1 delta

where dxnorm = enorm(D x).
"""
    a = np.asarray(a, dtype)
    b = np.asarray(b, dtype)
    ddiag = np.asarray(ddiag, dtype)

    n, m = a.shape
    assert m >= n
    assert b.shape == (m, )
    assert ddiag.shape == (n, )

    q, r, pmut = _qr_factor_full(a)
    bqt = np.dot(b, q.T)
    par, x = _lm_solve(r, pmut, ddiag, bqt, delta, par0,
                       enorm_mpfit_careful, np.finfo(dtype))
    dxnorm = enorm_mpfit_careful(ddiag * x, np.finfo(dtype))
    relnormdiff = (dxnorm - delta) / delta

    if par > 0:
        relnormdiff = abs(relnormdiff)

    return par, x, dxnorm, relnormdiff