def matrix_free_solve(matvec, b):
    return lax.custom_linear_solve(matvec, b, richardson_iteration, richardson_iteration)
def matrix_free_solve(matvec, b):
    return lax.custom_linear_solve(matvec, b, richardson_iteration, richardson_iteration)