diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 6006688fdffe..893634b7b435 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -37,6 +37,7 @@
 sys.modules["qiskit._accelerate.nlayout"] = qiskit._accelerate.nlayout
 sys.modules["qiskit._accelerate.stochastic_swap"] = qiskit._accelerate.stochastic_swap
 sys.modules["qiskit._accelerate.sabre_swap"] = qiskit._accelerate.sabre_swap
+sys.modules["qiskit._accelerate.sabre_layout"] = qiskit._accelerate.sabre_layout
 sys.modules["qiskit._accelerate.pauli_expval"] = qiskit._accelerate.pauli_expval
 sys.modules["qiskit._accelerate.dense_layout"] = qiskit._accelerate.dense_layout
 sys.modules["qiskit._accelerate.sparse_pauli_op"] = qiskit._accelerate.sparse_pauli_op
diff --git a/qiskit/transpiler/passes/layout/sabre_layout.py b/qiskit/transpiler/passes/layout/sabre_layout.py
index 3c3abdaf3dce..30c5f67a8861 100644
--- a/qiskit/transpiler/passes/layout/sabre_layout.py
+++ b/qiskit/transpiler/passes/layout/sabre_layout.py
@@ -13,35 +13,70 @@
 """Layout selection using the SABRE bidirectional search approach from Li et al.
 """
 
+from collections import defaultdict
+import copy
 import logging
+import time
+
 import numpy as np
+import rustworkx as rx
 
 from qiskit.converters import dag_to_circuit
 from qiskit.transpiler.passes.layout.set_layout import SetLayout
 from qiskit.transpiler.passes.layout.full_ancilla_allocation import FullAncillaAllocation
 from qiskit.transpiler.passes.layout.enlarge_with_ancilla import EnlargeWithAncilla
 from qiskit.transpiler.passes.layout.apply_layout import ApplyLayout
-from qiskit.transpiler.passes.routing import SabreSwap
+from qiskit.transpiler.passes.layout import vf2_utils
 from qiskit.transpiler.passmanager import PassManager
 from qiskit.transpiler.layout import Layout
-from qiskit.transpiler.basepasses import AnalysisPass
+from qiskit.transpiler.basepasses import TransformationPass
 from qiskit.transpiler.exceptions import TranspilerError
+from qiskit._accelerate.nlayout import NLayout
+from qiskit._accelerate.sabre_layout import sabre_layout_and_routing
+from qiskit._accelerate.sabre_swap import (
+    Heuristic,
+    NeighborTable,
+)
+from qiskit.transpiler.passes.routing.sabre_swap import process_swaps, apply_gate
+from qiskit.tools.parallel import CPU_COUNT
+from qiskit.circuit.controlflow import ControlFlowOp, ForLoopOp
+from qiskit.converters import circuit_to_dag
 
 logger = logging.getLogger(__name__)
 
 
-class SabreLayout(AnalysisPass):
+class SabreLayout(TransformationPass):
     """Choose a Layout via iterative bidirectional routing of the input circuit.
 
-    Starting with a random initial `Layout`, the algorithm does a full routing
-    of the circuit (via the `routing_pass` method) to end up with a
-    `final_layout`. This final_layout is then used as the initial_layout for
-    routing the reverse circuit. The algorithm iterates a number of times until
-    it finds an initial_layout that reduces full routing cost.
+    The algorithm does a full routing of the circuit (via the `routing_pass`
+    method) to end up with a `final_layout`. This final_layout is then used as
+    the initial_layout for routing the reverse circuit. The algorithm iterates a
+    number of times until it finds an initial_layout that reduces full routing cost.
+
+    Prior to running the SABRE algorithm this transpiler pass will try to find the layout
+    for deepest layer that is has an isomorphic subgraph in the coupling graph. This is
+    done by progressively using the algorithm from :class:`~.VF2Layout` on the circuit
+    until a mapping is not found. This partial layout is then used to seed the SABRE algorithm
+    and then random physical bits are selected for the remaining elements in the mapping.
 
     This method exploits the reversibility of quantum circuits, and tries to
     include global circuit information in the choice of initial_layout.
 
+    By default this pass will run both layout and routing and will transform the
+    circuit so that the layout is applied to the input dag (meaning that the output
+    circuit will have ancilla qubits allocated for unused qubits on the coupling map
+    and the qubits will be reordered to match the mapped physical qubits) and then
+    routing will be applied (inserting :class:`~.SwapGate`s to account for limited
+    connectivity). This is unlike most other layout passes which are :class:`~.AnalysisPass`
+    objects and just find an initial layout and set that on the property set. This is
+    done because by default the pass will run parallel seed trials with different random
+    seeds for selecting the random initial layout and then selecting the routed output
+    which results in the least number of swap gates needed.
+
+    You can use the ``routing_pass`` argument to have this pass operate as a typical
+    layout pass. When specified this will use the specified routing pass to select an
+    initial layout only and will not run multiple seed trials.
+
     **References:**
 
     [1] Li, Gushu, Yufei Ding, and Yuan Xie. "Tackling the qubit mapping problem
@@ -50,15 +85,31 @@ class SabreLayout(AnalysisPass):
     """
 
     def __init__(
-        self, coupling_map, routing_pass=None, seed=None, max_iterations=3, swap_trials=None
+        self,
+        coupling_map,
+        routing_pass=None,
+        seed=None,
+        max_iterations=3,
+        swap_trials=None,
+        layout_trials=None,
+        skip_routing=False,
+        target=None,
+        vf2_partial_layout=True,
+        vf2_call_limit=None,
+        vf2_time_limit=None,
+        vf2_max_trials=None,
     ):
         """SabreLayout initializer.
 
         Args:
             coupling_map (Coupling): directed graph representing a coupling map.
             routing_pass (BasePass): the routing pass to use while iterating.
-                This is mutually exclusive with the ``swap_trials`` argument and
-                if both are set an error will be raised.
+                If specified this pass operates as an :class:`~.AnalysisPass` and
+                will only populate the ``layout`` field in the property set and
+                the input dag is returned unmodified. This argument is mutually
+                exclusive with the ``swap_trials`` and the ``layout_trials``
+                arguments and if this is specified at the same time as either
+                argument an error will be raised.
             seed (int): seed for setting a random first trial layout.
             max_iterations (int): number of forward-backward iterations.
             swap_trials (int): The number of trials to run of
@@ -71,20 +122,69 @@ def __init__(
                 on the number of trials run. This option is mutually exclusive
                 with the ``routing_pass`` argument and an error will be raised
                 if both are used.
+            layout_trials (int): The number of random seed trials to run
+                layout with. When > 1 the trial that resuls in the output with
+                the fewest swap gates will be selected. If this is not specified
+                (and ``routing_pass`` is not set) then the number of local
+                physical CPUs will be used as the default value. This option is
+                mutually exclusive with the ``routing_pass`` argument and an error
+                will be raised if both are used.
+            skip_routing (bool): If this is set ``True`` and ``routing_pass`` is not used
+                then routing will not be applied to the output circuit.  Only the layout
+                will be returned in the property set. This is a tradeoff to run custom
+                routing with multiple layout trials, as using this option will cause
+                SabreLayout to run the routing stage internally but not use that result.
+            target (Target): A target representing the backend device to run ``SabreLayout`` on.
+                If specified it will supersede a set value for ``coupling_map``.
+            vf2_partial_layout (bool): Run vf2 partial layout
+            vf2_call_limit (int): The number of state visits to attempt in each execution of
+                VF2 to attempt to find a partial layout.
+            vf2_time_limit (float): The total time limit in seconds to run VF2 to find a partial
+                layout
+            vf2_max_trials (int): The maximum number of trials to run VF2 to find
+                a partial layout. If this is not specified the number of trials will be limited
+                based on the number of edges in the interaction graph or the coupling graph
+                (whichever is larger) if no other limits are set. If set to a value <= 0 no
+                limit on the number of trials will be set.
 
         Raises:
-            TranspilerError: If both ``routing_pass`` and ``swap_trials`` are
-                specified
+            TranspilerError: If both ``routing_pass`` and ``swap_trials`` or
+            both ``routing_pass`` and ``layout_trials`` are specified
         """
         super().__init__()
         self.coupling_map = coupling_map
-        if routing_pass is not None and swap_trials is not None:
+        self._neighbor_table = None
+        if routing_pass is not None and (swap_trials is not None or layout_trials is not None):
             raise TranspilerError("Both routing_pass and swap_trials can't be set at the same time")
         self.routing_pass = routing_pass
         self.seed = seed
         self.max_iterations = max_iterations
         self.trials = swap_trials
