def METHOD_NAME():
    """solve(a == L, out) should return the same as solving with the assembled objects.
    This relies on two different code paths agreeing on the same set of solver parameters."""
    mesh = UnitSquareMesh(50, 50)
    V = FunctionSpace(mesh, "CG", 1)
    f = Function(V)
    f.assign(1)
    f.vector()[:] = 1.
    t = TestFunction(V)
    q = TrialFunction(V)
    a = inner(q, t) * dx
    L = inner(f, t) * dx
    # Solve the system using forms
    sol = Function(V)
    solve(a == L, sol)
    # And again
    sol2 = Function(V)
    solve(a == L, sol2)
    assert np_norm(sol.vector()[:] - sol2.vector()[:]) == 0
    # Solve the system using preassembled objects
    sol3 = Function(V)
    solve(assemble(a), sol3, assemble(L))
    assert np_norm(sol.vector()[:] - sol3.vector()[:]) < 5e-14
    # Same, solving into vector
    sol4 = sol3.vector()
    solve(assemble(a), sol4, assemble(L))
    assert np_norm(sol.vector()[:] - sol4[:]) < 5e-14