def METHOD_NAME(psi, src, t):
    assert src != psi
    t("setup")
    p, mmp, r = g.lattice(src), g.lattice(src), g.lattice(src)
    if prec is not None:
        z = g.lattice(src)
    t("matrix")
    mat(mmp, psi)  # in, out
    t("setup")
    g.axpy(r, -1.0, mmp, src)
    if prec is not None:
        z[:] = 0
        prec(z, r)
        g.copy(p, z)
        cp = g.inner_product(r, z).real
    else:
        g.copy(p, r)
        cp = g.norm2(p)
    ssq = g.norm2(src)
    if ssq == 0.0:
        psi[:] = 0
        return
    rsq = self.eps**2.0 * ssq
    for k in range(self.maxiter):
        c = cp
        t("matrix")
        mat(mmp, p)
        t("inner_product")
        d = g.inner_product(p, mmp).real
        a = c / d
        t("axpy_norm2")
        if prec is not None:
            # c = <r,z>, d = <p,A p>
            g.axpy(r, -a, mmp, r)
            t("prec")
            z[:] = 0
            prec(z, r)
            t("axpy_norm2")
            cp = g.inner_product(r, z).real
        else:
            cp = g.axpy_norm2(r, -a, mmp, r)
        t("linear combination")
        b = cp / c
        psi += a * p
        if prec is not None:
            g.axpy(p, b, p, z)
        else:
            g.axpy(p, b, p, r)
        t("other")
        res = abs(cp)
        self.log_convergence(k, res, rsq)
        if k + 1 >= self.miniter:
            if self.eps_abs is not None and res <= self.eps_abs**2.0:
                self.log(f"converged in {k+1} iterations (absolute criterion)")
                return
            if res <= rsq:
                self.log(f"converged in {k+1} iterations")
                return
    self.log(f"NOT converged in {k+1} iterations;  squared residual {res:e} / {rsq:e}")