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