docs: 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