code: def _qrd_solve_full(a, b, ddiag, dtype=np.float):
      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)
      x, s = _manual_qrd_solve(r[:,:n], pmut, ddiag, bqt,
                               dtype=dtype, build_s=True)
      return x, s, pmut