-        self.swap_trials = swap_trials
+        if swap_trials is None:
+            self.swap_trials = CPU_COUNT
+        else:
+            self.swap_trials = swap_trials
+        if layout_trials is None:
+            self.layout_trials = CPU_COUNT
+        else:
+            self.layout_trials = layout_trials
+        self.skip_routing = skip_routing
+        self.target = target
+        if target is not None:
+            self.coupling_map = target.build_coupling_map()
+        self.avg_error_map = None
+        self.vf2_partial_layout = vf2_partial_layout
+        self.call_limit = vf2_call_limit
+        self.time_limit = vf2_time_limit
+        self.max_trials = vf2_max_trials
+        if self.coupling_map is not None:
+            if not self.coupling_map.is_symmetric:
+                # deepcopy is needed here to avoid modifications updating
+                # shared references in passes which require directional
+                # constraints
+                self.coupling_map = copy.deepcopy(self.coupling_map)
+                self.coupling_map.make_symmetric()
+            self._neighbor_table = NeighborTable(rx.adjacency_matrix(self.coupling_map.graph))
 
     def run(self, dag):
         """Run the SabreLayout pass on `dag`.
@@ -92,6 +192,10 @@ def run(self, dag):
         Args:
             dag (DAGCircuit): DAG to find layout for.
 
+        Returns:
+           DAGCircuit: The output dag if swap mapping was run
+            (otherwise the input dag is returned unmodified).
+
         Raises:
             TranspilerError: if dag wider than self.coupling_map
         """
@@ -99,46 +203,127 @@ def run(self, dag):
             raise TranspilerError("More virtual qubits exist than physical.")
 
         # Choose a random initial_layout.
-        if self.seed is None:
-            self.seed = np.random.randint(0, np.iinfo(np.int32).max)
-        rng = np.random.default_rng(self.seed)
+        if self.routing_pass is not None:
+            if self.seed is None:
+                seed = np.random.randint(0, np.iinfo(np.int32).max)
+            else:
+                seed = self.seed
+            rng = np.random.default_rng(seed)
 
-        physical_qubits = rng.choice(self.coupling_map.size(), len(dag.qubits), replace=False)
-        physical_qubits = rng.permutation(physical_qubits)
-        initial_layout = Layout({q: dag.qubits[i] for i, q in enumerate(physical_qubits)})
+            physical_qubits = rng.choice(self.coupling_map.size(), len(dag.qubits), replace=False)
+            physical_qubits = rng.permutation(physical_qubits)
+            initial_layout = Layout({q: dag.qubits[i] for i, q in enumerate(physical_qubits)})
 
-        if self.routing_pass is None:
-            self.routing_pass = SabreSwap(
-                self.coupling_map, "decay", seed=self.seed, fake_run=True, trials=self.swap_trials
-            )
-        else:
             self.routing_pass.fake_run = True
 
-        # Do forward-backward iterations.
-        circ = dag_to_circuit(dag)
-        rev_circ = circ.reverse_ops()
-        for _ in range(self.max_iterations):
-            for _ in ("forward", "backward"):
-                pm = self._layout_and_route_passmanager(initial_layout)
-                new_circ = pm.run(circ)
-
-                # Update initial layout and reverse the unmapped circuit.
-                pass_final_layout = pm.property_set["final_layout"]
-                final_layout = self._compose_layouts(
-                    initial_layout, pass_final_layout, new_circ.qregs
-                )
-                initial_layout = final_layout
-                circ, rev_circ = rev_circ, circ
+            # Do forward-backward iterations.
+            circ = dag_to_circuit(dag)
+            rev_circ = circ.reverse_ops()
+            for _ in range(self.max_iterations):
+                for _ in ("forward", "backward"):
+                    pm = self._layout_and_route_passmanager(initial_layout)
+                    new_circ = pm.run(circ)
 
-            # Diagnostics
-            logger.info("new initial layout")
-            logger.info(initial_layout)
+                    # Update initial layout and reverse the unmapped circuit.
+                    pass_final_layout = pm.property_set["final_layout"]
+                    final_layout = self._compose_layouts(
+                        initial_layout, pass_final_layout, new_circ.qregs
+                    )
+                    initial_layout = final_layout
+                    circ, rev_circ = rev_circ, circ
 
-        for qreg in dag.qregs.values():
-            initial_layout.add_register(qreg)
+                # Diagnostics
+                logger.info("new initial layout")
+                logger.info(initial_layout)
+
+            for qreg in dag.qregs.values():
+                initial_layout.add_register(qreg)
+            self.property_set["layout"] = initial_layout
+            self.routing_pass.fake_run = False
+            return dag
+        else:
+            dist_matrix = self.coupling_map.distance_matrix
+            original_qubit_indices = {bit: index for index, bit in enumerate(dag.qubits)}
+            original_clbit_indices = {bit: index for index, bit in enumerate(dag.clbits)}
+            partial_layout = None
+            if self.vf2_partial_layout:
+                partial_layout_virtual_bits = self._vf2_partial_layout(dag).get_virtual_bits()
+                partial_layout = [partial_layout_virtual_bits.get(i, None) for i in dag.qubits]
 
