diff --git a/halo2_proofs/Cargo.toml b/halo2_proofs/Cargo.toml
index 113e3e88..1d744357 100644
--- a/halo2_proofs/Cargo.toml
+++ b/halo2_proofs/Cargo.toml
@@ -77,6 +77,7 @@ getrandom = { version = "0.2", features = ["js"] }
 default = ["batch"]
 dev-graph = ["plotters", "tabbycat"]
 gadget-traces = ["backtrace"]
+thread-safe-region = []
 sanity-checks = []
 batch = ["rand_core/getrandom"]
 circuit-params = []
@@ -87,3 +88,5 @@ bench = false
 [[example]]
 name = "circuit-layout"
 required-features = ["dev-graph"]
+
+
diff --git a/halo2_proofs/examples/vector-mul.rs b/halo2_proofs/examples/vector-mul.rs
new file mode 100644
index 00000000..23dc6d27
--- /dev/null
+++ b/halo2_proofs/examples/vector-mul.rs
@@ -0,0 +1,350 @@
+use std::marker::PhantomData;
+
+use halo2_proofs::{
+    arithmetic::Field,
+    circuit::{AssignedCell, Chip, Layouter, Region, SimpleFloorPlanner, Value},
+    plonk::{Advice, Circuit, Column, ConstraintSystem, Error, Instance, Selector},
+    poly::Rotation,
+};
+
+// ANCHOR: instructions
+trait NumericInstructions<F: Field>: Chip<F> {
+    /// Variable representing a number.
+    type Num;
+
+    /// Loads a number into the circuit as a private input.
+    fn load_private(
+        &self,
+        layouter: impl Layouter<F>,
+        a: &[Value<F>],
+    ) -> Result<Vec<Self::Num>, Error>;
+
+    /// Returns `c = a * b`. The caller is responsible for ensuring that `a.len() == b.len()`.
+    fn mul(
+        &self,
+        layouter: impl Layouter<F>,
+        a: &[Self::Num],
+        b: &[Self::Num],
+    ) -> Result<Vec<Self::Num>, Error>;
+
+    /// Exposes a number as a public input to the circuit.
+    fn expose_public(
+        &self,
+        layouter: impl Layouter<F>,
+        num: &Self::Num,
+        row: usize,
+    ) -> Result<(), Error>;
+}
+// ANCHOR_END: instructions
+
+// ANCHOR: chip
+/// The chip that will implement our instructions! Chips store their own
+/// config, as well as type markers if necessary.
+struct FieldChip<F: Field> {
+    config: FieldConfig,
+    _marker: PhantomData<F>,
+}
+// ANCHOR_END: chip
+
+// ANCHOR: chip-config
+/// Chip state is stored in a config struct. This is generated by the chip
+/// during configuration, and then stored inside the chip.
+#[derive(Clone, Debug)]
+struct FieldConfig {
+    /// For this chip, we will use two advice columns to implement our instructions.
+    /// These are also the columns through which we communicate with other parts of
+    /// the circuit.
+    advice: [Column<Advice>; 3],
+
+    /// This is the public input (instance) column.
+    instance: Column<Instance>,
+
+    // We need a selector to enable the multiplication gate, so that we aren't placing
+    // any constraints on cells where `NumericInstructions::mul` is not being used.
+    // This is important when building larger circuits, where columns are used by
+    // multiple sets of instructions.
+    s_mul: Selector,
+}
+
+impl<F: Field> FieldChip<F> {
+    fn construct(config: <Self as Chip<F>>::Config) -> Self {
+        Self {
+            config,
+            _marker: PhantomData,
+        }
+    }
+
+    fn configure(
+        meta: &mut ConstraintSystem<F>,
+        advice: [Column<Advice>; 3],
+        instance: Column<Instance>,
+    ) -> <Self as Chip<F>>::Config {
+        meta.enable_equality(instance);
+        for column in &advice {
+            meta.enable_equality(*column);
+        }
+        let s_mul = meta.selector();
+
+        // Define our multiplication gate!
+        meta.create_gate("mul", |meta| {
+            // To implement multiplication, we need three advice cells and a selector
+            // cell. We arrange them like so:
+            //
+            // | a0  | a1  | a2  | s_mul |
+            // |-----|-----|-----|-------|
+            // | lhs | rhs | out | s_mul |
+            //
+            // Gates may refer to any relative offsets we want, but each distinct
+            // offset adds a cost to the proof. The most common offsets are 0 (the
+            // current row), 1 (the next row), and -1 (the previous row), for which
+            // `Rotation` has specific constructors.
+            let lhs = meta.query_advice(advice[0], Rotation::cur());
+            let rhs = meta.query_advice(advice[1], Rotation::cur());
+            let out = meta.query_advice(advice[2], Rotation::cur());
+            let s_mul = meta.query_selector(s_mul);
+
+            // Finally, we return the polynomial expressions that constrain this gate.
+            // For our multiplication gate, we only need a single polynomial constraint.
+            //
+            // The polynomial expressions returned from `create_gate` will be
+            // constrained by the proving system to equal zero. Our expression
+            // has the following properties:
+            // - When s_mul = 0, any value is allowed in lhs, rhs, and out.
+            // - When s_mul != 0, this constrains lhs * rhs = out.
+            vec![s_mul * (lhs * rhs - out)]
+        });
+
+        FieldConfig {
+            advice,
+            instance,
+            s_mul,
+        }
+    }
+}
+// ANCHOR_END: chip-config
+
+// ANCHOR: chip-impl
+impl<F: Field> Chip<F> for FieldChip<F> {
+    type Config = FieldConfig;
+    type Loaded = ();
+
+    fn config(&self) -> &Self::Config {
+        &self.config
+    }
+
+    fn loaded(&self) -> &Self::Loaded {
+        &()
+    }
+}
+// ANCHOR_END: chip-impl
+
+// ANCHOR: instructions-impl
+/// A variable representing a number.
+#[derive(Clone, Debug)]
+struct Number<F: Field>(AssignedCell<F, F>);
+
+impl<F: Field> NumericInstructions<F> for FieldChip<F> {
+    type Num = Number<F>;
+
+    fn load_private(
+        &self,
+        mut layouter: impl Layouter<F>,
+        values: &[Value<F>],
+    ) -> Result<Vec<Self::Num>, Error> {
+        let config = self.config();
+
+        layouter.assign_region(
+            || "load private",
+            |mut region| {
+                values
+                    .iter()
+                    .enumerate()
+                    .map(|(i, value)| {
+                        region
+                            .assign_advice(|| "private input", config.advice[0], i, || *value)
+                            .map(Number)
+                    })
+                    .collect()
+            },
+        )
+    }
+
+    fn mul(
+        &self,
+        mut layouter: impl Layouter<F>,
+        a: &[Self::Num],
+        b: &[Self::Num],
+    ) -> Result<Vec<Self::Num>, Error> {
+        let config = self.config();
+        assert_eq!(a.len(), b.len());
+
+        #[cfg(feature = "thread-safe-region")]
+        {
+            use rayon::prelude::{
+                IndexedParallelIterator, IntoParallelRefIterator, ParallelIterator,
+            };
+            layouter.assign_region(
+                || "mul",
+                |region: Region<'_, F>| {
+                    let thread_safe_region = std::sync::Mutex::new(region);
+                    a.par_iter()
+                        .zip(b.par_iter())
+                        .enumerate()
+                        .map(|(i, (a, b))| {
+                            let mut region = thread_safe_region.lock().unwrap();
+
+                            config.s_mul.enable(&mut region, i)?;
+
+                            a.0.copy_advice(|| "lhs", &mut region, config.advice[0], i)?;
+                            b.0.copy_advice(|| "rhs", &mut region, config.advice[1], i)?;
+
+                            let value = a.0.value().copied() * b.0.value();
+
+                            // Finally, we do the assignment to the output, returning a
+                            // variable to be used in another part of the circuit.
+                            region
+                                .assign_advice(|| "lhs * rhs", config.advice[2], i, || value)
+                                .map(Number)
+                        })
+                        .collect()
+                },
+            )
+        }
+
+        #[cfg(not(feature = "thread-safe-region"))]
+        layouter.assign_region(
+            || "mul",
+            |mut region: Region<'_, F>| {
+                a.iter()
+                    .zip(b.iter())
+                    .enumerate()
+                    .map(|(i, (a, b))| {
+                        config.s_mul.enable(&mut region, i)?;
+
+                        a.0.copy_advice(|| "lhs", &mut region, config.advice[0], i)?;
+                        b.0.copy_advice(|| "rhs", &mut region, config.advice[1], i)?;
+
+                        let value = a.0.value().copied() * b.0.value();
+
+                        // Finally, we do the assignment to the output, returning a
+                        // variable to be used in another part of the circuit.
+                        region
+                            .assign_advice(|| "lhs * rhs", config.advice[2], i, || value)
+                            .map(Number)
+                    })
+                    .collect()
+            },
+        )
+    }
+
+    fn expose_public(
+        &self,
+        mut layouter: impl Layouter<F>,
+        num: &Self::Num,
+        row: usize,
+    ) -> Result<(), Error> {
+        let config = self.config();
+
+        layouter.constrain_instance(num.0.cell(), config.instance, row)
+    }
+}
+// ANCHOR_END: instructions-impl
+
+// ANCHOR: circuit
+/// The full circuit implementation.
+///
+/// In this struct we store the private input variables. We use `Option<F>` because
+/// they won't have any value during key generation. During proving, if any of these
+/// were `None` we would get an error.
+#[derive(Default)]
+struct MyCircuit<F: Field> {
+    a: Vec<Value<F>>,
+    b: Vec<Value<F>>,
+}
+
+impl<F: Field> Circuit<F> for MyCircuit<F> {
+    // Since we are using a single chip for everything, we can just reuse its config.
+    type Config = FieldConfig;
+    type FloorPlanner = SimpleFloorPlanner;
+    #[cfg(feature = "circuit-params")]
+    type Params = ();
+
+    fn without_witnesses(&self) -> Self {
+        Self::default()
+    }
+
+    fn configure(meta: &mut ConstraintSystem<F>) -> Self::Config {
+        // We create the three advice columns that FieldChip uses for I/O.
+        let advice = [
+            meta.advice_column(),
+            meta.advice_column(),
+            meta.advice_column(),
+        ];
+
+        // We also need an instance column to store public inputs.
+        let instance = meta.instance_column();
+
+        FieldChip::configure(meta, advice, instance)
+    }
+
+    fn synthesize(
+        &self,
+        config: Self::Config,
+        mut layouter: impl Layouter<F>,
+    ) -> Result<(), Error> {
+        let field_chip = FieldChip::<F>::construct(config);
+
+        // Load our private values into the circuit.
+        let a = field_chip.load_private(layouter.namespace(|| "load a"), &self.a)?;
+        let b = field_chip.load_private(layouter.namespace(|| "load b"), &self.b)?;
+
+        let ab = field_chip.mul(layouter.namespace(|| "a * b"), &a, &b)?;
+
+        for (i, c) in ab.iter().enumerate() {
+            // Expose the result as a public input to the circuit.
+            field_chip.expose_public(layouter.namespace(|| "expose c"), c, i)?;
+        }
+        Ok(())
+    }
+}
+// ANCHOR_END: circuit
+
+fn main() {
+    use halo2_proofs::dev::MockProver;
+    use halo2curves::pasta::Fp;
+
+    const N: usize = 20000;
+    // ANCHOR: test-circuit
+    // The number of rows in our circuit cannot exceed 2^k. Since our example
+    // circuit is very small, we can pick a very small value here.
+    let k = 16;
+
+    // Prepare the private and public inputs to the circuit!
+    let a = [Fp::from(2); N];
+    let b = [Fp::from(3); N];
+    let c: Vec<Fp> = a.iter().zip(b).map(|(&a, b)| a * b).collect();
+
+    // Instantiate the circuit with the private inputs.
+    let circuit = MyCircuit {
+        a: a.iter().map(|&x| Value::known(x)).collect(),
+        b: b.iter().map(|&x| Value::known(x)).collect(),
+    };
+
+    // Arrange the public input. We expose the multiplication result in row 0
+    // of the instance column, so we position it there in our public inputs.
+    let mut public_inputs = c;
+
+    let start = std::time::Instant::now();
+    // Given the correct public input, our circuit will verify.
+    let prover = MockProver::run(k, &circuit, vec![public_inputs.clone()]).unwrap();
+    assert_eq!(prover.verify(), Ok(()));
+    println!("positive test took {:?}", start.elapsed());
+
+    // If we try some other public input, the proof will fail!
+    let start = std::time::Instant::now();
+    public_inputs[0] += Fp::one();
+    let prover = MockProver::run(k, &circuit, vec![public_inputs]).unwrap();
+    assert!(prover.verify().is_err());
+    println!("negative test took {:?}", start.elapsed());
+    // ANCHOR_END: test-circuit
+}
diff --git a/halo2_proofs/src/circuit/floor_planner/single_pass.rs b/halo2_proofs/src/circuit/floor_planner/single_pass.rs
index 6be22977..54268006 100644
--- a/halo2_proofs/src/circuit/floor_planner/single_pass.rs
+++ b/halo2_proofs/src/circuit/floor_planner/single_pass.rs
@@ -7,7 +7,7 @@ use ff::Field;
 
 use crate::{
     circuit::{
-        layouter::{RegionColumn, RegionLayouter, RegionShape, TableLayouter},
+        layouter::{RegionColumn, RegionLayouter, RegionShape, SyncDeps, TableLayouter},
         Cell, Layouter, Region, RegionIndex, RegionStart, Table, Value,
     },
     plonk::{
@@ -25,7 +25,7 @@ use crate::{
 pub struct SimpleFloorPlanner;
 
 impl FloorPlanner for SimpleFloorPlanner {
-    fn synthesize<F: Field, CS: Assignment<F>, C: Circuit<F>>(
+    fn synthesize<F: Field, CS: Assignment<F> + SyncDeps, C: Circuit<F>>(
         cs: &mut CS,
         circuit: &C,
         config: C::Config,
@@ -73,7 +73,9 @@ impl<'a, F: Field, CS: Assignment<F>> SingleChipLayouter<'a, F, CS> {
     }
 }
 
-impl<'a, F: Field, CS: Assignment<F> + 'a> Layouter<F> for SingleChipLayouter<'a, F, CS> {
+impl<'a, F: Field, CS: Assignment<F> + 'a + SyncDeps> Layouter<F>
+    for SingleChipLayouter<'a, F, CS>
+{
     type Root = Self;
 
     fn assign_region<A, AR, N, NR>(&mut self, name: N, mut assignment: A) -> Result<AR, Error>
@@ -263,7 +265,12 @@ impl<'r, 'a, F: Field, CS: Assignment<F> + 'a> SingleChipLayouterRegion<'r, 'a,
     }
 }
 
-impl<'r, 'a, F: Field, CS: Assignment<F> + 'a> RegionLayouter<F>
+impl<'r, 'a, F: Field, CS: Assignment<F> + 'a + SyncDeps> SyncDeps
+    for SingleChipLayouterRegion<'r, 'a, F, CS>
+{
+}
+
+impl<'r, 'a, F: Field, CS: Assignment<F> + 'a + SyncDeps> RegionLayouter<F>
     for SingleChipLayouterRegion<'r, 'a, F, CS>
 {
     fn enable_selector<'v>(
diff --git a/halo2_proofs/src/circuit/floor_planner/v1.rs b/halo2_proofs/src/circuit/floor_planner/v1.rs
index 52ca2f7f..3df8c181 100644
--- a/halo2_proofs/src/circuit/floor_planner/v1.rs
+++ b/halo2_proofs/src/circuit/floor_planner/v1.rs
@@ -5,7 +5,7 @@ use ff::Field;
 use crate::{
     circuit::{
         floor_planner::single_pass::SimpleTableLayouter,
-        layouter::{RegionColumn, RegionLayouter, RegionShape, TableLayouter},
+        layouter::{RegionColumn, RegionLayouter, RegionShape, SyncDeps, TableLayouter},
         Cell, Layouter, Region, RegionIndex, RegionStart, Table, Value,
     },
     plonk::{
@@ -43,7 +43,7 @@ impl<'a, F: Field, CS: Assignment<F> + 'a> fmt::Debug for V1Plan<'a, F, CS> {
     }
 }
 
-impl<'a, F: Field, CS: Assignment<F>> V1Plan<'a, F, CS> {
+impl<'a, F: Field, CS: Assignment<F> + SyncDeps> V1Plan<'a, F, CS> {
     /// Creates a new v1 layouter.
     pub fn new(cs: &'a mut CS) -> Result<Self, Error> {
         let ret = V1Plan {
@@ -57,7 +57,7 @@ impl<'a, F: Field, CS: Assignment<F>> V1Plan<'a, F, CS> {
 }
 
 impl FloorPlanner for V1 {
-    fn synthesize<F: Field, CS: Assignment<F>, C: Circuit<F>>(
+    fn synthesize<F: Field, CS: Assignment<F> + SyncDeps, C: Circuit<F>>(
         cs: &mut CS,
         circuit: &C,
         config: C::Config,
@@ -160,7 +160,7 @@ impl<'p, 'a, F: Field, CS: Assignment<F> + 'a> V1Pass<'p, 'a, F, CS> {
     }
 }
 
-impl<'p, 'a, F: Field, CS: Assignment<F> + 'a> Layouter<F> for V1Pass<'p, 'a, F, CS> {
+impl<'p, 'a, F: Field, CS: Assignment<F> + SyncDeps> Layouter<F> for V1Pass<'p, 'a, F, CS> {
     type Root = Self;
 
     fn assign_region<A, AR, N, NR>(&mut self, name: N, assignment: A) -> Result<AR, Error>
@@ -264,7 +264,7 @@ pub struct AssignmentPass<'p, 'a, F: Field, CS: Assignment<F> + 'a> {
     region_index: usize,
 }
 
-impl<'p, 'a, F: Field, CS: Assignment<F> + 'a> AssignmentPass<'p, 'a, F, CS> {
+impl<'p, 'a, F: Field, CS: Assignment<F> + SyncDeps> AssignmentPass<'p, 'a, F, CS> {
     fn new(plan: &'p mut V1Plan<'a, F, CS>) -> Self {
         AssignmentPass {
             plan,
@@ -384,7 +384,9 @@ impl<'r, 'a, F: Field, CS: Assignment<F> + 'a> V1Region<'r, 'a, F, CS> {
     }
 }
 
-impl<'r, 'a, F: Field, CS: Assignment<F> + 'a> RegionLayouter<F> for V1Region<'r, 'a, F, CS> {
+impl<'r, 'a, F: Field, CS: Assignment<F> + SyncDeps> SyncDeps for V1Region<'r, 'a, F, CS> {}
+
+impl<'r, 'a, F: Field, CS: Assignment<F> + SyncDeps> RegionLayouter<F> for V1Region<'r, 'a, F, CS> {
     fn enable_selector<'v>(
         &'v mut self,
         annotation: &'v (dyn Fn() -> String + 'v),
diff --git a/halo2_proofs/src/circuit/layouter.rs b/halo2_proofs/src/circuit/layouter.rs
index 4b9d74b6..e9eac92a 100644
--- a/halo2_proofs/src/circuit/layouter.rs
+++ b/halo2_proofs/src/circuit/layouter.rs
@@ -9,6 +9,14 @@ use ff::Field;
 use super::{Cell, RegionIndex, Value};
 use crate::plonk::{Advice, Any, Assigned, Column, Error, Fixed, Instance, Selector, TableColumn};
 
+/// Intermediate trait requirements for [`RegionLayouter`] when thread-safe regions are enabled.
+#[cfg(feature = "thread-safe-region")]
+pub trait SyncDeps: Send + Sync {}
+
+/// Intermediate trait requirements for [`RegionLayouter`].
+#[cfg(not(feature = "thread-safe-region"))]
+pub trait SyncDeps {}
+
 /// Helper trait for implementing a custom [`Layouter`].
 ///
 /// This trait is used for implementing region assignments:
@@ -39,7 +47,7 @@ use crate::plonk::{Advice, Any, Assigned, Column, Error, Fixed, Instance, Select
 /// `Chip::Config`).
 ///
 /// [`Layouter`]: super::Layouter
-pub trait RegionLayouter<F: Field>: fmt::Debug {
+pub trait RegionLayouter<F: Field>: fmt::Debug + SyncDeps {
     /// Enables a selector at the given offset.
     fn enable_selector<'v>(
         &'v mut self,
@@ -141,6 +149,8 @@ pub struct RegionShape {
     pub(super) row_count: usize,
 }
 
+impl SyncDeps for RegionShape {}
+
 /// The virtual column involved in a region. This includes concrete columns,
 /// as well as selectors that are not concrete columns at this stage.
 #[derive(Eq, PartialEq, Copy, Clone, Debug, Hash)]
diff --git a/halo2_proofs/src/dev.rs b/halo2_proofs/src/dev.rs
index 98b99012..a1bb0610 100644
--- a/halo2_proofs/src/dev.rs
+++ b/halo2_proofs/src/dev.rs
@@ -12,6 +12,7 @@ use ff::Field;
 use ff::FromUniformBytes;
 use group::Group;
 
+use crate::circuit::layouter::SyncDeps;
 use crate::plonk::permutation::keygen::Assembly;
 use crate::{
     circuit,
@@ -697,6 +698,9 @@ impl<F: FromUniformBytes<64> + Ord> MockProver<F> {
             v
         }));
 
+        #[cfg(feature = "thread-safe-region")]
+        prover.permutation.build_ordered_mapping();
+
         Ok(prover)
     }
 
@@ -996,7 +1000,7 @@ impl<F: FromUniformBytes<64> + Ord> MockProver<F> {
                         })
                         .collect::<Vec<_>>()
                 });
-
+        let mapping = self.permutation.mapping();
         // Check that permutations preserve the original values of the cells.
         let perm_errors = {
             // Original values of columns involved in the permutation.
@@ -1014,14 +1018,12 @@ impl<F: FromUniformBytes<64> + Ord> MockProver<F> {
             };
 
             // Iterate over each column of the permutation
-            self.permutation
-                .mapping()
-                .iter()
-                .enumerate()
-                .flat_map(move |(column, values)| {
-                    // Iterate over each row of the column to check that the cell's
-                    // value is preserved by the mapping.
-                    values.iter().enumerate().filter_map(move |(row, cell)| {
+            mapping.enumerate().flat_map(move |(column, values)| {
+                // Iterate over each row of the column to check that the cell's
+                // value is preserved by the mapping.
+                values
+                    .enumerate()
+                    .filter_map(move |(row, cell)| {
                         let original_cell = original(column, row);
                         let permuted_cell = original(cell.0, cell.1);
                         if original_cell == permuted_cell {
@@ -1039,7 +1041,8 @@ impl<F: FromUniformBytes<64> + Ord> MockProver<F> {
                             })
                         }
                     })
-                })
+                    .collect::<Vec<_>>()
+            })
         };
 
         let mut errors: Vec<_> = iter::empty()
@@ -1357,7 +1360,7 @@ impl<F: FromUniformBytes<64> + Ord> MockProver<F> {
                         })
                         .collect::<Vec<_>>()
                 });
-
+        let mapping = self.permutation.mapping();
         // Check that permutations preserve the original values of the cells.
         let perm_errors = {
             // Original values of columns involved in the permutation.
@@ -1375,36 +1378,31 @@ impl<F: FromUniformBytes<64> + Ord> MockProver<F> {
             };
 
             // Iterate over each column of the permutation
-            self.permutation
-                .mapping()
-                .iter()
-                .enumerate()
-                .flat_map(move |(column, values)| {
-                    // Iterate over each row of the column to check that the cell's
-                    // value is preserved by the mapping.
-                    values
-                        .par_iter()
-                        .enumerate()
-                        .filter_map(move |(row, cell)| {
-                            let original_cell = original(column, row);
-                            let permuted_cell = original(cell.0, cell.1);
-                            if original_cell == permuted_cell {
-                                None
-                            } else {
-                                let columns = self.cs.permutation.get_columns();
-                                let column = columns.get(column).unwrap();
-                                Some(VerifyFailure::Permutation {
-                                    column: (*column).into(),
-                                    location: FailureLocation::find(
-                                        &self.regions,
-                                        row,
-                                        Some(column).into_iter().cloned().collect(),
-                                    ),
-                                })
-                            }
-                        })
-                        .collect::<Vec<_>>()
-                })
+            mapping.enumerate().flat_map(move |(column, values)| {
+                // Iterate over each row of the column to check that the cell's
+                // value is preserved by the mapping.
+                values
+                    .enumerate()
+                    .filter_map(move |(row, cell)| {
+                        let original_cell = original(column, row);
+                        let permuted_cell = original(cell.0, cell.1);
+                        if original_cell == permuted_cell {
+                            None
+                        } else {
+                            let columns = self.cs.permutation.get_columns();
+                            let column = columns.get(column).unwrap();
+                            Some(VerifyFailure::Permutation {
+                                column: (*column).into(),
+                                location: FailureLocation::find(
+                                    &self.regions,
+                                    row,
+                                    Some(column).into_iter().cloned().collect(),
+                                ),
+                            })
+                        }
+                    })
+                    .collect::<Vec<_>>()
+            })
         };
 
         let mut errors: Vec<_> = iter::empty()
@@ -2096,3 +2094,5 @@ mod tests {
         )
     }
 }
+
+impl<F: Field> SyncDeps for MockProver<F> {}
diff --git a/halo2_proofs/src/dev/cost.rs b/halo2_proofs/src/dev/cost.rs
index 35e68540..f2c4e9f2 100644
--- a/halo2_proofs/src/dev/cost.rs
+++ b/halo2_proofs/src/dev/cost.rs
@@ -11,7 +11,7 @@ use ff::{Field, PrimeField};
 use group::prime::PrimeGroup;
 
 use crate::{
-    circuit::Value,
+    circuit::{layouter::SyncDeps, Value},
     plonk::{
         Advice, Any, Assigned, Assignment, Challenge, Circuit, Column, ConstraintSystem, Error,
         Fixed, FloorPlanner, Instance, Selector,
@@ -46,6 +46,8 @@ struct Assembly {
     selectors: Vec<Vec<bool>>,
 }
 
+impl SyncDeps for Assembly {}
+
 impl<F: Field> Assignment<F> for Assembly {
     fn enter_region<NR, N>(&mut self, _: N)
     where
diff --git a/halo2_proofs/src/dev/graph.rs b/halo2_proofs/src/dev/graph.rs
index b2739400..008d5840 100644
--- a/halo2_proofs/src/dev/graph.rs
+++ b/halo2_proofs/src/dev/graph.rs
@@ -80,6 +80,8 @@ struct Graph {
     current_namespace: Vec<usize>,
 }
 
+impl crate::dev::SyncDeps for Graph {}
+
 impl<F: Field> Assignment<F> for Graph {
     fn enter_region<NR, N>(&mut self, _: N)
     where
diff --git a/halo2_proofs/src/dev/graph/layout.rs b/halo2_proofs/src/dev/graph/layout.rs
index 59eced28..4782fd8a 100644
--- a/halo2_proofs/src/dev/graph/layout.rs
+++ b/halo2_proofs/src/dev/graph/layout.rs
@@ -355,6 +355,8 @@ struct Layout {
     selectors: Vec<Vec<bool>>,
 }
 
+impl crate::dev::SyncDeps for Layout {}
+
 impl Layout {
     fn new(k: u32, n: usize, num_selectors: usize) -> Self {
         Layout {
diff --git a/halo2_proofs/src/plonk/circuit.rs b/halo2_proofs/src/plonk/circuit.rs
index 97e28fc8..1c55eb5c 100644
--- a/halo2_proofs/src/plonk/circuit.rs
+++ b/halo2_proofs/src/plonk/circuit.rs
@@ -1,4 +1,5 @@
 use super::{lookup, permutation, Assigned, Error};
+use crate::circuit::layouter::SyncDeps;
 use crate::dev::metadata;
 use crate::{
     circuit::{Layouter, Region, Value},
@@ -736,7 +737,7 @@ pub trait FloorPlanner {
     /// - Perform any necessary setup or measurement tasks, which may involve one or more
     ///   calls to `Circuit::default().synthesize(config, &mut layouter)`.
     /// - Call `circuit.synthesize(config, &mut layouter)` exactly once.
-    fn synthesize<F: Field, CS: Assignment<F>, C: Circuit<F>>(
+    fn synthesize<F: Field, CS: Assignment<F> + SyncDeps, C: Circuit<F>>(
         cs: &mut CS,
         circuit: &C,
         config: C::Config,
diff --git a/halo2_proofs/src/plonk/keygen.rs b/halo2_proofs/src/plonk/keygen.rs
index 16ef4c6f..91922b69 100644
--- a/halo2_proofs/src/plonk/keygen.rs
+++ b/halo2_proofs/src/plonk/keygen.rs
@@ -16,7 +16,7 @@ use super::{
 };
 use crate::{
     arithmetic::{parallelize, CurveAffine},
-    circuit::Value,
+    circuit::{layouter::SyncDeps, Value},
     poly::{
         batch_invert_assigned,
         commitment::{Blind, Params, MSM},
@@ -61,6 +61,8 @@ struct Assembly<F: Field> {
     _marker: std::marker::PhantomData<F>,
 }
 
+impl<F: Field> SyncDeps for Assembly<F> {}
+
 impl<F: Field> Assignment<F> for Assembly<F> {
     fn enter_region<NR, N>(&mut self, _: N)
     where
diff --git a/halo2_proofs/src/plonk/permutation/keygen.rs b/halo2_proofs/src/plonk/permutation/keygen.rs
index bcc0d02b..57882fad 100644
--- a/halo2_proofs/src/plonk/permutation/keygen.rs
+++ b/halo2_proofs/src/plonk/permutation/keygen.rs
@@ -1,5 +1,11 @@
+use std::collections::{BTreeMap, BTreeSet, HashMap, HashSet};
+
 use ff::{Field, PrimeField};
 use group::Curve;
+use rayon::prelude::{
+    IndexedParallelIterator, IntoParallelIterator, IntoParallelRefIterator,
+    IntoParallelRefMutIterator, ParallelIterator, ParallelSliceMut,
+};
 
 use super::{Argument, ProvingKey, VerifyingKey};
 use crate::{
@@ -11,6 +17,7 @@ use crate::{
     },
 };
 
+#[cfg(not(feature = "thread-safe-region"))]
 /// Struct that accumulates all the necessary data in order to construct the permutation argument.
 #[derive(Clone, Debug, PartialEq, Eq)]
 pub struct Assembly {
@@ -24,6 +31,7 @@ pub struct Assembly {
     sizes: Vec<Vec<usize>>,
 }
 
+#[cfg(not(feature = "thread-safe-region"))]
 impl Assembly {
     pub(crate) fn new(n: usize, p: &Argument) -> Self {
         // Initialize the copy vector to keep track of copy constraints in all
@@ -108,145 +116,340 @@ impl Assembly {
         domain: &EvaluationDomain<C::Scalar>,
         p: &Argument,
     ) -> VerifyingKey<C> {
-        // Compute [omega^0, omega^1, ..., omega^{params.n - 1}]
-        let mut omega_powers = vec![C::Scalar::ZERO; params.n() as usize];
-        {
-            let omega = domain.get_omega();
-            parallelize(&mut omega_powers, |o, start| {
-                let mut cur = omega.pow_vartime(&[start as u64]);
-                for v in o.iter_mut() {
-                    *v = cur;
-                    cur *= &omega;
-                }
-            })
+        build_vk(params, domain, p, |i, j| self.mapping[i][j])
+    }
+
+    pub(crate) fn build_pk<'params, C: CurveAffine, P: Params<'params, C>>(
+        self,
+        params: &P,
+        domain: &EvaluationDomain<C::Scalar>,
+        p: &Argument,
+    ) -> ProvingKey<C> {
+        build_pk(params, domain, p, |i, j| self.mapping[i][j])
+    }
+
+    /// Returns columns that participate in the permutation argument.
+    pub fn columns(&self) -> &[Column<Any>] {
+        &self.columns
+    }
+
+    /// Returns mappings of the copies.
+    pub fn mapping(
+        &self,
+    ) -> impl Iterator<Item = impl IndexedParallelIterator<Item = (usize, usize)> + '_> {
+        self.mapping.iter().map(|c| c.par_iter().copied())
+    }
+}
+
+#[cfg(feature = "thread-safe-region")]
+/// Struct that accumulates all the necessary data in order to construct the permutation argument.
+#[derive(Clone, Debug, PartialEq, Eq)]
+pub struct Assembly {
+    /// Columns that participate on the copy permutation argument.
+    columns: Vec<Column<Any>>,
+    /// Mapping of the actual copies done.
+    cycles: Vec<Vec<(usize, usize)>>,
+    /// Mapping of the actual copies done.
+    ordered_cycles: Vec<BTreeSet<(usize, usize)>>,
+    /// Mapping of the actual copies done.
+    aux: HashMap<(usize, usize), usize>,
+    /// total length of a column
+    col_len: usize,
+    /// number of columns
+    num_cols: usize,
+}
+
+#[cfg(feature = "thread-safe-region")]
+impl Assembly {
+    pub(crate) fn new(n: usize, p: &Argument) -> Self {
+        Assembly {
+            columns: p.columns.clone(),
+            cycles: Vec::with_capacity(n),
+            ordered_cycles: Vec::with_capacity(n),
+            aux: HashMap::new(),
+            col_len: n,
+            num_cols: p.columns.len(),
         }
+    }
 
-        // Compute [omega_powers * \delta^0, omega_powers * \delta^1, ..., omega_powers * \delta^m]
-        let mut deltaomega = vec![omega_powers; p.columns.len()];
-        {
-            parallelize(&mut deltaomega, |o, start| {
-                let mut cur = C::Scalar::DELTA.pow_vartime(&[start as u64]);
-                for omega_powers in o.iter_mut() {
-                    for v in omega_powers {
-                        *v *= &cur;
-                    }
-                    cur *= &<C::Scalar as PrimeField>::DELTA;
-                }
-            });
+    pub(crate) fn copy(
+        &mut self,
+        left_column: Column<Any>,
+        left_row: usize,
+        right_column: Column<Any>,
+        right_row: usize,
+    ) -> Result<(), Error> {
+        let left_column = self
+            .columns
+            .iter()
+            .position(|c| c == &left_column)
+            .ok_or(Error::ColumnNotInPermutation(left_column))?;
+        let right_column = self
+            .columns
+            .iter()
+            .position(|c| c == &right_column)
+            .ok_or(Error::ColumnNotInPermutation(right_column))?;
+
+        // Check bounds
+        if left_row >= self.col_len || right_row >= self.col_len {
+            return Err(Error::BoundsFailure);
         }
 
-        // Computes the permutation polynomial based on the permutation
-        // description in the assembly.
-        let mut permutations = vec![domain.empty_lagrange(); p.columns.len()];
-        {
-            parallelize(&mut permutations, |o, start| {
-                for (x, permutation_poly) in o.iter_mut().enumerate() {
-                    let i = start + x;
-                    for (j, p) in permutation_poly.iter_mut().enumerate() {
-                        let (permuted_i, permuted_j) = self.mapping[i][j];
-                        *p = deltaomega[permuted_i][permuted_j];
-                    }
-                }
-            });
+        let left_cycle = self.aux.get(&(left_column, left_row));
+        let right_cycle = self.aux.get(&(right_column, right_row));
+
+        // extract cycle elements
+        let right_cycle_elems = match right_cycle {
+            Some(i) => {
+                let entry = self.cycles[*i].clone();
+                self.cycles[*i] = vec![];
+                entry
+            }
+            None => [(right_column, right_row)].into(),
+        };
+
+        assert!(right_cycle_elems.contains(&(right_column, right_row)));
+
+        // merge cycles
+        let cycle_idx = match left_cycle {
+            Some(i) => {
+                let entry = &mut self.cycles[*i];
+                entry.extend(right_cycle_elems.clone());
+                *i
+            }
+            // if they were singletons -- create a new cycle entry
+            None => {
+                let mut set: Vec<(usize, usize)> = right_cycle_elems.clone();
+                set.push((left_column, left_row));
+                self.cycles.push(set);
+                let cycle_idx = self.cycles.len() - 1;
+                self.aux.insert((left_column, left_row), cycle_idx);
+                cycle_idx
+            }
+        };
+
+        let index_updates = vec![cycle_idx; right_cycle_elems.len()].into_iter();
+        let updates = right_cycle_elems.into_iter().zip(index_updates);
+
+        self.aux.extend(updates);
+
+        Ok(())
+    }
+
+    /// Builds the ordered mapping of the cycles.
+    /// This will only get executed once.
+    pub fn build_ordered_mapping(&mut self) {
+        // will only get called once
+        if self.ordered_cycles.is_empty() && !self.cycles.is_empty() {
+            self.ordered_cycles = self
+                .cycles
+                .par_iter_mut()
+                .map(|col| {
+                    let mut set = BTreeSet::new();
+                    set.extend(col.clone());
+                    // free up memory
+                    *col = vec![];
+                    set
+                })
+                .collect();
         }
+    }
 
-        // Pre-compute commitments for the URS.
-        let mut commitments = Vec::with_capacity(p.columns.len());
-        for permutation in &permutations {
-            // Compute commitment to permutation polynomial
-            commitments.push(
-                params
-                    .commit_lagrange(permutation, Blind::default())
-                    .to_affine(),
-            );
+    fn mapping_at_idx(&self, col: usize, row: usize) -> (usize, usize) {
+        assert!(
+            !self.ordered_cycles.is_empty() || self.cycles.is_empty(),
+            "cycles have not been ordered"
+        );
+
+        if let Some(cycle_idx) = self.aux.get(&(col, row)) {
+            let cycle = &self.ordered_cycles[*cycle_idx];
+            let mut cycle_iter = cycle.range((
+                std::ops::Bound::Excluded((col, row)),
+                std::ops::Bound::Unbounded,
+            ));
+            // point to the next node in the cycle
+            match cycle_iter.next() {
+                Some((i, j)) => (*i, *j),
+                // wrap back around to the first element which SHOULD exist
+                None => *(cycle.iter().next().unwrap()),
+            }
+        // is a singleton
+        } else {
+            (col, row)
         }
+    }
 
-        VerifyingKey { commitments }
+    pub(crate) fn build_vk<'params, C: CurveAffine, P: Params<'params, C>>(
+        &mut self,
+        params: &P,
+        domain: &EvaluationDomain<C::Scalar>,
+        p: &Argument,
+    ) -> VerifyingKey<C> {
+        self.build_ordered_mapping();
+        build_vk(params, domain, p, |i, j| self.mapping_at_idx(i, j))
     }
 
     pub(crate) fn build_pk<'params, C: CurveAffine, P: Params<'params, C>>(
-        self,
+        &mut self,
         params: &P,
         domain: &EvaluationDomain<C::Scalar>,
         p: &Argument,
     ) -> ProvingKey<C> {
-        // Compute [omega^0, omega^1, ..., omega^{params.n - 1}]
-        let mut omega_powers = vec![C::Scalar::ZERO; params.n() as usize];
-        {
-            let omega = domain.get_omega();
-            parallelize(&mut omega_powers, |o, start| {
-                let mut cur = omega.pow_vartime(&[start as u64]);
-                for v in o.iter_mut() {
-                    *v = cur;
-                    cur *= &omega;
-                }
-            })
-        }
+        self.build_ordered_mapping();
+        build_pk(params, domain, p, |i, j| self.mapping_at_idx(i, j))
+    }
 
-        // Compute [omega_powers * \delta^0, omega_powers * \delta^1, ..., omega_powers * \delta^m]
-        let mut deltaomega = vec![omega_powers; p.columns.len()];
-        {
-            parallelize(&mut deltaomega, |o, start| {
-                let mut cur = C::Scalar::DELTA.pow_vartime(&[start as u64]);
-                for omega_powers in o.iter_mut() {
-                    for v in omega_powers {
-                        *v *= &cur;
-                    }
-                    cur *= &C::Scalar::DELTA;
-                }
-            });
-        }
+    /// Returns columns that participate in the permutation argument.
+    pub fn columns(&self) -> &[Column<Any>] {
+        &self.columns
+    }
 
-        // Compute permutation polynomials, convert to coset form.
-        let mut permutations = vec![domain.empty_lagrange(); p.columns.len()];
-        {
-            parallelize(&mut permutations, |o, start| {
-                for (x, permutation_poly) in o.iter_mut().enumerate() {
-                    let i = start + x;
-                    for (j, p) in permutation_poly.iter_mut().enumerate() {
-                        let (permuted_i, permuted_j) = self.mapping[i][j];
-                        *p = deltaomega[permuted_i][permuted_j];
-                    }
-                }
-            });
-        }
+    /// Returns mappings of the copies.
+    pub fn mapping(
+        &self,
+    ) -> impl Iterator<Item = impl IndexedParallelIterator<Item = (usize, usize)> + '_> {
+        (0..self.num_cols).map(move |i| {
+            (0..self.col_len)
+                .into_par_iter()
+                .map(move |j| self.mapping_at_idx(i, j))
+        })
+    }
+}
 
-        let mut polys = vec![domain.empty_coeff(); p.columns.len()];
-        {
-            parallelize(&mut polys, |o, start| {
-                for (x, poly) in o.iter_mut().enumerate() {
-                    let i = start + x;
-                    let permutation_poly = permutations[i].clone();
-                    *poly = domain.lagrange_to_coeff(permutation_poly);
+pub(crate) fn build_pk<'params, C: CurveAffine, P: Params<'params, C>>(
+    params: &P,
+    domain: &EvaluationDomain<C::Scalar>,
+    p: &Argument,
+    mapping: impl Fn(usize, usize) -> (usize, usize) + Sync,
+) -> ProvingKey<C> {
+    // Compute [omega^0, omega^1, ..., omega^{params.n - 1}]
+    let mut omega_powers = vec![C::Scalar::ZERO; params.n() as usize];
+    {
+        let omega = domain.get_omega();
+        parallelize(&mut omega_powers, |o, start| {
+            let mut cur = omega.pow_vartime(&[start as u64]);
+            for v in o.iter_mut() {
+                *v = cur;
+                cur *= &omega;
+            }
+        })
+    }
+
+    // Compute [omega_powers * \delta^0, omega_powers * \delta^1, ..., omega_powers * \delta^m]
+    let mut deltaomega = vec![omega_powers; p.columns.len()];
+    {
+        parallelize(&mut deltaomega, |o, start| {
+            let mut cur = C::Scalar::DELTA.pow_vartime(&[start as u64]);
+            for omega_powers in o.iter_mut() {
+                for v in omega_powers {
+                    *v *= &cur;
                 }
-            });
-        }
+                cur *= &C::Scalar::DELTA;
+            }
+        });
+    }
 
-        let mut cosets = vec![domain.empty_extended(); p.columns.len()];
-        {
-            parallelize(&mut cosets, |o, start| {
-                for (x, coset) in o.iter_mut().enumerate() {
-                    let i = start + x;
-                    let poly = polys[i].clone();
-                    *coset = domain.coeff_to_extended(poly);
+    // Compute permutation polynomials, convert to coset form.
+    let mut permutations = vec![domain.empty_lagrange(); p.columns.len()];
+    {
+        parallelize(&mut permutations, |o, start| {
+            for (x, permutation_poly) in o.iter_mut().enumerate() {
+                let i = start + x;
+                for (j, p) in permutation_poly.iter_mut().enumerate() {
+                    let (permuted_i, permuted_j) = mapping(i, j);
+                    *p = deltaomega[permuted_i][permuted_j];
                 }
-            });
-        }
+            }
+        });
+    }
 
-        ProvingKey {
-            permutations,
-            polys,
-            cosets,
-        }
+    let mut polys = vec![domain.empty_coeff(); p.columns.len()];
+    {
+        parallelize(&mut polys, |o, start| {
+            for (x, poly) in o.iter_mut().enumerate() {
+                let i = start + x;
+                let permutation_poly = permutations[i].clone();
+                *poly = domain.lagrange_to_coeff(permutation_poly);
+            }
+        });
     }
 
-    /// Returns columns that participate on the permutation argument.
-    pub fn columns(&self) -> &[Column<Any>] {
-        &self.columns
+    let mut cosets = vec![domain.empty_extended(); p.columns.len()];
+    {
+        parallelize(&mut cosets, |o, start| {
+            for (x, coset) in o.iter_mut().enumerate() {
+                let i = start + x;
+                let poly = polys[i].clone();
+                *coset = domain.coeff_to_extended(poly);
+            }
+        });
     }
 
-    /// Returns mappings of the copies.
-    pub fn mapping(&self) -> &[Vec<(usize, usize)>] {
-        &self.mapping
+    ProvingKey {
+        permutations,
+        polys,
+        cosets,
+    }
+}
+
+pub(crate) fn build_vk<'params, C: CurveAffine, P: Params<'params, C>>(
+    params: &P,
+    domain: &EvaluationDomain<C::Scalar>,
+    p: &Argument,
+    mapping: impl Fn(usize, usize) -> (usize, usize) + Sync,
+) -> VerifyingKey<C> {
+    // Compute [omega^0, omega^1, ..., omega^{params.n - 1}]
+    let mut omega_powers = vec![C::Scalar::ZERO; params.n() as usize];
+    {
+        let omega = domain.get_omega();
+        parallelize(&mut omega_powers, |o, start| {
+            let mut cur = omega.pow_vartime(&[start as u64]);
+            for v in o.iter_mut() {
+                *v = cur;
+                cur *= &omega;
+            }
+        })
+    }
+
+    // Compute [omega_powers * \delta^0, omega_powers * \delta^1, ..., omega_powers * \delta^m]
+    let mut deltaomega = vec![omega_powers; p.columns.len()];
+    {
+        parallelize(&mut deltaomega, |o, start| {
+            let mut cur = C::Scalar::DELTA.pow_vartime(&[start as u64]);
+            for omega_powers in o.iter_mut() {
+                for v in omega_powers {
+                    *v *= &cur;
+                }
+                cur *= &<C::Scalar as PrimeField>::DELTA;
+            }
+        });
+    }
+
+    // Computes the permutation polynomial based on the permutation
+    // description in the assembly.
+    let mut permutations = vec![domain.empty_lagrange(); p.columns.len()];
+    {
+        parallelize(&mut permutations, |o, start| {
+            for (x, permutation_poly) in o.iter_mut().enumerate() {
+                let i = start + x;
+                for (j, p) in permutation_poly.iter_mut().enumerate() {
+                    let (permuted_i, permuted_j) = mapping(i, j);
+                    *p = deltaomega[permuted_i][permuted_j];
+                }
+            }
+        });
+    }
+
+    // Pre-compute commitments for the URS.
+    let mut commitments = Vec::with_capacity(p.columns.len());
+    for permutation in &permutations {
+        // Compute commitment to permutation polynomial
+        commitments.push(
+            params
+                .commit_lagrange(permutation, Blind::default())
+                .to_affine(),
+        );
     }
+
+    VerifyingKey { commitments }
 }
diff --git a/halo2_proofs/src/plonk/prover.rs b/halo2_proofs/src/plonk/prover.rs
index 1258f09f..6314481d 100644
--- a/halo2_proofs/src/plonk/prover.rs
+++ b/halo2_proofs/src/plonk/prover.rs
@@ -18,6 +18,7 @@ use super::{
     lookup, permutation, vanishing, ChallengeBeta, ChallengeGamma, ChallengeTheta, ChallengeX,
     ChallengeY, Error, Expression, ProvingKey,
 };
+use crate::circuit::layouter::SyncDeps;
 use crate::{
     arithmetic::{eval_polynomial, CurveAffine},
     circuit::Value,
@@ -153,6 +154,8 @@ where
         _marker: std::marker::PhantomData<F>,
     }
 
+    impl<'a, F: Field> SyncDeps for WitnessCollection<'a, F> {}
+
     impl<'a, F: Field> Assignment<F> for WitnessCollection<'a, F> {
         fn enter_region<NR, N>(&mut self, _: N)
         where