def METHOD_NAME(backend):
    from qibo.quantum_info import random_density_matrix
    initial_rho = random_density_matrix(2**3, backend=backend)
    c = Circuit(3, density_matrix=True)
    c.add(gates.H(0))
    c.add(gates.H(1))
    c.add(gates.CNOT(0, 1))
    c.add(gates.H(2))
    final_rho = backend.execute_circuit(c, np.copy(initial_rho))
    h = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
    cnot = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
    m1 = np.kron(np.kron(h, h), np.eye(2))
    m2 = np.kron(cnot, np.eye(2))
    m3 = np.kron(np.eye(4), h)
    m1 = backend.cast(m1, dtype=m1.dtype)
    m2 = backend.cast(m2, dtype=m2.dtype)
    m3 = backend.cast(m3, dtype=m3.dtype)
    target_rho = np.dot(m1, np.dot(initial_rho, np.transpose(np.conj(m1))))
    target_rho = np.dot(m2, np.dot(target_rho, np.transpose(np.conj(m2))))
    target_rho = np.dot(m3, np.dot(target_rho, np.transpose(np.conj(m3))))
    backend.assert_allclose(final_rho, target_rho)