docs: def _qrd_solve_full(a, b, ddiag, dtype=np.float):
      """Solve the equation A^T x = B, D x = 0.
  
  Parameters:
  a     - an n-by-m array, m >= n
  b     - an m-vector
  ddiag - an n-vector giving the diagonal of D. (The rest of D is 0.)
  
  Returns:
  x    - n-vector solving the equation.
  s    - the n-by-n supplementary matrix s.
  pmut - n-element permutation vector defining the permutation matrix P.
  
  The equations are solved in a least-squares sense if the system is
  rank-deficient.  D is a diagonal matrix and hence only its diagonal is
  in fact supplied as an argument. The matrix s is full lower triangular
  and solves the equation
  
  P^T (A A^T + D D) P = S^T S (needs transposition?)
  
  where P is the permutation matrix defined by the vector pmut; it puts
  the rows of 'a' in order of nonincreasing rank, so that a[pmut]
  has its rows sorted that way.
  """
  
      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, )
  
      # The computation is straightforward.
  
      q, r, pmut = _qr_factor_full(a)
      bqt = np.dot(b, q.T)
      x, s = _manual_qrd_solve(r[:,:n], pmut, ddiag, bqt,
                               dtype=dtype, build_s=True)
  
      return x, s, pmut