-        self.property_set["layout"] = initial_layout
-        self.routing_pass.fake_run = False
+            dag_list = []
+            for node in dag.topological_op_nodes():
+                cargs = {original_clbit_indices[x] for x in node.cargs}
+                if node.op.condition is not None:
+                    for clbit in dag._bits_in_condition(node.op.condition):
+                        cargs.add(original_clbit_indices[clbit])
+
+                dag_list.append(
+                    (
+                        node._node_id,
+                        [original_qubit_indices[x] for x in node.qargs],
+                        cargs,
+                    )
+                )
+            ((initial_layout, final_layout), swap_map, gate_order) = sabre_layout_and_routing(
+                len(dag.clbits),
+                dag_list,
+                self._neighbor_table,
+                dist_matrix,
+                Heuristic.Decay,
+                self.seed,
+                self.max_iterations,
+                self.swap_trials,
+                self.layout_trials,
+                partial_layout,
+            )
+            # Apply initial layout selected.
+            # this is a pseudo-pass manager to avoid the repeated round trip between
+            # dag and circuit and just use a dag
+            original_dag = dag
+            layout_dict = {}
+            num_qubits = len(dag.qubits)
+            for k, v in initial_layout.layout_mapping():
+                if k < num_qubits:
+                    layout_dict[dag.qubits[k]] = v
+            initital_layout = Layout(layout_dict)
+            self.property_set["layout"] = initital_layout
+            # If skip_routing is set then return the layout in the property set
+            # and throwaway the extra work we did to compute the swap map
+            if self.skip_routing:
+                return dag
+            ancilla_pass = FullAncillaAllocation(self.coupling_map)
+            ancilla_pass.property_set = self.property_set
+            dag = ancilla_pass.run(dag)
+            enlarge_pass = EnlargeWithAncilla()
+            enlarge_pass.property_set = ancilla_pass.property_set
+            dag = enlarge_pass.run(dag)
+            apply_pass = ApplyLayout()
+            apply_pass.property_set = enlarge_pass.property_set
+            dag = apply_pass.run(dag)
+            # Apply sabre swap ontop of circuit with sabre layout
+            final_layout_mapping = final_layout.layout_mapping()
+            self.property_set["final_layout"] = Layout(
+                {dag.qubits[k]: v for (k, v) in final_layout_mapping}
+            )
+            mapped_dag = dag.copy_empty_like()
+            canonical_register = dag.qregs["q"]
+            qubit_indices = {bit: idx for idx, bit in enumerate(canonical_register)}
+            original_layout = NLayout.generate_trivial_layout(self.coupling_map.size())
+            for node_id in gate_order:
+                node = original_dag._multi_graph[node_id]
+                process_swaps(
+                    swap_map,
+                    node,
+                    mapped_dag,
+                    original_layout,
+                    canonical_register,
+                    False,
+                    qubit_indices,
+                )
+                apply_gate(
+                    mapped_dag, node, original_layout, canonical_register, False, layout_dict
+                )
+            return mapped_dag
 
     def _layout_and_route_passmanager(self, initial_layout):
         """Return a passmanager for a full layout and routing.
@@ -168,3 +353,175 @@ def _compose_layouts(self, initial_layout, pass_final_layout, qregs):
         qubit_map = Layout.combine_into_edge_map(initial_layout, trivial_layout)
         final_layout = {v: pass_final_layout._v2p[qubit_map[v]] for v in initial_layout._v2p}
         return Layout(final_layout)
+
+    # TODO: Migrate this to rust as part of sabre_layout.rs after
+    # https://github.com/Qiskit/rustworkx/issues/741 is implemented and released
+    def _vf2_partial_layout(self, dag):
+        """Find a partial layout using vf2 on the deepest subgraph that is isomorphic to
+        the coupling graph."""
+        im_graph_node_map = {}
+        reverse_im_graph_node_map = {}
+        im_graph = rx.PyGraph(multigraph=False)
+        logger.debug("Buidling interaction graphs")
+        largest_im_graph = None
+        best_mapping = None
+        first_mapping = None
+        if self.avg_error_map is None:
+            self.avg_error_map = vf2_utils.build_average_error_map(
+                self.target, None, self.coupling_map
+            )
+
+        cm_graph, cm_nodes = vf2_utils.shuffle_coupling_graph(self.coupling_map, self.seed, False)
+        # To avoid trying to over optimize the result by default limit the number
+        # of trials based on the size of the graphs. For circuits with simple layouts
+        # like an all 1q circuit we don't want to sit forever trying every possible
+        # mapping in the search space if no other limits are set
+        if self.max_trials is None and self.call_limit is None and self.time_limit is None:
+            im_graph_edge_count = len(im_graph.edge_list())
+            cm_graph_edge_count = len(self.coupling_map.graph.edge_list())
+            self.max_trials = max(im_graph_edge_count, cm_graph_edge_count) + 15
+
+        start_time = time.time()
+
+        # A more efficient search pattern would be to do a binary search
+        # and find, but to conserve memory and avoid a large number of
+        # unecessary graphs this searchs from the beginning and continues
+        # until there is no vf2 match
+        def _visit(dag, weight, wire_map):
+            for node in dag.topological_op_nodes():
+                nonlocal largest_im_graph
+                largest_im_graph = im_graph.copy()
+                if getattr(node.op, "_directive", False):
+                    continue
+                if isinstance(node.op, ControlFlowOp):
+                    if isinstance(node.op, ForLoopOp):
+                        inner_weight = len(node.op.params[0]) * weight
+                    else:
+                        inner_weight = weight
+                    for block in node.op.blocks:
+                        inner_wire_map = {
+                            inner: wire_map[outer] for outer, inner in zip(node.qargs, block.qubits)
+                        }
+                        _visit(circuit_to_dag(block), inner_weight, inner_wire_map)
+                    continue
+                len_args = len(node.qargs)
+                qargs = [wire_map[q] for q in node.qargs]
+                if len_args == 1:
+                    if qargs[0] not in im_graph_node_map:
+                        weights = defaultdict(int)
+                        weights[node.name] += weight
+                        im_graph_node_map[qargs[0]] = im_graph.add_node(weights)
+                        reverse_im_graph_node_map[im_graph_node_map[qargs[0]]] = qargs[0]
+                    else:
+                        im_graph[im_graph_node_map[qargs[0]]][node.op.name] += weight
+                if len_args == 2:
+                    if qargs[0] not in im_graph_node_map:
+                        im_graph_node_map[qargs[0]] = im_graph.add_node(defaultdict(int))
+                        reverse_im_graph_node_map[im_graph_node_map[qargs[0]]] = qargs[0]
+                    if qargs[1] not in im_graph_node_map:
+                        im_graph_node_map[qargs[1]] = im_graph.add_node(defaultdict(int))
+                        reverse_im_graph_node_map[im_graph_node_map[qargs[1]]] = qargs[1]
+                    edge = (im_graph_node_map[qargs[0]], im_graph_node_map[qargs[1]])
+                    if im_graph.has_edge(*edge):
+                        im_graph.get_edge_data(*edge)[node.name] += weight
+                    else:
+                        weights = defaultdict(int)
+                        weights[node.name] += weight
+                        im_graph.add_edge(*edge, weights)
+                if len_args > 2:
+                    raise TranspilerError(
+                        "Encountered an instruction operating on more than 2 qubits, this pass "
+                        "only functions with 1 or 2 qubit operations."
+                    )
+                vf2_mapping = rx.vf2_mapping(
+                    cm_graph,
+                    im_graph,
+                    subgraph=True,
+                    id_order=False,
+                    induced=False,
+                    call_limit=self.call_limit,
+                )
+                try:
+                    nonlocal first_mapping
+                    first_mapping = next(vf2_mapping)
+                except StopIteration:
+                    break
+                nonlocal best_mapping
+                best_mapping = vf2_mapping
+                elapsed_time = time.time() - start_time
+                if (
+                    self.time_limit is not None
+                    and best_mapping is not None
+                    and elapsed_time >= self.time_limit
+                ):
+                    logger.debug(
+                        "SabreLayout VF2 heuristic has taken %s which exceeds configured max time: %s",
+                        elapsed_time,
+                        self.time_limit,
+                    )
+                    break
+
+        _visit(dag, 1, {bit: bit for bit in dag.qubits})
+        logger.debug("Finding best mappings of largest partial subgraph")
+        im_graph = largest_im_graph
+
+        def mapping_to_layout(layout_mapping):
+            return Layout({reverse_im_graph_node_map[k]: v for k, v in layout_mapping.items()})
+
+        layout_mapping = {im_i: cm_nodes[cm_i] for cm_i, im_i in first_mapping.items()}
+        chosen_layout = mapping_to_layout(layout_mapping)
+        chosen_layout_score = vf2_utils.score_layout(
+            self.avg_error_map,
+            layout_mapping,
+            im_graph_node_map,
+            reverse_im_graph_node_map,
+            im_graph,
+            False,
+        )
+        trials = 1
+        for mapping in best_mapping:  # pylint: disable=not-an-iterable
+            trials += 1
+            logger.debug("Running trial: %s", trials)
+            layout_mapping = {im_i: cm_nodes[cm_i] for cm_i, im_i in mapping.items()}
+            # If the graphs have the same number of nodes we don't need to score or do multiple
+            # trials as the score heuristic currently doesn't weigh nodes based on gates on a
+            # qubit so the scores will always all be the same
+            if len(cm_graph) == len(im_graph):
+                break
+            layout_score = vf2_utils.score_layout(
+                self.avg_error_map,
+                layout_mapping,
+                im_graph_node_map,
+                reverse_im_graph_node_map,
+                im_graph,
+                False,
+            )
+            logger.debug("Trial %s has score %s", trials, layout_score)
+            if chosen_layout is None:
+                chosen_layout = mapping_to_layout(layout_mapping)
+                chosen_layout_score = layout_score
+            elif layout_score < chosen_layout_score:
+                layout = mapping_to_layout(layout_mapping)
+                logger.debug(
+                    "Found layout %s has a lower score (%s) than previous best %s (%s)",
+                    layout,
+                    layout_score,
+                    chosen_layout,
+                    chosen_layout_score,
+                )
+                chosen_layout = layout
+                chosen_layout_score = layout_score
+            if self.max_trials and trials >= self.max_trials:
+                logger.debug("Trial %s is >= configured max trials %s", trials, self.max_trials)
+                break
+            elapsed_time = time.time() - start_time
+            if self.time_limit is not None and elapsed_time >= self.time_limit:
+                logger.debug(
+                    "VF2Layout has taken %s which exceeds configured max time: %s",
+                    elapsed_time,
+                    self.time_limit,
+                )
+                break
+        for reg in dag.qregs.values():
+            chosen_layout.add_register(reg)
+        return chosen_layout
diff --git a/qiskit/transpiler/passes/layout/vf2_utils.py b/qiskit/transpiler/passes/layout/vf2_utils.py
index 1cb868d1d247..80a09338355b 100644
--- a/qiskit/transpiler/passes/layout/vf2_utils.py
+++ b/qiskit/transpiler/passes/layout/vf2_utils.py
@@ -98,8 +98,9 @@ def score_layout(
         size = 0
     nlayout = NLayout(layout_mapping, size + 1, size + 1)
     bit_list = np.zeros(len(im_graph), dtype=np.int32)
-    for node_index in bit_map.values():
-        bit_list[node_index] = sum(im_graph[node_index].values())
+    for qubit, node_index in bit_map.items():
+        if qubit in layout_mapping:
+            bit_list[node_index] = sum(im_graph[node_index].values())
     edge_list = {
         (edge[0], edge[1]): sum(edge[2].values()) for edge in im_graph.edge_index_map().values()
     }
diff --git a/qiskit/transpiler/passes/routing/sabre_swap.py b/qiskit/transpiler/passes/routing/sabre_swap.py
index 374967fec01a..73c711d5477a 100644
--- a/qiskit/transpiler/passes/routing/sabre_swap.py
+++ b/qiskit/transpiler/passes/routing/sabre_swap.py
@@ -142,6 +142,9 @@ def __init__(self, coupling_map, heuristic="basic", seed=None, fake_run=False, t
         if coupling_map is None or coupling_map.is_symmetric:
             self.coupling_map = coupling_map
         else:
+            # A deepcopy is needed here to avoid modifications updating
+            # shared references in passes which require directional
+            # constraints
             self.coupling_map = deepcopy(coupling_map)
             self.coupling_map.make_symmetric()
         self._neighbor_table = None
@@ -237,33 +240,63 @@ def run(self, dag):
         if not self.fake_run:
             for node_id in gate_order:
                 node = dag._multi_graph[node_id]
-                self._process_swaps(swap_map, node, mapped_dag, original_layout, canonical_register)
-                self._apply_gate(mapped_dag, node, original_layout, canonical_register)
-            return mapped_dag
-        return dag
-
-    def _process_swaps(self, swap_map, node, mapped_dag, current_layout, canonical_register):
-        if node._node_id in swap_map:
-            for swap in swap_map[node._node_id]:
-                swap_qargs = [canonical_register[swap[0]], canonical_register[swap[1]]]
-                self._apply_gate(
+                process_swaps(
+                    swap_map,
+                    node,
                     mapped_dag,
-                    DAGOpNode(op=SwapGate(), qargs=swap_qargs),
-                    current_layout,
+                    original_layout,
                     canonical_register,
+                    self.fake_run,
+                    self._qubit_indices,
                 )
-                current_layout.swap_logical(*swap)
-
-    def _apply_gate(self, mapped_dag, node, current_layout, canonical_register):
-        new_node = self._transform_gate_for_layout(node, current_layout, canonical_register)
-        if self.fake_run:
-            return new_node
-        return mapped_dag.apply_operation_back(new_node.op, new_node.qargs, new_node.cargs)
-
-    def _transform_gate_for_layout(self, op_node, layout, device_qreg):
-        """Return node implementing a virtual op on given layout."""
-        mapped_op_node = copy(op_node)
-        mapped_op_node.qargs = tuple(
-            device_qreg[layout.logical_to_physical(self._qubit_indices[x])] for x in op_node.qargs
-        )
-        return mapped_op_node
+                apply_gate(
+                    mapped_dag,
+                    node,
+                    original_layout,
+                    canonical_register,
+                    self.fake_run,
+                    self._qubit_indices,
+                )
+            return mapped_dag
+        return dag
+
+
+def process_swaps(
+    swap_map,
+    node,
+    mapped_dag,
+    current_layout,
+    canonical_register,
+    fake_run,
+    qubit_indices,
+):
+    """Process swaps from SwapMap."""
+    if node._node_id in swap_map:
+        for swap in swap_map[node._node_id]:
+            swap_qargs = [canonical_register[swap[0]], canonical_register[swap[1]]]
+            apply_gate(
+                mapped_dag,
+                DAGOpNode(op=SwapGate(), qargs=swap_qargs),
+                current_layout,
+                canonical_register,
+                fake_run,
+                qubit_indices,
+            )
+            current_layout.swap_logical(*swap)
+
+
+def apply_gate(mapped_dag, node, current_layout, canonical_register, fake_run, qubit_indices):
+    """Apply gate given the current layout."""
+    new_node = transform_gate_for_layout(node, current_layout, canonical_register, qubit_indices)
+    if fake_run:
+        return new_node
+    return mapped_dag.apply_operation_back(new_node.op, new_node.qargs, new_node.cargs)
+
+
+def transform_gate_for_layout(op_node, layout, device_qreg, qubit_indices):
+    """Return node implementing a virtual op on given layout."""
+    mapped_op_node = copy(op_node)
+    mapped_op_node.qargs = tuple(
+        device_qreg[layout.logical_to_physical(qubit_indices[x])] for x in op_node.qargs
+    )
+    return mapped_op_node
diff --git a/qiskit/transpiler/preset_passmanagers/level0.py b/qiskit/transpiler/preset_passmanagers/level0.py
index 78cffd88f24a..e74e91307753 100644
--- a/qiskit/transpiler/preset_passmanagers/level0.py
+++ b/qiskit/transpiler/preset_passmanagers/level0.py
@@ -84,8 +84,15 @@ def _choose_layout_condition(property_set):
     elif layout_method == "noise_adaptive":
         _choose_layout = NoiseAdaptiveLayout(backend_properties)
     elif layout_method == "sabre":
+        skip_routing = pass_manager_config.routing_method is not None and routing_method != "sabre"
         _choose_layout = SabreLayout(
-            coupling_map, max_iterations=1, seed=seed_transpiler, swap_trials=5
+            coupling_map,
+            max_iterations=1,
+            seed=seed_transpiler,
+            swap_trials=5,
+            layout_trials=5,
+            skip_routing=skip_routing,
+            target=target,
         )
 
     # Choose routing pass
@@ -109,10 +116,17 @@ def _choose_layout_condition(property_set):
                 "layout", layout_method, pass_manager_config, optimization_level=0
             )
         else:
+
+            def _swap_mapped(property_set):
+                return property_set["final_layout"] is None
+
             layout = PassManager()
             layout.append(_given_layout)
             layout.append(_choose_layout, condition=_choose_layout_condition)
-            layout += common.generate_embed_passmanager(coupling_map)
+            embed = common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [pass_ for x in embed.passes() for pass_ in x["passes"]], condition=_swap_mapped
+            )
         routing = routing_pm
     else:
         layout = None
diff --git a/qiskit/transpiler/preset_passmanagers/level1.py b/qiskit/transpiler/preset_passmanagers/level1.py
index e3b1cdc88735..f98b7b3dcdde 100644
--- a/qiskit/transpiler/preset_passmanagers/level1.py
+++ b/qiskit/transpiler/preset_passmanagers/level1.py
@@ -34,6 +34,7 @@
 from qiskit.transpiler.passes import Optimize1qGatesDecomposition
 from qiskit.transpiler.passes import CheckMap
 from qiskit.transpiler.passes import GatesInBasis
+from qiskit.transpiler.passes import BarrierBeforeFinalMeasurements
 from qiskit.transpiler.preset_passmanagers import common
 from qiskit.transpiler.passes.layout.vf2_layout import VF2LayoutStopReason
 
@@ -138,12 +139,29 @@ def _vf2_match_not_found(property_set):
         _improve_layout = NoiseAdaptiveLayout(backend_properties)
     elif layout_method == "sabre":
         _improve_layout = SabreLayout(
-            coupling_map, max_iterations=2, seed=seed_transpiler, swap_trials=5
+            coupling_map,
+            max_iterations=2,
+            seed=seed_transpiler,
+            swap_trials=5,
+            layout_trials=5,
+            skip_routing=pass_manager_config.routing_method is not None
+            and routing_method != "sabre",
+            target=target,
+            vf2_call_limit=int(5e4),
         )
     elif layout_method is None:
         _improve_layout = common.if_has_control_flow_else(
             DenseLayout(coupling_map, backend_properties, target=target),
-            SabreLayout(coupling_map, max_iterations=2, seed=seed_transpiler, swap_trials=5),
+            SabreLayout(
+                coupling_map,
+                max_iterations=2,
+                seed=seed_transpiler,
+                swap_trials=5,
+                layout_trials=5,
+                skip_routing=pass_manager_config.routing_method is not None
+                and routing_method != "sabre",
+                target=target,
+            ),
         ).to_flow_controller()
 
     # Choose routing pass
@@ -196,12 +214,21 @@ def _opt_control(property_set):
                 "layout", layout_method, pass_manager_config, optimization_level=1
             )
         else:
+
+            def _swap_mapped(property_set):
+                return property_set["final_layout"] is None
+
             layout = PassManager()
             layout.append(_given_layout)
             layout.append(_choose_layout_0, condition=_choose_layout_condition)
             layout.append(_choose_layout_1, condition=_layout_not_perfect)
-            layout.append(_improve_layout, condition=_vf2_match_not_found)
-            layout += common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [BarrierBeforeFinalMeasurements(), _improve_layout], condition=_vf2_match_not_found
+            )
+            embed = common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [pass_ for x in embed.passes() for pass_ in x["passes"]], condition=_swap_mapped
+            )
 
         routing = routing_pm
 
diff --git a/qiskit/transpiler/preset_passmanagers/level2.py b/qiskit/transpiler/preset_passmanagers/level2.py
index adbb75d86b55..4dc43db32dce 100644
--- a/qiskit/transpiler/preset_passmanagers/level2.py
+++ b/qiskit/transpiler/preset_passmanagers/level2.py
@@ -34,6 +34,7 @@
 from qiskit.transpiler.passes import Optimize1qGatesDecomposition
 from qiskit.transpiler.passes import CommutativeCancellation
 from qiskit.transpiler.passes import GatesInBasis
+from qiskit.transpiler.passes import BarrierBeforeFinalMeasurements
 from qiskit.transpiler.preset_passmanagers import common
 from qiskit.transpiler.passes.layout.vf2_layout import VF2LayoutStopReason
 
@@ -128,7 +129,15 @@ def _vf2_match_not_found(property_set):
         _choose_layout_1 = NoiseAdaptiveLayout(backend_properties)
     elif layout_method == "sabre":
         _choose_layout_1 = SabreLayout(
-            coupling_map, max_iterations=2, seed=seed_transpiler, swap_trials=10
+            coupling_map,
+            max_iterations=2,
+            seed=seed_transpiler,
+            swap_trials=10,
+            layout_trials=10,
+            skip_routing=pass_manager_config.routing_method is not None
+            and routing_method != "sabre",
+            target=target,
+            vf2_call_limit=int(5e6),
         )
 
     # Choose routing pass
@@ -165,11 +174,20 @@ def _opt_control(property_set):
                 "layout", layout_method, pass_manager_config, optimization_level=2
             )
         else:
+
+            def _swap_mapped(property_set):
+                return property_set["final_layout"] is None
+
             layout = PassManager()
             layout.append(_given_layout)
             layout.append(_choose_layout_0, condition=_choose_layout_condition)
-            layout.append(_choose_layout_1, condition=_vf2_match_not_found)
-            layout += common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [BarrierBeforeFinalMeasurements(), _choose_layout_1], condition=_vf2_match_not_found
+            )
+            embed = common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [pass_ for x in embed.passes() for pass_ in x["passes"]], condition=_swap_mapped
+            )
         routing = routing_pm
     else:
         layout = None
diff --git a/qiskit/transpiler/preset_passmanagers/level3.py b/qiskit/transpiler/preset_passmanagers/level3.py
index 5abb44fa1c81..5f786c5d69ca 100644
--- a/qiskit/transpiler/preset_passmanagers/level3.py
+++ b/qiskit/transpiler/preset_passmanagers/level3.py
@@ -40,6 +40,7 @@
 from qiskit.transpiler.passes import ConsolidateBlocks
 from qiskit.transpiler.passes import UnitarySynthesis
 from qiskit.transpiler.passes import GatesInBasis
+from qiskit.transpiler.passes import BarrierBeforeFinalMeasurements
 from qiskit.transpiler.runningpassmanager import ConditionalController
 from qiskit.transpiler.preset_passmanagers import common
 from qiskit.transpiler.passes.layout.vf2_layout import VF2LayoutStopReason
@@ -134,7 +135,15 @@ def _vf2_match_not_found(property_set):
         _choose_layout_1 = NoiseAdaptiveLayout(backend_properties)
     elif layout_method == "sabre":
         _choose_layout_1 = SabreLayout(
-            coupling_map, max_iterations=4, seed=seed_transpiler, swap_trials=20
+            coupling_map,
+            max_iterations=4,
+            seed=seed_transpiler,
+            swap_trials=20,
+            layout_trials=20,
+            skip_routing=pass_manager_config.routing_method is not None
+            and routing_method != "sabre",
+            target=target,
+            vf2_call_limit=int(3e7),
         )
 
     # Choose routing pass
@@ -192,11 +201,20 @@ def _opt_control(property_set):
                 "layout", layout_method, pass_manager_config, optimization_level=3
             )
         else:
+
+            def _swap_mapped(property_set):
+                return property_set["final_layout"] is None
+
             layout = PassManager()
             layout.append(_given_layout)
             layout.append(_choose_layout_0, condition=_choose_layout_condition)
-            layout.append(_choose_layout_1, condition=_vf2_match_not_found)
-            layout += common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [BarrierBeforeFinalMeasurements(), _choose_layout_1], condition=_vf2_match_not_found
+            )
+            embed = common.generate_embed_passmanager(coupling_map)
+            layout.append(
+                [pass_ for x in embed.passes() for pass_ in x["passes"]], condition=_swap_mapped
+            )
         routing = routing_pm
     else:
         layout = None
diff --git a/releasenotes/notes/rusty-sabre-layout-2e1ca05d1902dcb5.yaml b/releasenotes/notes/rusty-sabre-layout-2e1ca05d1902dcb5.yaml
new file mode 100644
index 000000000000..5a835d1e00bb
--- /dev/null
+++ b/releasenotes/notes/rusty-sabre-layout-2e1ca05d1902dcb5.yaml
@@ -0,0 +1,102 @@
+---
+features:
+  - |
+    The :class:`~.SabreLayout` transpiler pass has greatly improved performance
+    as it has been re-written in Rust. As part of this rewrite the pass has been
+    transformed from an analysis pass to a transformation pass that will run both
+    layout and routing. This was done to not only improve the runtime performance
+    but also improve the quality of the results. The previous functionality of the
+    pass as an analysis pass can be retained by manually setting the ``routing_pass``
+    argument or using the new ``skip_routing`` argument.
+  - |
+    The :class:`~.SabreLayout` transpiler pass has a new constructor argument
+    ``layout_trials``. This argument is used to control how many random number
+    generator seeds will be attempted  to run :class:`~.SabreLayout` with. When
+    set the SABRE layout algorithm is run ``layout_trials`` number of times and
+    the best quality output (measured in the lowest number of swap gates added)
+    is selected. These seed trials are executed in parallel using multithreading
+    to minimize the potential performance overhead of running layout multiple
+    times. By default if this is not specified the :class:`~.SabreLayout`
+    pass will default to using the number of physical CPUs are available on the
+    local system.
+upgrade:
+  - |
+    The default behavior of the :class:`~.SabreLayout` compiler pass has
+    changed. The pass is no longer an :class:`~.AnalysisPass` and by default
+    will compute the initital layout, apply it to the circuit, and will
+    also run :class:`~.SabreSwap` internally and apply the swap mapping
+    and set the ``final_layout`` property set with the permutation caused
+    by swap insertions. This means for users running :class:`~.SabreLayout`
+    as part of a custom :class:`~.PassManager` will need to adjust the pass
+    manager to account for this (unless they were setting the ``routing_pass``
+    argument for :class:`~.SabreLayout`). This change was made in the interest
+    of improving the quality output, the layout and routing quality are highly
+    coupled and :class:`~.SabreLayout` will now run multiple parallel seed
+    trials and to calculate which seed provides the best results it needs to
+    perform both the layout and routing together. There are three ways you can
+    adjust the usage in your custom pass manager. The first is to avoid using
+    embedding in your preset pass manager. If you were previously running something
+    like::
+
+        from qiskit.transpiler import PassManager
+        from qiskit.transpiler.preset_passmanagers import common
+        from qiskit.transpiler.passes.SabreLayout
+
+        pm = PassManager()
+        pm.append(SabreLayout(coupling_map)
+        pm += common.generate_embed_passmanager(coupling_map)
+
+    to compute the layout and then apply it (which was typically followed by routing)
+    you can adjust the usage to just simply be::
+
+
+        from qiskit.transpiler import PassManager
+        from qiskit.transpiler.preset_passmanagers import common
+        from qiskit.transpiler.passes.SabreLayout
+
+        pm = PassManager()
+        pm.append(SabreLayout(coupling_map)
+
+    as :class:`~.SabreLayout` will apply the layout and you no longer need the embedding
+    stage. Alternatively, you can specify the ``routing_pass`` argument which will revert
+    :class:`~.SabreLayout` to its previous behavior. For example, if you want to run
+    :class:`~.SabreLayout` as it was run in previous releases you can do something like::
+
+        from qiskit.transpiler.passes import SabreSwap, SabreLayout
+        routing_pass = SabreSwap(
+            coupling_map, "decay", seed=seed, fake_run=True
+        )
+        layout_pass = SabreLayout(coupling_map, routing_pass=routing_pass, seed=seed)
+
+    which will have :class:`~.SabreLayout` run as an analysis pass and just set
+    the ``layout`` property set. The final approach is to leverage the ``skip_routing``
+    argument on :class:`~SabreLayout`, when this argument is set to ``True`` it will
+    skip applying the found layout and inserting the swap gates from routing. However,
+    doing this has a runtime penalty as :class:`~SabreLayout` will still be computing
+    the routing and just does not use this data. The first two approaches outlined do
+    not have additional overhead associated with them.
+  - |
+    The layouts computed by the :class:`~.SabreLayout` pass (when run without
+    the ``routing_pass`` argument) with a fixed seed value may change from
+    previous releases. This is caused by a new random number generator being
+    used as part of the rewrite of the :class:`~.SabreLayout` pass in Rust which
+    significantly improved the performance. If you rely on having consistent
+    output you can run the pass in an earlier version of Qiskit and leverage
+    :mod:`qiskit.qpy` to save the circuit and then load it using the current
+    version. Alternatively you can explicitly set the ``routing_pass`` argument
+    to an instance of :class:`~.SabreSwap` to mirror the previous behavior
+    of :class:`~.SabreLayout`::
+
+        from qiskit.transpiler.passes import SabreSwap, SabreLayout
+
+
+        routing_pass = SabreSwap(
+            coupling_map, "decay", seed=seed, fake_run=True
+        )
+        layout_pass = SabreLayout(coupling_map, routing_pass=routing_pass, seed=seed)
+
+    which will mirror the behavior of the pass in the previous release. Note, that if you
+    were using the ``swap_trials`` argument on :class:`~.SabreLayout` in previous releases
+    when adjusting the usage to this form that you will need to set ``trials`` argument
+    on the :class:`~.SabreSwap` constructor if you want to retain the previous output with
+    a fixed seed.
diff --git a/src/lib.rs b/src/lib.rs
index 3f5c5dbea612..f560fc96a452 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -23,6 +23,7 @@ mod nlayout;
 mod optimize_1q_gates;
 mod pauli_exp_val;
 mod results;
+mod sabre_layout;
 mod sabre_swap;
 mod sampled_exp_val;
 mod sparse_pauli_op;
@@ -54,6 +55,7 @@ fn _accelerate(_py: Python<'_>, m: &PyModule) -> PyResult<()> {
     m.add_wrapped(wrap_pymodule!(results::results))?;
     m.add_wrapped(wrap_pymodule!(optimize_1q_gates::optimize_1q_gates))?;
     m.add_wrapped(wrap_pymodule!(sampled_exp_val::sampled_exp_val))?;
+    m.add_wrapped(wrap_pymodule!(sabre_layout::sabre_layout))?;
     m.add_wrapped(wrap_pymodule!(vf2_layout::vf2_layout))?;
     Ok(())
 }
diff --git a/src/nlayout.rs b/src/nlayout.rs
index ae5fe007a7e6..70e2767aaed5 100644
--- a/src/nlayout.rs
+++ b/src/nlayout.rs
@@ -116,6 +116,14 @@ impl NLayout {
     pub fn copy(&self) -> NLayout {
         self.clone()
     }
+
+    #[staticmethod]
+    pub fn generate_trivial_layout(num_qubits: usize) -> Self {
+        NLayout {
+            logic_to_phys: (0..num_qubits).collect(),
+            phys_to_logic: (0..num_qubits).collect(),
+        }
+    }
 }
 
 #[pymodule]
diff --git a/src/sabre_layout.rs b/src/sabre_layout.rs
new file mode 100644
index 000000000000..9824b3072d3f
--- /dev/null
+++ b/src/sabre_layout.rs
@@ -0,0 +1,235 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2022
+//
+// This code is licensed under the Apache License, Version 2.0. You may
+// obtain a copy of this license in the LICENSE.txt file in the root directory
+// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+//
+// Any modifications or derivative works of this code must retain this
+// copyright notice, and modified files need to carry a notice indicating
+// that they have been altered from the originals.
+#![allow(clippy::too_many_arguments)]
+
+use hashbrown::HashSet;
+use ndarray::prelude::*;
+use numpy::IntoPyArray;
+use numpy::PyReadonlyArray2;
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+use pyo3::Python;
+use rand::prelude::*;
+use rand_pcg::Pcg64Mcg;
+use rayon::prelude::*;
+
+use crate::getenv_use_multiple_threads;
+use crate::nlayout::NLayout;
+use crate::sabre_swap::neighbor_table::NeighborTable;
+use crate::sabre_swap::sabre_dag::SabreDAG;
+use crate::sabre_swap::swap_map::SwapMap;
+use crate::sabre_swap::{build_swap_map_inner, Heuristic};
+
+#[pyfunction]
+pub fn sabre_layout_and_routing(
+    py: Python,
+    num_clbits: usize,
+    mut dag_nodes: Vec<(usize, Vec<usize>, HashSet<usize>)>,
+    neighbor_table: &NeighborTable,
+    distance_matrix: PyReadonlyArray2<f64>,
+    heuristic: &Heuristic,
+    seed: Option<u64>,
+    max_iterations: usize,
+    num_swap_trials: usize,
+    num_layout_trials: usize,
+    partial_layout: Option<Vec<Option<usize>>>,
+) -> ([NLayout; 2], SwapMap, PyObject) {
+    let run_in_parallel = getenv_use_multiple_threads();
+    let outer_rng = match seed {
+        Some(seed) => Pcg64Mcg::seed_from_u64(seed),
+        None => Pcg64Mcg::from_entropy(),
+    };
+    let seed_vec: Vec<u64> = outer_rng
+        .sample_iter(&rand::distributions::Standard)
+        .take(num_layout_trials)
+        .collect();
+    let dist = distance_matrix.as_array();
+    let result = if run_in_parallel && num_layout_trials > 1 {
+        seed_vec
+            .into_par_iter()
+            .enumerate()
+            .map(|(index, seed_trial)| {
+                (
+                    index,
+                    layout_trial(
+                        num_clbits,
+                        &mut dag_nodes.clone(),
+                        neighbor_table,
+                        &dist,
+                        heuristic,
+                        seed_trial,
+                        max_iterations,
+                        num_swap_trials,
+                        run_in_parallel,
+                        partial_layout.clone(),
+                    ),
+                )
+            })
+            .min_by_key(|(index, result)| {
+                (
+                    result.1.map.values().map(|x| x.len()).sum::<usize>(),
+                    *index,
+                )
+            })
+            .unwrap()
+            .1
+    } else {
+        seed_vec
+            .into_iter()
+            .map(|seed_trial| {
+                layout_trial(
+                    num_clbits,
+                    &mut dag_nodes,
+                    neighbor_table,
+                    &dist,
+                    heuristic,
+                    seed_trial,
+                    max_iterations,
+                    num_swap_trials,
+                    run_in_parallel,
+                    partial_layout.clone(),
+                )
+            })
+            .min_by_key(|result| result.1.map.values().map(|x| x.len()).sum::<usize>())
+            .unwrap()
+    };
+    (result.0, result.1, result.2.into_pyarray(py).into())
+}
+
+fn layout_trial(
+    num_clbits: usize,
+    dag_nodes: &mut Vec<(usize, Vec<usize>, HashSet<usize>)>,
+    neighbor_table: &NeighborTable,
+    distance_matrix: &ArrayView2<f64>,
+    heuristic: &Heuristic,
+    seed: u64,
+    max_iterations: usize,
+    num_swap_trials: usize,
+    run_swap_in_parallel: bool,
+    partial_layout: Option<Vec<Option<usize>>>,
+) -> ([NLayout; 2], SwapMap, Vec<usize>) {
+    // Pick a random initial layout and fully populate ancillas in that layout too
+    let num_physical_qubits = distance_matrix.shape()[0];
+    let mut rng = Pcg64Mcg::seed_from_u64(seed);
+    let mut physical_qubits: Vec<usize>;
+    match partial_layout {
+        Some(partial_layout_bits) => {
+            let used_bits: HashSet<usize> = partial_layout_bits
+                .iter()
+                .filter_map(|x| x.as_ref())
+                .copied()
+                .collect();
+            let mut free_bits: Vec<usize> = (0..num_physical_qubits)
+                .filter(|x| !used_bits.contains(x))
+                .collect();
+            free_bits.shuffle(&mut rng);
+            physical_qubits = partial_layout_bits
+                .iter()
+                .map(|x| match x {
+                    Some(phys) => *phys,
+                    None => free_bits.pop().unwrap(),
+                })
+                .collect();
+        }
+        None => {
+            physical_qubits = (0..num_physical_qubits).collect();
+            physical_qubits.shuffle(&mut rng);
+        }
+    };
+    let mut phys_to_logic = vec![0; num_physical_qubits];
+    physical_qubits
+        .iter()
+        .enumerate()
+        .for_each(|(logic, phys)| phys_to_logic[*phys] = logic);
+
+    let mut initial_layout = NLayout {
+        logic_to_phys: physical_qubits,
+        phys_to_logic,
+    };
+    let mut rev_dag_nodes: Vec<(usize, Vec<usize>, HashSet<usize>)> =
+        dag_nodes.iter().rev().cloned().collect();
+    for _iter in 0..max_iterations {
+        // forward and reverse
+        for _direction in 0..2 {
+            let dag = apply_layout(dag_nodes, &initial_layout, num_physical_qubits, num_clbits);
+            let mut pass_final_layout = NLayout::generate_trivial_layout(num_physical_qubits);
+            build_swap_map_inner(
+                num_physical_qubits,
+                &dag,
+                neighbor_table,
+                distance_matrix,
+                heuristic,
+                Some(seed),
+                &mut pass_final_layout,
+                num_swap_trials,
+                Some(run_swap_in_parallel),
+            );
+            let final_layout = compose_layout(&initial_layout, &pass_final_layout);
+            initial_layout = final_layout;
+            std::mem::swap(dag_nodes, &mut rev_dag_nodes);
+        }
+    }
+    let layout_dag = apply_layout(dag_nodes, &initial_layout, num_physical_qubits, num_clbits);
+    let mut final_layout = NLayout::generate_trivial_layout(num_physical_qubits);
+    let (swap_map, gate_order) = build_swap_map_inner(
+        num_physical_qubits,
+        &layout_dag,
+        neighbor_table,
+        distance_matrix,
+        heuristic,
+        Some(seed),
+        &mut final_layout,
+        num_swap_trials,
+        Some(run_swap_in_parallel),
+    );
+    ([initial_layout, final_layout], swap_map, gate_order)
+}
+
+fn apply_layout(
+    dag_nodes: &[(usize, Vec<usize>, HashSet<usize>)],
+    layout: &NLayout,
+    num_qubits: usize,
+    num_clbits: usize,
+) -> SabreDAG {
+    let layout_dag_nodes: Vec<(usize, Vec<usize>, HashSet<usize>)> = dag_nodes
+        .iter()
+        .map(|(node_index, qargs, cargs)| {
+            let new_qargs: Vec<usize> = qargs.iter().map(|n| layout.logic_to_phys[*n]).collect();
+            (*node_index, new_qargs, cargs.clone())
+        })
+        .collect();
+    SabreDAG::new(num_qubits, num_clbits, layout_dag_nodes).unwrap()
+}
+
+fn compose_layout(initial_layout: &NLayout, final_layout: &NLayout) -> NLayout {
+    let logic_to_phys = initial_layout
+        .logic_to_phys
+        .iter()
+        .map(|n| final_layout.logic_to_phys[*n])
+        .collect();
+    let phys_to_logic = final_layout
+        .phys_to_logic
+        .iter()
+        .map(|n| initial_layout.phys_to_logic[*n])
+        .collect();
+
+    NLayout {
+        logic_to_phys,
+        phys_to_logic,
+    }
+}
+
+#[pymodule]
+pub fn sabre_layout(_py: Python, m: &PyModule) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(sabre_layout_and_routing))?;
+    Ok(())
+}
diff --git a/src/sabre_swap/mod.rs b/src/sabre_swap/mod.rs
index 7ab6aff816bb..d10e719c570d 100644
--- a/src/sabre_swap/mod.rs
+++ b/src/sabre_swap/mod.rs
@@ -155,9 +155,38 @@ pub fn build_swap_map(
     seed: Option<u64>,
     layout: &mut NLayout,
     num_trials: usize,
+    run_in_parallel: Option<bool>,
 ) -> (SwapMap, PyObject) {
-    let run_in_parallel = getenv_use_multiple_threads() && num_trials > 1;
     let dist = distance_matrix.as_array();
+    let (swap_map, gate_order) = build_swap_map_inner(
+        num_qubits,
+        dag,
+        neighbor_table,
+        &dist,
+        heuristic,
+        seed,
+        layout,
+        num_trials,
+        run_in_parallel,
+    );
+    (swap_map, gate_order.into_pyarray(py).into())
+}
+
+pub fn build_swap_map_inner(
+    num_qubits: usize,
+    dag: &SabreDAG,
+    neighbor_table: &NeighborTable,
+    dist: &ArrayView2<f64>,
+    heuristic: &Heuristic,
+    seed: Option<u64>,
+    layout: &mut NLayout,
+    num_trials: usize,
+    run_in_parallel: Option<bool>,
+) -> (SwapMap, Vec<usize>) {
+    let run_in_parallel = match run_in_parallel {
+        Some(run_in_parallel) => run_in_parallel,
+        None => getenv_use_multiple_threads() && num_trials > 1,
+    };
     let coupling_graph: DiGraph<(), ()> = cmap_from_neighor_table(neighbor_table);
     let outer_rng = match seed {
         Some(seed) => Pcg64Mcg::seed_from_u64(seed),
@@ -178,7 +207,7 @@ pub fn build_swap_map(
                         num_qubits,
                         dag,
                         neighbor_table,
-                        &dist,
+                        dist,
                         &coupling_graph,
                         heuristic,
                         seed_trial,
@@ -202,7 +231,7 @@ pub fn build_swap_map(
                     num_qubits,
                     dag,
                     neighbor_table,
-                    &dist,
+                    dist,
                     &coupling_graph,
                     heuristic,
                     seed_trial,
@@ -217,7 +246,7 @@ pub fn build_swap_map(
         SwapMap {
             map: result.out_map,
         },
-        result.gate_order.into_pyarray(py).into(),
+        result.gate_order,
     )
 }
 
diff --git a/test/python/transpiler/test_preset_passmanagers.py b/test/python/transpiler/test_preset_passmanagers.py
index 148488edabf4..9184bd996583 100644
--- a/test/python/transpiler/test_preset_passmanagers.py
+++ b/test/python/transpiler/test_preset_passmanagers.py
@@ -716,32 +716,78 @@ def test_layout_tokyo_fully_connected_cx(self, level):
         }
 
         sabre_layout = {
-            11: qr[0],
-            17: qr[1],
-            16: qr[2],
-            6: qr[3],
-            18: qr[4],
             0: ancilla[0],
             1: ancilla[1],
             2: ancilla[2],
             3: ancilla[3],
             4: ancilla[4],
-            5: ancilla[5],
-            7: ancilla[6],
-            8: ancilla[7],
-            9: ancilla[8],
-            10: ancilla[9],
-            12: ancilla[10],
-            13: ancilla[11],
-            14: ancilla[12],
-            15: ancilla[13],
+            5: qr[1],
+            6: qr[2],
+            7: ancilla[5],
+            8: ancilla[6],
+            9: ancilla[7],
+            10: qr[3],
+            11: qr[0],
+            12: qr[4],
+            13: ancilla[8],
+            14: ancilla[9],
+            15: ancilla[10],
+            16: ancilla[11],
+            17: ancilla[12],
+            18: ancilla[13],
+            19: ancilla[14],
+        }
+
+        sabre_layout_lvl_2 = {
+            0: ancilla[0],
+            1: qr[4],
+            2: ancilla[1],
+            3: ancilla[2],
+            4: ancilla[3],
+            5: qr[1],
+            6: qr[0],
+            7: ancilla[4],
+            8: ancilla[5],
+            9: ancilla[6],
+            10: qr[2],
+            11: qr[3],
+            12: ancilla[7],
+            13: ancilla[8],
+            14: ancilla[9],
+            15: ancilla[10],
+            16: ancilla[11],
+            17: ancilla[12],
+            18: ancilla[13],
+            19: ancilla[14],
+        }
+
+        sabre_layout_lvl_3 = {
+            0: ancilla[0],
+            1: qr[4],
+            2: ancilla[1],
+            3: ancilla[2],
+            4: ancilla[3],
+            5: qr[1],
+            6: qr[0],
+            7: ancilla[4],
+            8: ancilla[5],
+            9: ancilla[6],
+            10: qr[2],
+            11: qr[3],
+            12: ancilla[7],
+            13: ancilla[8],
+            14: ancilla[9],
+            15: ancilla[10],
+            16: ancilla[11],
+            17: ancilla[12],
+            18: ancilla[13],
             19: ancilla[14],
         }
 
         expected_layout_level0 = trivial_layout
         expected_layout_level1 = sabre_layout
-        expected_layout_level2 = sabre_layout
-        expected_layout_level3 = sabre_layout
+        expected_layout_level2 = sabre_layout_lvl_2
+        expected_layout_level3 = sabre_layout_lvl_3
 
         expected_layouts = [
             expected_layout_level0,
@@ -910,8 +956,11 @@ def test_1(self, circuit, level):
         resulting_basis = {node.name for node in circuit_to_dag(result).op_nodes()}
         self.assertIn("swap", resulting_basis)
 
+    # Skipping optimization level 3 because the swap gates get absorbed into
+    # a unitary block as part of the KAK decompostion optimization passes and
+    # optimized away.
     @combine(
-        level=[0, 1, 2, 3],
+        level=[0, 1, 2],
         dsc="If swap in basis, do not decompose it. level: {level}",
         name="level{level}",
     )
@@ -931,7 +980,7 @@ def test_2(self, level):
             optimization_level=level,
             basis_gates=basis,
             coupling_map=coupling_map,
-            seed_transpiler=42123,
+            seed_transpiler=421234242,
         )
         self.assertIsInstance(result, QuantumCircuit)
         resulting_basis = {node.name for node in circuit_to_dag(result).op_nodes()}
diff --git a/test/python/transpiler/test_sabre_layout.py b/test/python/transpiler/test_sabre_layout.py
index 7aabc5ef1a83..f1757a955f56 100644
--- a/test/python/transpiler/test_sabre_layout.py
+++ b/test/python/transpiler/test_sabre_layout.py
@@ -55,15 +55,15 @@ def test_5q_circuit_20q_coupling(self):
         circuit.cx(qr[1], qr[2])
 
         dag = circuit_to_dag(circuit)
-        pass_ = SabreLayout(CouplingMap(self.cmap20), seed=0, swap_trials=32)
+        pass_ = SabreLayout(CouplingMap(self.cmap20), seed=0, swap_trials=32, layout_trials=32)
         pass_.run(dag)
 
         layout = pass_.property_set["layout"]
-        self.assertEqual(layout[qr[0]], 10)
-        self.assertEqual(layout[qr[1]], 12)
-        self.assertEqual(layout[qr[2]], 7)
-        self.assertEqual(layout[qr[3]], 11)
-        self.assertEqual(layout[qr[4]], 13)
+        self.assertEqual(layout[qr[0]], 18)
+        self.assertEqual(layout[qr[1]], 11)
+        self.assertEqual(layout[qr[2]], 13)
+        self.assertEqual(layout[qr[3]], 12)
+        self.assertEqual(layout[qr[4]], 14)
 
     def test_6q_circuit_20q_coupling(self):
         """Test finds layout for 6q circuit on 20q device."""
@@ -91,16 +91,16 @@ def test_6q_circuit_20q_coupling(self):
         circuit.cx(qr1[1], qr0[0])
 
         dag = circuit_to_dag(circuit)
-        pass_ = SabreLayout(CouplingMap(self.cmap20), seed=0)
+        pass_ = SabreLayout(CouplingMap(self.cmap20), seed=0, swap_trials=32, layout_trials=32)
         pass_.run(dag)
 
         layout = pass_.property_set["layout"]
-        self.assertEqual(layout[qr0[0]], 2)
-        self.assertEqual(layout[qr0[1]], 3)
-        self.assertEqual(layout[qr0[2]], 10)
-        self.assertEqual(layout[qr1[0]], 1)
-        self.assertEqual(layout[qr1[1]], 7)
-        self.assertEqual(layout[qr1[2]], 5)
+        self.assertEqual(layout[qr0[0]], 12)
+        self.assertEqual(layout[qr0[1]], 7)
+        self.assertEqual(layout[qr0[2]], 14)
+        self.assertEqual(layout[qr1[0]], 11)
+        self.assertEqual(layout[qr1[1]], 18)
+        self.assertEqual(layout[qr1[2]], 13)
 
     def test_layout_with_classical_bits(self):
         """Test sabre layout with classical bits recreate from issue #8635."""
