def onecons_qcqp(z, f, tol=1e-6):
    """ Solves a nonconvex problem
      minimize ||x-z||_2^2
      subject to f(x) = x^T P x + q^T x + r ~ 0
      where the relation ~ is given by f.relop (either <= or ==)
    """

    # if constraint is ineq and z is feasible: z is the solution
    if f.relop == '<=' and f.eval(z) <= 0:
        return z

    if f.eigh is None:
        Psymm = (f.P + f.P.T)/2.
        f.eigh = LA.eigh(np.asarray(Psymm.todense()))

    lmb, Q = f.eigh
    zhat = Q.T.dot(z)
    qhat = Q.T.dot(f.qarray)

    # now solve a transformed problem
    # minimize ||xhat - zhat||_2^2
    # subject to sum(lmb_i xhat_i^2) + qhat^T xhat + r = 0
    # constraint is now equality from
    # complementary slackness
    xhat = lambda nu: -np.divide(nu*qhat-2*zhat, 2*(1+nu*lmb))
    phi = lambda xhat: lmb.dot(np.power(xhat, 2)) + qhat.dot(xhat) + f.r

    s = -np.inf
    e = np.inf
    for l in lmb:
        if l > 0: s = max(s, -1./l)
        if l < 0: e = min(e, -1./l)
    if s == -np.inf:
        s = -1.
        while phi(xhat(s)) <= 0: s *= 2.
    if e == np.inf:
        e = 1.
        while phi(xhat(e)) >= 0: e *= 2.
    while e-s > tol:
        m = (s+e)/2.
        p = phi(xhat(m))
        if p > 0: s = m
        elif p < 0: e = m
        else:
            s = e = m
            break
    nu = (s+e)/2.
    return Q.dot(xhat(nu))