def cg_steihaug(vs, Avp, g, z0, Delta, rtol, maxiter=1000, monitor=None, C=None
    ):
    """ best effort solving for y = A^{-1} g with cg,
        given the trust-region constraint;

        This is roughly ported from Jeff Regier's

            https://github.com/jeff-regier/Celeste.jl/blob/316cae12efcc394a5253ee799056c5513c2efc67/src/cg_trust_region.jl

        the algorithm is almost identical to the one on NWU wiki,

            https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods

        if given B, update approximated Hessian, see Morale and Nocedal, 2000

            http://epubs.siam.org/doi/abs/10.1137/S1052623497327854


        C(v, direction) the preconditioner operator. C^{direction}. 
        Usually, C is close to A, such that C^{-1} A has a lower condition number.

        One particular case is if C is the diagonal of A.

        See Steihaug's paper. https://epubs.siam.org/doi/pdf/10.1137/0720042

    """
    if C is None:
        C = lambda x, direction: x
    dot = vs.dot
    mul = vs.mul
    addmul = vs.addmul
    if z0 is None:
        z0 = vs.zeros_like(g)
    r0 = addmul(Avp(z0), g, -1)
    mr0 = C(r0, -1)
    d0 = mul(mr0, 1)
    j = 0
    rho_init = dot(mr0, r0)
    rho0 = rho_init
    while True:
        Bd0 = Avp(d0)
        dBd0 = dot(d0, Bd0)
        alpha = rho0 / dBd0
        p0 = addmul(z0, d0, -alpha)
        message = ''
        if dBd0 == 0:
            rho1 = 0
            z1 = z0
            r1 = r0
            mr1 = mr0
            d1 = d0
            message = 'zero hessian'
        elif dBd0 <= 0 or dot(p0, C(p0, 1)) ** 0.5 >= Delta:
            a_ = dot(d0, C(d0, 1))
            b_ = 2 * dot(z0, C(d0, 1))
            c_ = dot(z0, C(z0, 1)) - Delta ** 2
            cond = b_ ** 2 - 4 * a_ * c_
            if a_ == 0 or cond < 0:
                rho1 = -1
                tau = 1.0
                if c_ > 0:
                    tau = Delta / c_ ** 0.5
                z1 = mul(z0, tau)
                if a_ == 0:
                    message = 'already at the right direction'
                    rho1 = -1
                else:
                    message = (
                        'no solution to second order equation, restarting ')
                    z0 = vs.zeros_like(g)
                    r0 = addmul(Avp(z0), g, -1)
                    mr0 = C(r0, -1)
                    d0 = mul(mr0, 1)
                    rho1 = dot(mr0, r0)
            else:
                tau = (-b_ - cond ** 0.5) / (2 * a_)
                z1 = addmul(z0, d0, tau)
                if dBd0 <= 0:
                    rho1 = -1
                else:
                    rho1 = 0
                if dBd0 <= 0:
                    message = 'negative curvature '
                else:
                    message = 'truncation'
            r1 = r0
            mr1 = mr0
            d1 = d0
        else:
            z1 = addmul(z0, d0, -alpha)
            r1 = addmul(r0, Bd0, -alpha)
            mr1 = C(r1, -1)
            rho1 = dot(mr1, r1)
            d1 = mul(mr1, 1)
            d1 = addmul(d1, d0, rho1 / rho0)
            message = 'regular iteration'
        r0 = r1
        mr0 = mr1
        d0 = d1
        z0 = z1
        rho0 = rho1
        if monitor is not None:
            z0norm = dot(z0, z0)
            zgnorm = dot(z0, g)
            gnorm = dot(g, g)
            monitor(j, message, rho0, rho_init, rtol, zgnorm / z0norm ** 
                0.5 / gnorm ** 0.5)
        if rho1 / rho_init < rtol ** 2:
            break
        if j >= maxiter:
            break
        j = j + 1
    return z0