@@ -132,14 +132,14 @@ def test_layout_with_classical_bits(self):
         res = transpile(qc, FakeKolkata(), layout_method="sabre", seed_transpiler=1234)
         self.assertIsInstance(res, QuantumCircuit)
         layout = res._layout.initial_layout
-        self.assertEqual(layout[qc.qubits[0]], 14)
-        self.assertEqual(layout[qc.qubits[1]], 19)
-        self.assertEqual(layout[qc.qubits[2]], 7)
-        self.assertEqual(layout[qc.qubits[3]], 13)
-        self.assertEqual(layout[qc.qubits[4]], 6)
-        self.assertEqual(layout[qc.qubits[5]], 16)
-        self.assertEqual(layout[qc.qubits[6]], 18)
-        self.assertEqual(layout[qc.qubits[7]], 26)
+        self.assertEqual(layout[qc.qubits[0]], 11)
+        self.assertEqual(layout[qc.qubits[1]], 22)
+        self.assertEqual(layout[qc.qubits[2]], 21)
+        self.assertEqual(layout[qc.qubits[3]], 19)
+        self.assertEqual(layout[qc.qubits[4]], 26)
+        self.assertEqual(layout[qc.qubits[5]], 8)
+        self.assertEqual(layout[qc.qubits[6]], 17)
+        self.assertEqual(layout[qc.qubits[7]], 1)
 
     # pylint: disable=line-too-long
     def test_layout_many_search_trials(self):
