def METHOD_NAME(A, l, b, complex_evp):
    """Compute an Arnoldi factorization."""
    v = b * (1 / b.norm()[0])
    H = np.zeros((l, l), dtype=np.complex_ if complex_evp else np.float_)
    V = A.source.empty(reserve=l)
    V.append(v)
    for i in range(l):
        v = A.apply(v)
        V.append(v)
        _, R = gram_schmidt(V, return_R=True, atol=0, rtol=0, offset=len(V) - 1, copy=False)
        H[:i + 2, i] = R[:l, i + 1]
        v = V[-1]
    return V[:l], H, v * R[l, l]