diff --git a/apps/qsim_amplitudes.cc b/apps/qsim_amplitudes.cc
index 2816016f..70851c13 100644
--- a/apps/qsim_amplitudes.cc
+++ b/apps/qsim_amplitudes.cc
@@ -23,10 +23,10 @@
 
 #include "../lib/bitstring.h"
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io_file.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsim.h"
 #include "../lib/simmux.h"
 #include "../lib/util.h"
@@ -153,7 +153,7 @@ int main(int argc, char* argv[]) {
     return 1;
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
 
diff --git a/apps/qsim_base.cc b/apps/qsim_base.cc
index 84c63ccb..17ef2e6a 100644
--- a/apps/qsim_base.cc
+++ b/apps/qsim_base.cc
@@ -20,10 +20,10 @@
 #include <string>
 
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io_file.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsim.h"
 #include "../lib/simmux.h"
 
@@ -105,7 +105,7 @@ int main(int argc, char* argv[]) {
     return 1;
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
   using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
diff --git a/apps/qsim_von_neumann.cc b/apps/qsim_von_neumann.cc
index e22de4cf..96d0c917 100644
--- a/apps/qsim_von_neumann.cc
+++ b/apps/qsim_von_neumann.cc
@@ -22,10 +22,10 @@
 #include <string>
 
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io_file.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsim.h"
 #include "../lib/simmux.h"
 
@@ -90,7 +90,7 @@ int main(int argc, char* argv[]) {
     return 1;
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
 
@@ -105,8 +105,8 @@ int main(int argc, char* argv[]) {
     };
 
     auto size = state_space.Size(state);
-    double entropy = -ParallelFor::RunReduce(opt.num_threads, size, f, Op(),
-                                             state_space, state);
+    double entropy = -For::RunReduce(opt.num_threads, size, f, Op(),
+                                     state_space, state);
     IO::messagef("entropy=%g\n", entropy);
   };
 
diff --git a/apps/qsimh_amplitudes.cc b/apps/qsimh_amplitudes.cc
index 2d2c9ee6..e28a18c9 100644
--- a/apps/qsimh_amplitudes.cc
+++ b/apps/qsimh_amplitudes.cc
@@ -23,10 +23,10 @@
 
 #include "../lib/bitstring.h"
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io_file.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsimh.h"
 #include "../lib/simmux.h"
 #include "../lib/util.h"
@@ -186,9 +186,9 @@ int main(int argc, char* argv[]) {
     return 1;
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using HybridSimulator = HybridSimulator<IO, GateQSim<float>, BasicGateFuser,
-                                          Simulator, ParallelFor>;
+                                          Simulator, For>;
   using Runner = QSimHRunner<IO, HybridSimulator>;
 
   Runner::Parameter param;
diff --git a/apps/qsimh_base.cc b/apps/qsimh_base.cc
index 99ce15a1..c285d039 100644
--- a/apps/qsimh_base.cc
+++ b/apps/qsimh_base.cc
@@ -23,10 +23,10 @@
 
 #include "../lib/bitstring.h"
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io_file.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsimh.h"
 #include "../lib/simmux.h"
 #include "../lib/util.h"
@@ -151,9 +151,9 @@ int main(int argc, char* argv[]) {
     bitstrings.push_back(i);
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using HybridSimulator = HybridSimulator<IO, GateQSim<float>, BasicGateFuser,
-                                          Simulator, ParallelFor>;
+                                          Simulator, For>;
   using Runner = QSimHRunner<IO, HybridSimulator>;
 
   Runner::Parameter param;
diff --git a/docs/type_reference.md b/docs/type_reference.md
index e99bbfc7..c8485f9e 100644
--- a/docs/type_reference.md
+++ b/docs/type_reference.md
@@ -30,6 +30,7 @@ they fulfill the same expectations.
 | Op (util.h)             | [`to_int` (in `Options`)](apps/qsim_amplitudes.cc)    |
 | ParallelFor             | [`ParallelFor`](lib/parfor.h)                     |
 | Params                  | Vector of `fp_type`.                              |
+| SequentialFor           | [`SequentialFor`](lib/seqfor.h)                   |
 | Simulator               | [`SimulatorAVX`](lib/simulator_avx.h)             |
 | State                   | Unique pointer to `fp_type`.                      |
 | StateSpace              | [`StateSpace`](lib/statespace.h)                  |
diff --git a/docs/usage.md b/docs/usage.md
index c0d8fc14..54a1030f 100644
--- a/docs/usage.md
+++ b/docs/usage.md
@@ -135,10 +135,3 @@ Example:
 ./qsimh_amplitudes.x -c ../circuits/circuit_q40 -d 47 -k 0,1,2,3,4,5,6,7,8,9,10,13,14,15,16,17,23,24 -t 8 -w 1 -p 9 -r 4 -i ../circuits/bitstrings_q40_s1 -o ampl_q40_s1_w1 -v 1
 ...
 ```
-
-# Running without AVX
-
-SimulatorBasic can be used on older systems that do not support AVX
-instructions. In this case, simulator_basic.h should be included instead of
-simulator_avx.h and SimulatorAVX should be replaced by
-SimulatorBasic<ParallelFor, float>.
diff --git a/lib/BUILD b/lib/BUILD
index daf03b3b..bb2297b9 100644
--- a/lib/BUILD
+++ b/lib/BUILD
@@ -9,6 +9,7 @@ cc_library(
         "bitstring.h",
         "circuit_qsim_parser.h",
         "circuit.h",
+        "formux.h",
         "fuser_basic.h",
         "fuser.h",
         "gate.h",
@@ -41,6 +42,7 @@ cc_library(
     hdrs = [
         "circuit.h",
         "circuit_qsim_parser.h",
+        "formux.h",
         "fuser.h",
         "fuser_basic.h",
         "gate.h",
@@ -51,6 +53,7 @@ cc_library(
         "matrix.h",
         "parfor.h",
         "run_qsim.h",
+        "seqfor.h",
         "simmux.h",
         "simulator_avx.h",
         "simulator_basic.h",
@@ -69,6 +72,7 @@ cc_library(
     hdrs = [
         "circuit.h",
         "circuit_qsim_parser.h",
+        "formux.h",
         "fuser.h",
         "fuser_basic.h",
         "gate.h",
@@ -80,6 +84,7 @@ cc_library(
         "matrix.h",
         "parfor.h",
         "run_qsimh.h",
+        "seqfor.h",
         "simmux.h",
         "simulator_avx.h",
         "simulator_basic.h",
@@ -139,6 +144,16 @@ cc_library(
     hdrs = ["seqfor.h"],
 )
 
+# All three state-vector simulators with multiplexer
+cc_library(
+    name = "formux",
+    hdrs = ["formux.h"],
+    deps = [
+        ":parfor",
+        ":seqfor",
+    ],
+)
+
 ### Gate libraries ###
 
 cc_library(
diff --git a/lib/formux.h b/lib/formux.h
new file mode 100644
index 00000000..4401e9bb
--- /dev/null
+++ b/lib/formux.h
@@ -0,0 +1,30 @@
+// Copyright 2019 Google LLC. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef FORMUX_H_
+#define FORMUX_H_
+
+#ifdef _OPENMP
+# include "parfor.h"
+  namespace qsim {
+    using For = ParallelFor;
+  }
+#else
+# include "seqfor.h"
+  namespace qsim {
+    using For = SequentialFor;
+  }
+#endif
+
+#endif  // FORMUX_H_
diff --git a/lib/hybrid.h b/lib/hybrid.h
index 5bd99447..4910582c 100644
--- a/lib/hybrid.h
+++ b/lib/hybrid.h
@@ -27,7 +27,7 @@ namespace qsim {
 
 // Hybrid Feynmann-Schrodiner simulator.
 template <typename IO, typename GateT, template <typename> class FuserT,
-          typename Simulator, typename ParallelFor>
+          typename Simulator, typename For>
 struct HybridSimulator final {
  public:
   using Gate = GateT;
@@ -368,9 +368,8 @@ struct HybridSimulator final {
         };
 
         // Collect results.
-        ParallelFor::Run(param.num_threads, results.size(), f,
-                         sspace0, sspace1, *rstate0, *rstate1, indices,
-                         results);
+        For::Run(param.num_threads, results.size(), f, sspace0, sspace1,
+                 *rstate0, *rstate1, indices, results);
       }
     }
 
diff --git a/lib/simmux.h b/lib/simmux.h
index 067ae83e..626eab1c 100644
--- a/lib/simmux.h
+++ b/lib/simmux.h
@@ -18,20 +18,20 @@
 #ifdef __AVX2__
 # include "simulator_avx.h"
   namespace qsim {
-    template <typename ParallelFor>
-    using Simulator = SimulatorAVX<ParallelFor>;
+    template <typename For>
+    using Simulator = SimulatorAVX<For>;
   }
 #elif __SSE4_1__
 # include "simulator_sse.h"
   namespace qsim {
-    template <typename ParallelFor>
-    using Simulator = SimulatorSSE<ParallelFor>;
+    template <typename For>
+    using Simulator = SimulatorSSE<For>;
   }
 #else
 # include "simulator_basic.h"
   namespace qsim {
-    template <typename ParallelFor>
-    using Simulator = SimulatorBasic<ParallelFor>;
+    template <typename For>
+    using Simulator = SimulatorBasic<For>;
   }
 #endif
 
diff --git a/lib/simulator_avx.h b/lib/simulator_avx.h
index dc68510b..f3e2e943 100644
--- a/lib/simulator_avx.h
+++ b/lib/simulator_avx.h
@@ -25,10 +25,10 @@
 namespace qsim {
 
 // Quantum circuit simulator with AVX vectorization.
-template <typename ParallelFor>
+template <typename For>
 class SimulatorAVX final {
  public:
-  using StateSpace = StateSpaceAVX<ParallelFor>;
+  using StateSpace = StateSpaceAVX<For>;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
 
@@ -130,8 +130,7 @@ class SimulatorAVX final {
       _mm256_store_ps(rstate + p + 8, in);
     };
 
-    ParallelFor::Run(num_threads_, sizei / 16, f,
-                     sizek, mask0, mask1, matrix, rstate);
+    For::Run(num_threads_, sizei / 16, f, sizek, mask0, mask1, matrix, rstate);
   }
 
   // Applies a single-qubit gate for qubit <= 2.
@@ -224,8 +223,8 @@ class SimulatorAVX final {
       _mm256_store_ps(rstate + p + 8, in);
     };
 
-    ParallelFor::Run(num_threads_, std::max(uint64_t{1}, sizei / 16), f,
-                     q0, ml, matrix, rstate);
+    For::Run(num_threads_, std::max(uint64_t{1}, sizei / 16), f, q0, ml,
+             matrix, rstate);
   }
 
   // Applies two-qubit gate for qubit0 > 2 and qubit1 > 2.
@@ -377,8 +376,8 @@ class SimulatorAVX final {
       _mm256_store_ps(rstate + p + 8, in);
     };
 
-    ParallelFor::Run(num_threads_, sizei / 16, f,
-                     sizej, sizek, mask0, mask1, mask2, matrix, rstate);
+    For::Run(num_threads_, sizei / 16, f, sizej, sizek, mask0, mask1, mask2,
+             matrix, rstate);
   }
 
   // Applies a two-qubit gate for qubit0 <= 2 and qubit1 > 2.
@@ -587,8 +586,8 @@ class SimulatorAVX final {
       _mm256_store_ps(rstate + p + 8, in);
     };
 
-    ParallelFor::Run(num_threads_, sizei / 16, f,
-                     sizej, mask0, mask1, q0, ml, matrix, rstate);
+    For::Run(num_threads_, sizei / 16, f, sizej, mask0, mask1, q0, ml,
+             matrix, rstate);
   }
 
   // Applies a two-qubit gate for qubit0 <= 2 and qubit1 <= 2.
@@ -815,8 +814,8 @@ class SimulatorAVX final {
       _mm256_store_ps(rstate + p + 8, in);
     };
 
-    ParallelFor::Run(num_threads_, std::max(uint64_t{1}, sizei / 16), f,
-                     q, ml1, ml2, ml3, matrix, rstate);
+    For::Run(num_threads_, std::max(uint64_t{1}, sizei / 16), f, q,
+             ml1, ml2, ml3, matrix, rstate);
   }
 
   unsigned num_qubits_;
diff --git a/lib/simulator_basic.h b/lib/simulator_basic.h
index 5b98c985..f678ee76 100644
--- a/lib/simulator_basic.h
+++ b/lib/simulator_basic.h
@@ -22,10 +22,10 @@
 namespace qsim {
 
 // Quantim circuit simulator without vectorization.
-template <typename ParallelFor, typename FP = float>
+template <typename For, typename FP = float>
 class SimulatorBasic final {
  public:
-  using StateSpace = StateSpaceBasic<ParallelFor, FP>;
+  using StateSpace = StateSpaceBasic<For, FP>;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
 
@@ -66,8 +66,7 @@ class SimulatorBasic final {
       rstate[si1 + 1] = s0r * u[5] + s0i * u[4] + s1r * u[7] + s1i * u[6];
     };
 
-    ParallelFor::Run(num_threads_, sizei / 2, f,
-                     sizek, mask0, mask1, matrix, rstate);
+    For::Run(num_threads_, sizei / 2, f, sizek, mask0, mask1, matrix, rstate);
   }
 
   /**
@@ -130,8 +129,8 @@ class SimulatorBasic final {
           + s2r * u[29] + s2i * u[28] + s3r * u[31] + s3i * u[30];
     };
 
-    ParallelFor::Run(num_threads_, sizei / 2, f,
-                     sizej, sizek, mask0, mask1, mask2, matrix, rstate);
+    For::Run(num_threads_, sizei / 2, f, sizej, sizek, mask0, mask1, mask2,
+             matrix, rstate);
   }
 
  private:
diff --git a/lib/simulator_sse.h b/lib/simulator_sse.h
index 446499cf..f70a51c2 100644
--- a/lib/simulator_sse.h
+++ b/lib/simulator_sse.h
@@ -26,10 +26,10 @@
 namespace qsim {
 
 // Quantum circuit simulator with SSE vectorization.
-template <typename ParallelFor>
+template <typename For>
 class SimulatorSSE final {
  public:
-  using StateSpace = StateSpaceSSE<ParallelFor>;
+  using StateSpace = StateSpaceSSE<For>;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
 
@@ -131,8 +131,7 @@ class SimulatorSSE final {
       _mm_store_ps(rstate + p + 4, in);
     };
 
-    ParallelFor::Run(num_threads_, sizei / 8, f,
-                     sizek, mask0, mask1, matrix, rstate);
+    For::Run(num_threads_, sizei / 8, f, sizek, mask0, mask1, matrix, rstate);
   }
 
   // Applies a single-qubit gate for qubit <= 1.
@@ -213,8 +212,8 @@ class SimulatorSSE final {
       _mm_store_ps(rstate + p + 4, in);
     };
 
-    ParallelFor::Run(num_threads_, std::max(uint64_t{1}, sizei / 8), f,
-                     q0, matrix, rstate);
+    For::Run(num_threads_, std::max(uint64_t{1}, sizei / 8), f, q0,
+             matrix, rstate);
   }
 
   // Applies two-qubit gate for qubit0 > 1 and qubit1 > 1.
@@ -366,8 +365,8 @@ class SimulatorSSE final {
       _mm_store_ps(rstate + p + 4, in);
     };
 
-    ParallelFor::Run(num_threads_, sizei / 8, f,
-                     sizej, sizek, mask0, mask1, mask2, matrix, rstate);
+    For::Run(num_threads_, sizei / 8, f, sizej, sizek, mask0, mask1, mask2,
+             matrix, rstate);
   }
 
   // Applies a two-qubit gate for qubit0 <= 1 and qubit1 > 1.
@@ -558,8 +557,8 @@ class SimulatorSSE final {
       _mm_store_ps(rstate + p + 4, in);
     };
 
-    ParallelFor::Run(num_threads_, sizei / 8, f,
-                     sizej, mask0, mask1, q0, matrix, rstate);
+    For::Run(num_threads_, sizei / 8, f, sizej, mask0, mask1, q0,
+             matrix, rstate);
   }
 
   // Applies a two-qubit gate for qubit0 = 0 and qubit1 = 1.
@@ -602,7 +601,7 @@ class SimulatorSSE final {
           + s2r * u[29] + s2i * u[28] + s3r * u[31] + s3i * u[30];
     };
 
-    ParallelFor::Run(num_threads_, sizei / 8, f, matrix, rstate);
+    For::Run(num_threads_, sizei / 8, f, matrix, rstate);
   }
 
   unsigned num_qubits_;
diff --git a/lib/statespace.h b/lib/statespace.h
index fcf5974f..8cd68b56 100644
--- a/lib/statespace.h
+++ b/lib/statespace.h
@@ -32,7 +32,7 @@ inline void do_not_free(void*) noexcept {}
 }  // namespace detail
 
 // Routines for state-vector manipulations.
-template <typename ParallelFor, typename FP>
+template <typename For, typename FP>
 class StateSpace {
  public:
   using fp_type = FP;
@@ -90,7 +90,7 @@ class StateSpace {
       dest.get()[i] = src.get()[i];
     };
 
-    ParallelFor::Run(num_threads_, raw_size_, f, src, dest);
+    For::Run(num_threads_, raw_size_, f, src, dest);
   }
 
  protected:
diff --git a/lib/statespace_avx.h b/lib/statespace_avx.h
index 0893072f..36d94508 100644
--- a/lib/statespace_avx.h
+++ b/lib/statespace_avx.h
@@ -32,9 +32,9 @@ namespace qsim {
 // State is a vectorized sequence of eight real components followed by eight
 // imaginary components. Eight single-precison floating numbers can be loaded
 // into an AVX register.
-template <typename ParallelFor>
-struct StateSpaceAVX : public StateSpace<ParallelFor, float> {
-  using Base = StateSpace<ParallelFor, float>;
+template <typename For>
+struct StateSpaceAVX : public StateSpace<For, float> {
+  using Base = StateSpace<For, float>;
   using State = typename Base::State;
   using fp_type = typename Base::fp_type;
 
@@ -52,7 +52,7 @@ struct StateSpaceAVX : public StateSpace<ParallelFor, float> {
       _mm256_store_ps(state.get() + 16 * i + 8, val0);
     };
 
-    ParallelFor::Run(Base::num_threads_, Base::raw_size_ / 16, f, val0, state);
+    For::Run(Base::num_threads_, Base::raw_size_ / 16, f, val0, state);
   }
 
   // Uniform superposition.
@@ -82,8 +82,7 @@ struct StateSpaceAVX : public StateSpace<ParallelFor, float> {
       _mm256_store_ps(state.get() + 16 * i + 8, val0);
     };
 
-    ParallelFor::Run(Base::num_threads_, Base::raw_size_ / 16, f,
-                     val0, valu, state);
+    For::Run(Base::num_threads_, Base::raw_size_ / 16, f, val0, valu, state);
   }
 
   // |0> state.
@@ -125,7 +124,7 @@ struct StateSpaceAVX : public StateSpace<ParallelFor, float> {
            + buffer[4] + buffer[5] + buffer[6] + buffer[7];
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::raw_size_ / 16, f, Op(), state);
   }
 
@@ -156,7 +155,7 @@ struct StateSpaceAVX : public StateSpace<ParallelFor, float> {
       return std::complex<double>{re, im};
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::raw_size_ / 16, f, Op(), state1, state2);
   }
 
@@ -179,7 +178,7 @@ struct StateSpaceAVX : public StateSpace<ParallelFor, float> {
           + bre[4] + bre[5] + bre[6] + bre[7];
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::raw_size_ / 16, f, Op(), state1, state2);
   }
 
diff --git a/lib/statespace_basic.h b/lib/statespace_basic.h
index ff07c0d2..31663025 100644
--- a/lib/statespace_basic.h
+++ b/lib/statespace_basic.h
@@ -28,9 +28,9 @@ namespace qsim {
 // Routines for state-vector manipulations.
 // State is a non-vectorized sequence of one real amplitude followed by
 // one imaginary amplitude.
-template <typename ParallelFor, typename FP>
-struct StateSpaceBasic : public StateSpace<ParallelFor, FP> {
-  using Base = StateSpace<ParallelFor, FP>;
+template <typename For, typename FP>
+struct StateSpaceBasic : public StateSpace<For, FP> {
+  using Base = StateSpace<For, FP>;
   using State = typename Base::State;
   using fp_type = typename Base::fp_type;
 
@@ -43,7 +43,7 @@ struct StateSpaceBasic : public StateSpace<ParallelFor, FP> {
       state.get()[2 * i + 1] = 0;
     };
 
-    ParallelFor::Run(Base::num_threads_, Base::size_, f, state);
+    For::Run(Base::num_threads_, Base::size_, f, state);
   }
 
   // Uniform superposition.
@@ -56,8 +56,7 @@ struct StateSpaceBasic : public StateSpace<ParallelFor, FP> {
       state.get()[2 * i + 1] = 0;
     };
 
-    ParallelFor::Run(
-        Base::num_threads_, Base::size_, f, val, state);
+    For::Run(Base::num_threads_, Base::size_, f, val, state);
   }
 
   // |0> state.
@@ -93,8 +92,7 @@ struct StateSpaceBasic : public StateSpace<ParallelFor, FP> {
       return s[0] * s[0] + s[1] * s[1];
     };
 
-    return ParallelFor::RunReduce(
-        Base::num_threads_, Base::size_, f, Op(), state);
+    return For::RunReduce(Base::num_threads_, Base::size_, f, Op(), state);
   }
 
   std::complex<double> InnerProduct(
@@ -112,7 +110,7 @@ struct StateSpaceBasic : public StateSpace<ParallelFor, FP> {
       return std::complex<double>{re, im};
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::size_, f, Op(), state1, state2);
   }
 
@@ -127,7 +125,7 @@ struct StateSpaceBasic : public StateSpace<ParallelFor, FP> {
       return s1[0] * s2[0] + s1[1] * s2[1];
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::size_, f, Op(), state1, state2);
   }
 
diff --git a/lib/statespace_sse.h b/lib/statespace_sse.h
index 544f440f..14f6a4de 100644
--- a/lib/statespace_sse.h
+++ b/lib/statespace_sse.h
@@ -32,9 +32,9 @@ namespace qsim {
 // State is a vectorized sequence of four real components followed by four
 // imaginary components. Four single-precison floating numbers can be loaded
 // into an SSE register.
-template <typename ParallelFor>
-struct StateSpaceSSE : public StateSpace<ParallelFor, float> {
-  using Base = StateSpace<ParallelFor, float>;
+template <typename For>
+struct StateSpaceSSE : public StateSpace<For, float> {
+  using Base = StateSpace<For, float>;
   using State = typename Base::State;
   using fp_type = typename Base::fp_type;
 
@@ -52,7 +52,7 @@ struct StateSpaceSSE : public StateSpace<ParallelFor, float> {
       _mm_store_ps(state.get() + 8 * i + 4, val0);
     };
 
-    ParallelFor::Run(Base::num_threads_, Base::raw_size_ / 8, f, val0, state);
+    For::Run(Base::num_threads_, Base::raw_size_ / 8, f, val0, state);
   }
 
   // Uniform superposition.
@@ -76,8 +76,7 @@ struct StateSpaceSSE : public StateSpace<ParallelFor, float> {
       _mm_store_ps(state.get() + 8 * i + 4, val0);
     };
 
-    ParallelFor::Run(Base::num_threads_, Base::raw_size_ / 8, f,
-                     val0, valu, state);
+    For::Run(Base::num_threads_, Base::raw_size_ / 8, f, val0, valu, state);
   }
 
   // |0> state.
@@ -118,7 +117,7 @@ struct StateSpaceSSE : public StateSpace<ParallelFor, float> {
       return buffer[0] + buffer[1] + buffer[2] + buffer[3];
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::raw_size_ / 8, f, Op(), state);
   }
 
@@ -147,7 +146,7 @@ struct StateSpaceSSE : public StateSpace<ParallelFor, float> {
       return std::complex<double>{re, im};
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::raw_size_ / 8, f, Op(), state1, state2);
   }
 
@@ -169,7 +168,7 @@ struct StateSpaceSSE : public StateSpace<ParallelFor, float> {
       return bre[0] + bre[1] + bre[2] + bre[3];
     };
 
-    return ParallelFor::RunReduce(
+    return For::RunReduce(
         Base::num_threads_, Base::raw_size_ / 8, f, Op(), state1, state2);
   }
 
diff --git a/pybind_interface/pybind_main.cpp b/pybind_interface/pybind_main.cpp
index e3238f2e..ef26be64 100644
--- a/pybind_interface/pybind_main.cpp
+++ b/pybind_interface/pybind_main.cpp
@@ -24,10 +24,10 @@
 
 #include "../lib/bitstring.h"
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsim.h"
 #include "../lib/run_qsimh.h"
 #include "../lib/simmux.h"
@@ -128,7 +128,7 @@ std::vector<std::complex<float>> qsim_simulate(const py::dict &options) {
     return {};
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
 
@@ -175,7 +175,7 @@ py::array_t<float> qsim_simulate_fullstate(const py::dict &options) {
     return {};
   }
 
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
   using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
@@ -213,7 +213,7 @@ py::array_t<float> qsim_simulate_fullstate(const py::dict &options) {
     fsv[5] = fsv[10];
     fsv[7] = fsv[11];
   } else {
-    ParallelFor::Run(param.num_threads, buff_size / 16, reorder_fsv, fsv);
+    For::Run(param.num_threads, buff_size / 16, reorder_fsv, fsv);
   }
 
   auto capsule = py::capsule(
@@ -222,9 +222,9 @@ py::array_t<float> qsim_simulate_fullstate(const py::dict &options) {
 }
 
 std::vector<std::complex<float>> qsimh_simulate(const py::dict &options) {
-  using Simulator = qsim::Simulator<ParallelFor>;
+  using Simulator = qsim::Simulator<For>;
   using HybridSimulator = HybridSimulator<IO, GateQSim<float>, BasicGateFuser,
-                                          Simulator, ParallelFor>;
+                                          Simulator, For>;
   using Runner = QSimHRunner<IO, HybridSimulator>;
 
   Circuit<GateQSim<float>> circuit;
diff --git a/tests/BUILD b/tests/BUILD
index c4a059d6..846dd3c3 100644
--- a/tests/BUILD
+++ b/tests/BUILD
@@ -87,11 +87,11 @@ cc_test(
     deps = [
         "@com_google_googletest//:gtest_main",
         "//lib:circuit_qsim_parser",
+        "//lib:formux",
         "//lib:fuser_basic",
         "//lib:gates_qsim",
         "//lib:hybrid",
         "//lib:io",
-        "//lib:parfor",
         "//lib:simulator_avx",
     ],
     copts = select({
@@ -166,7 +166,7 @@ cc_test(
     deps = [
         ":simulator_testfixture",
         "@com_google_googletest//:gtest_main",
-        "//lib:parfor",
+        "//lib:formux",
         "//lib:simulator_avx",
     ],
     copts = select({
@@ -182,7 +182,7 @@ cc_test(
     deps = [
         ":simulator_testfixture",
         "@com_google_googletest//:gtest_main",
-        "//lib:parfor",
+        "//lib:formux",
         "//lib:simulator_basic",
     ],
     copts = select({
@@ -198,7 +198,7 @@ cc_test(
     deps = [
         ":simulator_testfixture",
         "@com_google_googletest//:gtest_main",
-        "//lib:parfor",
+        "//lib:formux",
         "//lib:simulator_sse",
     ],
     copts = select({
@@ -232,7 +232,7 @@ cc_test(
     deps = [
         ":statespace_testfixture",
         "@com_google_googletest//:gtest_main",
-        "//lib:parfor",
+        "//lib:formux",
         "//lib:simulator_avx",
         "//lib:statespace_avx",
     ],
@@ -249,7 +249,7 @@ cc_test(
     deps = [
         ":statespace_testfixture",
         "@com_google_googletest//:gtest_main",
-        "//lib:parfor",
+        "//lib:formux",
         "//lib:simulator_basic",
         "//lib:statespace_basic",
     ],
@@ -266,7 +266,7 @@ cc_test(
     deps = [
         ":statespace_testfixture",
         "@com_google_googletest//:gtest_main",
-        "//lib:parfor",
+        "//lib:formux",
         "//lib:simulator_sse",
         "//lib:statespace_sse",
     ],
diff --git a/tests/hybrid_test.cc b/tests/hybrid_test.cc
index ee589db2..e7fac2e7 100644
--- a/tests/hybrid_test.cc
+++ b/tests/hybrid_test.cc
@@ -20,11 +20,11 @@
 #include "gtest/gtest.h"
 
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/hybrid.h"
 #include "../lib/io.h"
-#include "../lib/parfor.h"
 #include "../lib/simulator_avx.h"
 
 namespace qsim {
@@ -66,9 +66,9 @@ R"(2
   EXPECT_EQ(circuit.num_qubits, 2);
   EXPECT_EQ(circuit.gates.size(), 23);
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using HybridSimulator = HybridSimulator<IO, GateQSim<float>, BasicGateFuser,
-                                          Simulator, ParallelFor>;
+                                          Simulator, For>;
   using Fuser = HybridSimulator::Fuser;
 
   std::vector<unsigned> parts = {0, 1};
@@ -245,9 +245,9 @@ R"(4
   EXPECT_EQ(circuit.num_qubits, 4);
   EXPECT_EQ(circuit.gates.size(), 63);
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using HybridSimulator = HybridSimulator<IO, GateQSim<float>, BasicGateFuser,
-                                          Simulator, ParallelFor>;
+                                          Simulator, For>;
   using Fuser = HybridSimulator::Fuser;
 
   std::vector<unsigned> parts = {0, 0, 1, 1};
diff --git a/tests/run_qsim_test.cc b/tests/run_qsim_test.cc
index 1ce7ca37..104c4ef7 100644
--- a/tests/run_qsim_test.cc
+++ b/tests/run_qsim_test.cc
@@ -22,10 +22,10 @@
 #include "gtest/gtest.h"
 
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsim.h"
 #include "../lib/simulator_avx.h"
 
@@ -72,7 +72,7 @@ TEST(RunQSimTest, QSimRunner1) {
   EXPECT_EQ(circuit.num_qubits, 4);
   EXPECT_EQ(circuit.gates.size(), 27);
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
   using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
@@ -109,7 +109,7 @@ TEST(RunQSimTest, QSimRunner2) {
   EXPECT_EQ(circuit.num_qubits, 4);
   EXPECT_EQ(circuit.gates.size(), 27);
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
   using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
@@ -144,7 +144,7 @@ TEST(RunQSimTest, CirqGates) {
   auto circuit = CirqCircuit1::GetCircuit<float>();
   const auto& expected_results = CirqCircuit1::expected_results;
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
   using Runner = QSimRunner<IO, BasicGateFuser<Cirq::GateCirq<float>>,
diff --git a/tests/run_qsimh_test.cc b/tests/run_qsimh_test.cc
index 4c0ec2a1..29d5aa29 100644
--- a/tests/run_qsimh_test.cc
+++ b/tests/run_qsimh_test.cc
@@ -22,10 +22,10 @@
 #include "gtest/gtest.h"
 
 #include "../lib/circuit_qsim_parser.h"
+#include "../lib/formux.h"
 #include "../lib/fuser_basic.h"
 #include "../lib/gates_qsim.h"
 #include "../lib/io.h"
-#include "../lib/parfor.h"
 #include "../lib/run_qsimh.h"
 #include "../lib/simulator_avx.h"
 
@@ -108,9 +108,9 @@ TEST(RunQSimHTest, QSimHRunner) {
   EXPECT_EQ(circuit.num_qubits, 4);
   EXPECT_EQ(circuit.gates.size(), 63);
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using HybridSimulator = HybridSimulator<IO, GateQSim<float>, BasicGateFuser,
-                                          Simulator, ParallelFor>;
+                                          Simulator, For>;
   using Runner = QSimHRunner<IO, HybridSimulator>;
 
   Runner::Parameter param;
@@ -165,10 +165,9 @@ TEST(RunQSimHTest, CirqGates) {
   auto circuit = CirqCircuit1::GetCircuit<float>();
   const auto& expected_results = CirqCircuit1::expected_results;
 
-  using Simulator = SimulatorAVX<ParallelFor>;
+  using Simulator = SimulatorAVX<For>;
   using HybridSimulator = HybridSimulator<IO, Cirq::GateCirq<float>,
-                                          BasicGateFuser, Simulator,
-                                          ParallelFor>;
+                                          BasicGateFuser, Simulator, For>;
   using Runner = QSimHRunner<IO, HybridSimulator>;
 
   Runner::Parameter param;
diff --git a/tests/simulator_avx_test.cc b/tests/simulator_avx_test.cc
index ba1787f3..d7843207 100644
--- a/tests/simulator_avx_test.cc
+++ b/tests/simulator_avx_test.cc
@@ -16,25 +16,25 @@
 
 #include "gtest/gtest.h"
 
-#include "../lib/parfor.h"
+#include "../lib/formux.h"
 #include "../lib/simulator_avx.h"
 
 namespace qsim {
 
 TEST(SimulatorAVXTest, ApplyGate1) {
-  TestApplyGate1<SimulatorAVX<ParallelFor>>();
+  TestApplyGate1<SimulatorAVX<For>>();
 }
 
 TEST(SimulatorAVXTest, ApplyGate2) {
-  TestApplyGate2<SimulatorAVX<ParallelFor>>();
+  TestApplyGate2<SimulatorAVX<For>>();
 }
 
 TEST(SimulatorAVXTest, ApplyGate3) {
-  TestApplyGate3<SimulatorAVX<ParallelFor>>();
+  TestApplyGate3<SimulatorAVX<For>>();
 }
 
 TEST(SimulatorAVXTest, ApplyGate5) {
-  TestApplyGate5<SimulatorAVX<ParallelFor>>();
+  TestApplyGate5<SimulatorAVX<For>>();
 }
 
 }  // namespace qsim
diff --git a/tests/simulator_basic_test.cc b/tests/simulator_basic_test.cc
index 32596a89..4347f58d 100644
--- a/tests/simulator_basic_test.cc
+++ b/tests/simulator_basic_test.cc
@@ -16,25 +16,25 @@
 
 #include "gtest/gtest.h"
 
-#include "../lib/parfor.h"
+#include "../lib/formux.h"
 #include "../lib/simulator_basic.h"
 
 namespace qsim {
 
 TEST(SimulatorBasicTest, ApplyGate1) {
-  TestApplyGate1<SimulatorBasic<ParallelFor, double>>();
+  TestApplyGate1<SimulatorBasic<For, double>>();
 }
 
 TEST(SimulatorBasicTest, ApplyGate2) {
-  TestApplyGate2<SimulatorBasic<ParallelFor, double>>();
+  TestApplyGate2<SimulatorBasic<For, double>>();
 }
 
 TEST(SimulatorBasicTest, ApplyGate3) {
-  TestApplyGate3<SimulatorBasic<ParallelFor, double>>();
+  TestApplyGate3<SimulatorBasic<For, double>>();
 }
 
 TEST(SimulatorBasicTest, ApplyGate5) {
-  TestApplyGate5<SimulatorBasic<ParallelFor, double>>();
+  TestApplyGate5<SimulatorBasic<For, double>>();
 }
 
 }  // namespace qsim
diff --git a/tests/simulator_sse_test.cc b/tests/simulator_sse_test.cc
index fbeaa7dc..5ff59f80 100644
--- a/tests/simulator_sse_test.cc
+++ b/tests/simulator_sse_test.cc
@@ -16,25 +16,25 @@
 
 #include "gtest/gtest.h"
 
-#include "../lib/parfor.h"
+#include "../lib/formux.h"
 #include "../lib/simulator_sse.h"
 
 namespace qsim {
 
 TEST(SimulatorSSETest, ApplyGate1) {
-  TestApplyGate1<SimulatorSSE<ParallelFor>>();
+  TestApplyGate1<SimulatorSSE<For>>();
 }
 
 TEST(SimulatorSSETest, ApplyGate2) {
-  TestApplyGate2<SimulatorSSE<ParallelFor>>();
+  TestApplyGate2<SimulatorSSE<For>>();
 }
 
 TEST(SimulatorSSETest, ApplyGate3) {
-  TestApplyGate3<SimulatorSSE<ParallelFor>>();
+  TestApplyGate3<SimulatorSSE<For>>();
 }
 
 TEST(SimulatorSSETest, ApplyGate5) {
-  TestApplyGate5<SimulatorSSE<ParallelFor>>();
+  TestApplyGate5<SimulatorSSE<For>>();
 }
 
 }  // namespace qsim
diff --git a/tests/statespace_avx_test.cc b/tests/statespace_avx_test.cc
index e05372d4..41af0d93 100644
--- a/tests/statespace_avx_test.cc
+++ b/tests/statespace_avx_test.cc
@@ -16,30 +16,30 @@
 
 #include "gtest/gtest.h"
 
-#include "../lib/parfor.h"
+#include "../lib/formux.h"
 #include "../lib/simulator_avx.h"
 #include "../lib/statespace_avx.h"
 
 namespace qsim {
 
 TEST(StateSpaceAVXTest, NormSmall) {
-  TestNormSmall<StateSpaceAVX<ParallelFor>>();
+  TestNormSmall<StateSpaceAVX<For>>();
 }
 
 TEST(StateSpaceAVXTest, NormAndInnerProductSmall) {
-  TestNormAndInnerProductSmall<StateSpaceAVX<ParallelFor>>();
+  TestNormAndInnerProductSmall<StateSpaceAVX<For>>();
 }
 
 TEST(StateSpaceAVXTest, NormAndInnerProduct) {
-  TestNormAndInnerProduct<SimulatorAVX<ParallelFor>>();
+  TestNormAndInnerProduct<SimulatorAVX<For>>();
 }
 
 TEST(StateSpaceAVXTest, SamplingSmall) {
-  TestSamplingSmall<StateSpaceAVX<ParallelFor>>();
+  TestSamplingSmall<StateSpaceAVX<For>>();
 }
 
 TEST(StateSpaceAVXTest, SamplingCrossEntropyDifference) {
-  TestSamplingCrossEntropyDifference<SimulatorAVX<ParallelFor>>();
+  TestSamplingCrossEntropyDifference<SimulatorAVX<For>>();
 }
 
 }  // namespace qsim
diff --git a/tests/statespace_basic_test.cc b/tests/statespace_basic_test.cc
index fa2ca702..6579d230 100644
--- a/tests/statespace_basic_test.cc
+++ b/tests/statespace_basic_test.cc
@@ -16,30 +16,30 @@
 
 #include "gtest/gtest.h"
 
-#include "../lib/parfor.h"
+#include "../lib/formux.h"
 #include "../lib/simulator_basic.h"
 #include "../lib/statespace_basic.h"
 
 namespace qsim {
 
 TEST(StateSpaceBasicTest, NormSmall) {
-  TestNormSmall<StateSpaceBasic<ParallelFor, float>>();
+  TestNormSmall<StateSpaceBasic<For, float>>();
 }
 
 TEST(StateSpaceBasicTest, NormAndInnerProductSmall) {
-  TestNormAndInnerProductSmall<StateSpaceBasic<ParallelFor, float>>();
+  TestNormAndInnerProductSmall<StateSpaceBasic<For, float>>();
 }
 
 TEST(StateSpaceBasicTest, NormAndInnerProduct) {
-  TestNormAndInnerProduct<SimulatorBasic<ParallelFor, float>>();
+  TestNormAndInnerProduct<SimulatorBasic<For, float>>();
 }
 
 TEST(StateSpaceBasicTest, SamplingSmall) {
-  TestSamplingSmall<StateSpaceBasic<ParallelFor, float>>();
+  TestSamplingSmall<StateSpaceBasic<For, float>>();
 }
 
 TEST(StateSpaceBasicTest, SamplingCrossEntropyDifference) {
-  TestSamplingCrossEntropyDifference<SimulatorBasic<ParallelFor, float>>();
+  TestSamplingCrossEntropyDifference<SimulatorBasic<For, float>>();
 }
 
 }  // namespace qsim
diff --git a/tests/statespace_sse_test.cc b/tests/statespace_sse_test.cc
index defb99dd..560df953 100644
--- a/tests/statespace_sse_test.cc
+++ b/tests/statespace_sse_test.cc
@@ -16,30 +16,30 @@
 
 #include "gtest/gtest.h"
 
-#include "../lib/parfor.h"
+#include "../lib/formux.h"
 #include "../lib/simulator_sse.h"
 #include "../lib/statespace_sse.h"
 
 namespace qsim {
 
 TEST(StateSpaceSSETest, NormSmall) {
-  TestNormSmall<StateSpaceSSE<ParallelFor>>();
+  TestNormSmall<StateSpaceSSE<For>>();
 }
 
 TEST(StateSpaceSSETest, NormAndInnerProductSmall) {
-  TestNormAndInnerProductSmall<StateSpaceSSE<ParallelFor>>();
+  TestNormAndInnerProductSmall<StateSpaceSSE<For>>();
 }
 
 TEST(StateSpaceSSETest, NormAndInnerProduct) {
-  TestNormAndInnerProduct<SimulatorSSE<ParallelFor>>();
+  TestNormAndInnerProduct<SimulatorSSE<For>>();
 }
 
 TEST(StateSpaceSSETest, SamplingSmall) {
-  TestSamplingSmall<StateSpaceSSE<ParallelFor>>();
+  TestSamplingSmall<StateSpaceSSE<For>>();
 }
 
 TEST(StateSpaceSSETest, SamplingCrossEntropyDifference) {
-  TestSamplingCrossEntropyDifference<SimulatorSSE<ParallelFor>>();
+  TestSamplingCrossEntropyDifference<SimulatorSSE<For>>();
 }
 
 }  // namespace qsim