@@ -194,19 +194,19 @@ def test_layout_many_search_trials(self):
         self.assertIsInstance(res, QuantumCircuit)
         layout = res._layout.initial_layout
         self.assertEqual(layout[qc.qubits[0]], 19)
-        self.assertEqual(layout[qc.qubits[1]], 22)
-        self.assertEqual(layout[qc.qubits[2]], 17)
+        self.assertEqual(layout[qc.qubits[1]], 10)
+        self.assertEqual(layout[qc.qubits[2]], 1)
         self.assertEqual(layout[qc.qubits[3]], 14)
-        self.assertEqual(layout[qc.qubits[4]], 18)
-        self.assertEqual(layout[qc.qubits[5]], 9)
-        self.assertEqual(layout[qc.qubits[6]], 11)
-        self.assertEqual(layout[qc.qubits[7]], 25)
-        self.assertEqual(layout[qc.qubits[8]], 16)
-        self.assertEqual(layout[qc.qubits[9]], 3)
-        self.assertEqual(layout[qc.qubits[10]], 12)
+        self.assertEqual(layout[qc.qubits[4]], 4)
+        self.assertEqual(layout[qc.qubits[5]], 11)
+        self.assertEqual(layout[qc.qubits[6]], 8)
+        self.assertEqual(layout[qc.qubits[7]], 7)
+        self.assertEqual(layout[qc.qubits[8]], 9)
+        self.assertEqual(layout[qc.qubits[9]], 0)
+        self.assertEqual(layout[qc.qubits[10]], 5)
         self.assertEqual(layout[qc.qubits[11]], 13)
-        self.assertEqual(layout[qc.qubits[12]], 20)
-        self.assertEqual(layout[qc.qubits[13]], 8)
+        self.assertEqual(layout[qc.qubits[12]], 2)
+        self.assertEqual(layout[qc.qubits[13]], 3)
 
 
 if __name__ == "__main__":
diff --git a/test/python/transpiler/test_unitary_synthesis.py b/test/python/transpiler/test_unitary_synthesis.py
index da4e87dcd85b..89af400e63cb 100644
--- a/test/python/transpiler/test_unitary_synthesis.py
+++ b/test/python/transpiler/test_unitary_synthesis.py
@@ -42,9 +42,6 @@
     SabreLayout,
     Depth,
     FixedPoint,
-    FullAncillaAllocation,
-    EnlargeWithAncilla,
-    ApplyLayout,
     Unroll3qOrMore,
     CheckMap,
     BarrierBeforeFinalMeasurements,
@@ -486,7 +483,6 @@ def _repeat_condition(property_set):
 
             seed = 2
             _map = [SabreLayout(coupling_map, max_iterations=2, seed=seed)]
-            _embed = [FullAncillaAllocation(coupling_map), EnlargeWithAncilla(), ApplyLayout()]
             _unroll3q = Unroll3qOrMore()
             _swap_check = CheckMap(coupling_map)
             _swap = [
@@ -509,7 +505,6 @@ def _repeat_condition(property_set):
 
             pm = PassManager()
             pm.append(_map)  # map to hardware by inserting swaps
-            pm.append(_embed)
             pm.append(_unroll3q)
             pm.append(_swap_check)
             pm.append(_swap)
diff --git a/tox.ini b/tox.ini
index 59cc9f6748b6..8dbd5b2d32c2 100644
--- a/tox.ini
+++ b/tox.ini
@@ -14,7 +14,7 @@ setenv =
   QISKIT_SUPRESS_PACKAGING_WARNINGS=Y
   QISKIT_TEST_CAPTURE_STREAMS=1
   QISKIT_PARALLEL=FALSE
-passenv = RAYON_NUM_THREADS OMP_NUM_THREADS QISKIT_PARALLEL RUST_BACKTRACE SETUPTOOLS_ENABLE_FEATURES
+passenv = RAYON_NUM_THREADS OMP_NUM_THREADS QISKIT_PARALLEL RUST_BACKTRACE SETUPTOOLS_ENABLE_FEATURES QISKIT_IN_PARALLEL
 deps = setuptools_rust  # This is work around for the bug of tox 3 (see #8606 for more details.)
        -r{toxinidir}/requirements.txt
        -r{toxinidir}/requirements-dev.txt