diff --git a/.gitmodules b/.gitmodules
index fceb099..ba4ca7f 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -2,6 +2,9 @@
 	path = testscpp/Catch2
 	url = https://github.com/catchorg/Catch2.git
 	branch = v2.x
+[submodule "extern/taskflow"]
+	path = extern/taskflow
+	url = https://github.com/taskflow/taskflow.git
 [submodule "extern/spdlog"]
 	path = extern/spdlog
 	url = https://github.com/gabime/spdlog.git
diff --git a/dwave/preprocessing/include/dwave/presolve.h b/dwave/preprocessing/include/dwave/presolve.h
index 08f4c7b..1f0dd85 100644
--- a/dwave/preprocessing/include/dwave/presolve.h
+++ b/dwave/preprocessing/include/dwave/presolve.h
@@ -15,12 +15,15 @@
 #pragma once
 
 #include <algorithm>
+#include <cstdint>
 #include <unordered_map>
 #include <utility>
 #include <vector>
 
 #include "spdlog/spdlog.h"
 #include "dimod/constrained_quadratic_model.h"
+#include "taskflow/core/taskflow.hpp"
+#include "taskflow/taskflow.hpp"
 
 namespace dwave {
 namespace presolve {
@@ -115,6 +118,23 @@ void Postsolver<bias_type, index_type, assignment_type>::substitute_variable(ind
     transforms_.back().offset = offset;
 }
 
+class PresolverTechniques {
+ public:
+    const static uint64_t SPIN_TO_BINARY = 1 << 0;
+    const static uint64_t REMOVE_OFFSETS = 1 << 1;
+    const static uint64_t FLIP_CONSTRAINTS = 1 << 2;
+    const static uint64_t REMOVE_SELF_LOOPS = 1 << 3;
+    const static uint64_t REMOVE_INVALID_MARKERS = 1 << 4;
+    const static uint64_t CHECK_FOR_NAN = 1 << 5;
+    const static uint64_t REMOVE_SINGLE_VAR_CONSTRAINTS = 1 << 6;
+    const static uint64_t REMOVE_ZERO_BIASES = 1 << 7;
+    const static uint64_t TIGHTEN_BOUNDS = 1 << 8;
+    const static uint64_t REMOVE_FIXED_VARS = 1 << 9;
+
+    const static uint64_t ALL = -1;
+    const static uint64_t NONE = 0;
+};
+
 template <class Bias, class Index, class Assignment>
 class Presolver {
  public:
@@ -140,7 +160,10 @@ class Presolver {
     model_type detach_model();
 
     /// Load the default presolve techniques.
-    void load_default_presolvers();
+    void load_default_presolvers(uint32_t max_rounds = 100) {
+        load_presolvers(PresolverTechniques::ALL, max_rounds);
+    }
+    void load_presolvers(uint64_t presolvers = PresolverTechniques::ALL, uint32_t max_rounds = 100);
 
     /// Return a const reference to the held constrained quadratic model.
     const model_type& model() const;
@@ -149,11 +172,28 @@ class Presolver {
     const Postsolver<bias_type, index_type, assignment_type>& postsolver() const;
 
  private:
+    struct TfHelper {
+        tf::Executor executor;
+        tf::Taskflow taskflow_onetime;
+        tf::Taskflow taskflow_trivial;
+        tf::Taskflow taskflow_cleanup;
+        uint32_t loop_counter;
+        bool loop_changed;
+        bool model_feasible = true;
+
+        bool operator=(const struct TfHelper& that) { return true; }
+    };
+
+    struct TfHelper tf_helper_;
+
+    void load_taskflow_onetime();
+    void load_taskflow_trivial(int max_rounds = 100);
+    void load_taskflow_cleanup();
+
     model_type model_;
     Postsolver<bias_type, index_type, assignment_type> postsolver_;
 
-    // todo: replace this with a vector of pointers or similar
-    bool default_techniques_;
+    uint64_t presolvers_ = PresolverTechniques::NONE;
 
     bool detached_;
 
@@ -271,6 +311,10 @@ class Presolver {
         return false;
     }
     bool technique_remove_single_variable_constraints() {
+        if (!tf_helper_.model_feasible) {
+            return false;
+        }
+
         bool ret = false;
         size_type c = 0;
         while (c < model_.num_constraints()) {
@@ -282,20 +326,17 @@ class Presolver {
                     switch (constraint.sense()) {
                         case dimod::Sense::EQ:
                             if (constraint.offset() != constraint.rhs()) {
-                                // need this exact message for Python
-                                throw std::logic_error("infeasible");
+                                tf_helper_.model_feasible = false;
                             }
                             break;
                         case dimod::Sense::LE:
                             if (constraint.offset() > constraint.rhs()) {
-                                // need this exact message for Python
-                                throw std::logic_error("infeasible");
+                                tf_helper_.model_feasible = false;
                             }
                             break;
                         case dimod::Sense::GE:
                             if (constraint.offset() < constraint.rhs()) {
-                                // need this exact message for Python
-                                throw std::logic_error("infeasible");
+                                tf_helper_.model_feasible = false;
                             }
                             break;
                     }
@@ -341,16 +382,22 @@ class Presolver {
         return ret;
     }
     bool technique_remove_zero_biases() {
-        bool ret = false;
+        if (!tf_helper_.model_feasible) {
+            return false;
+        }
 
+        bool ret = false;
         ret |= remove_zero_biases(model_.objective);
         for (size_t c = 0; c < model_.num_constraints(); ++c) {
             ret |= remove_zero_biases(model_.constraint_ref(c));
         }
-
         return ret;
     }
     bool technique_tighten_bounds() {
+        if (!tf_helper_.model_feasible) {
+            return false;
+        }
+
         bool ret = false;
         bias_type lb;
         bias_type ub;
@@ -377,7 +424,11 @@ class Presolver {
         return ret;
     }
     bool technique_remove_fixed_variables() {
-        bool ret = false; 
+        if (!tf_helper_.model_feasible) {
+            return false;
+        }
+
+        bool ret = false;
         size_type v = 0;
         while (v < model_.num_variables()) {
             if (model_.lower_bound(v) == model_.upper_bound(v)) {
@@ -419,54 +470,25 @@ class Presolver {
 
 template <class bias_type, class index_type, class assignment_type>
 Presolver<bias_type, index_type, assignment_type>::Presolver()
-        : model_(), postsolver_(), default_techniques_(false), detached_(false) {}
+        : model_(), postsolver_(), detached_(false) {}
 
 template <class bias_type, class index_type, class assignment_type>
 Presolver<bias_type, index_type, assignment_type>::Presolver(model_type model)
-        : model_(std::move(model)), postsolver_(), default_techniques_(), detached_(false) {}
+        : model_(std::move(model)), postsolver_(), detached_(false) {}
 
 template <class bias_type, class index_type, class assignment_type>
 void Presolver<bias_type, index_type, assignment_type>::apply() {
     if (detached_) throw std::logic_error("model has been detached, presolver is no longer valid");
-
-    // If no techniques have been loaded, return early.
-    if (!default_techniques_) return;
-
-    // One time techniques ----------------------------------------------------
-
-    // *-- spin-to-binary
-    technique_spin_to_binary();
-    // *-- remove offsets
-    technique_remove_offsets();
-    // *-- flip >= constraints
-    technique_flip_constraints();
-    // *-- remove self-loops
-    technique_remove_self_loops();
-
-    // Trivial techniques -----------------------------------------------------
-
-    bool changes = true;
-    const index_type max_num_rounds = 100;  // todo: make configurable
-    for (index_type num_rounds = 0; num_rounds < max_num_rounds; ++num_rounds) {
-        if (!changes) break;
-        changes = false;
-
-        // *-- clear out 0 variables/interactions in the constraints and objective
-        changes |= technique_remove_zero_biases();
-        // *-- todo: check for NAN
-        changes |= technique_check_for_nan();
-        // *-- remove single variable constraints
-        changes |= technique_remove_single_variable_constraints();
-        // *-- tighten bounds based on vartype
-        changes |= technique_tighten_bounds();
-        // *-- remove variables that are fixed by bounds
-        changes |= technique_remove_fixed_variables();
-   }
-
-    // Cleanup
-
-    // *-- remove any invalid discrete markers
-    technique_remove_invalid_markers();
+    if (presolvers_ == PresolverTechniques::NONE) return;
+
+    tf_helper_.executor.run(tf_helper_.taskflow_onetime).wait();
+    tf_helper_.executor.run(tf_helper_.taskflow_trivial).wait();
+    if (tf_helper_.model_feasible) {
+        tf_helper_.executor.run(tf_helper_.taskflow_cleanup).wait();
+    } else {
+        // need this exact message for Python
+        throw std::logic_error("infeasible");
+    }
 }
 
 template <class bias_type, class index_type, class assignment_type>
@@ -483,8 +505,117 @@ Presolver<bias_type, index_type, assignment_type>::detach_model() {
 }
 
 template <class bias_type, class index_type, class assignment_type>
-void Presolver<bias_type, index_type, assignment_type>::load_default_presolvers() {
-    default_techniques_ = true;
+void Presolver<bias_type, index_type, assignment_type>::load_presolvers(uint64_t presolvers,
+                                                                        uint32_t max_rounds) {
+    presolvers_ = presolvers;
+
+    load_taskflow_onetime();
+    load_taskflow_trivial(max_rounds);
+    load_taskflow_cleanup();
+}
+
+template <class bias_type, class index_type, class assignment_type>
+void Presolver<bias_type, index_type, assignment_type>::load_taskflow_onetime() {
+    std::vector<tf::Task> tasks;
+
+    if (presolvers_ | PresolverTechniques::SPIN_TO_BINARY) {
+        auto task = tf_helper_.taskflow_onetime.emplace([&]() { technique_spin_to_binary(); });
+        task.name("spin_to_binary");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::REMOVE_OFFSETS) {
+        auto task = tf_helper_.taskflow_onetime.emplace([&]() { technique_remove_offsets(); });
+        task.name("remove_offsets");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::FLIP_CONSTRAINTS) {
+        auto task = tf_helper_.taskflow_onetime.emplace([&]() { technique_flip_constraints(); });
+        task.name("flip_constraints");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::REMOVE_SELF_LOOPS) {
+        auto task = tf_helper_.taskflow_onetime.emplace([&]() { technique_remove_self_loops(); });
+        task.name("remove_self_loops");
+        tasks.emplace_back(task);
+    }
+
+    for (auto i = 0; i < tasks.size() - 1; i++) {
+        tasks[i].precede(tasks[i + 1]);
+    }
+}
+
+template <class bias_type, class index_type, class assignment_type>
+void Presolver<bias_type, index_type, assignment_type>::load_taskflow_trivial(int max_rounds) {
+    auto alpha = tf_helper_.taskflow_trivial.emplace([&]() {
+        tf_helper_.loop_changed = false;
+        tf_helper_.loop_counter = 0;
+    });
+    alpha.name("initialize");
+
+    auto omega = tf_helper_.taskflow_trivial.emplace([&]() {
+        if (tf_helper_.model_feasible && tf_helper_.loop_changed &&
+            ++tf_helper_.loop_counter < max_rounds) {
+            tf_helper_.loop_changed = false;
+            return 0;  // This will loop back to the first task in tasks, below
+        }
+        return 1;  // This will exit the loop
+    });
+    omega.name("conditional");
+
+    std::vector<tf::Task> tasks;
+    if (presolvers_ | PresolverTechniques::REMOVE_ZERO_BIASES) {
+        auto task = tf_helper_.taskflow_trivial.emplace(
+                [&]() { tf_helper_.loop_changed |= technique_remove_zero_biases(); });
+        task.name("remove_zero_biases");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::CHECK_FOR_NAN) {
+        auto task = tf_helper_.taskflow_trivial.emplace(
+                [&]() { tf_helper_.loop_changed |= technique_check_for_nan(); });
+        task.name("check_for_nan");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::REMOVE_SINGLE_VAR_CONSTRAINTS) {
+        auto task = tf_helper_.taskflow_trivial.emplace([&]() {
+            tf_helper_.loop_changed |= technique_remove_single_variable_constraints();
+        });
+        task.name("remove_single_variable_constraints");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::TIGHTEN_BOUNDS) {
+        auto task = tf_helper_.taskflow_trivial.emplace(
+                [&]() { tf_helper_.loop_changed |= technique_tighten_bounds(); });
+        task.name("tighten_bounds");
+        tasks.emplace_back(task);
+    }
+    if (presolvers_ | PresolverTechniques::REMOVE_FIXED_VARS) {
+        auto task = tf_helper_.taskflow_trivial.emplace(
+                [&]() { tf_helper_.loop_changed |= technique_remove_fixed_variables(); });
+        task.name("remove_fixed_variables");
+        tasks.emplace_back(task);
+    }
+
+    if (tasks.size() > 0) {
+        alpha.precede(tasks[0]);
+
+        for (auto i = 0; i < tasks.size() - 1; i++) {
+            tasks[i].precede(tasks[i + 1]);
+        }
+
+        tasks[tasks.size() - 1].precede(omega);
+        omega.precede(tasks[0]);
+    } else {
+        alpha.precede(omega);
+    }
+}
+
+template <class bias_type, class index_type, class assignment_type>
+void Presolver<bias_type, index_type, assignment_type>::load_taskflow_cleanup() {
+    if (presolvers_ | PresolverTechniques::REMOVE_INVALID_MARKERS) {
+        auto task =
+                tf_helper_.taskflow_cleanup.emplace([&]() { technique_remove_invalid_markers(); });
+        task.name("remove_invalid_markers");
+    }
 }
 
 template <class bias_type, class index_type, class assignment_type>
diff --git a/extern/taskflow b/extern/taskflow
new file mode 160000
index 0000000..f1490ff
--- /dev/null
+++ b/extern/taskflow
@@ -0,0 +1 @@
+Subproject commit f1490ffc286eba418107c4aa9b8913e8ca76443c
diff --git a/setup.py b/setup.py
index be30156..3d9005e 100644
--- a/setup.py
+++ b/setup.py
@@ -26,6 +26,7 @@
     'msvc': ['/std:c++17', '/EHsc'],
     'unix': [
         '-std=c++17', 
+        '-pthread'
        ],
 }
 
@@ -62,6 +63,7 @@ def build_extensions(self):
     include_dirs=[
         numpy.get_include(),
         dimod.get_include(),
+        'extern/taskflow',
         'extern/spdlog/include/',
         ],
     install_requires=[
diff --git a/testscpp/Makefile b/testscpp/Makefile
index b22be20..0979428 100644
--- a/testscpp/Makefile
+++ b/testscpp/Makefile
@@ -2,8 +2,9 @@ ROOT := ..
 SRC := $(ROOT)/dwave/preprocessing/
 CATCH2 := $(ROOT)/testscpp/Catch2/single_include/
 DIMOD := $(shell python -c 'import dimod; print(dimod.get_include())')
+TASKFLOW := $(ROOT)/extern/taskflow/
 SPDLOG := $(ROOT)/extern/spdlog/include/
-INCLUDES := -I $(SRC)/include/ -I $(DIMOD) -I $(CATCH2) -I $(SPDLOG)
+INCLUDES := -I $(SRC)/include/ -I $(DIMOD) -I $(CATCH2) -I $(SPDLOG) -I $(TASKFLOW)
 FLAGS := -std=c++17 -Wall -Wno-unknown-pragmas -Wno-sign-compare -Wno-deprecated-declarations -fcompare-debug-second -O3 
 
 all: update test_main test_main_parallel tests tests_parallel
@@ -15,8 +16,8 @@ tests_parallel: test_main_parallel.out
 	./test_main_parallel
 
 test_main: test_main.cpp
-	g++ $(FLAGS) -c test_main.cpp
-	g++ $(FLAGS) test_main.o tests/*.cpp -o test_main $(INCLUDES)
+	g++ $(FLAGS) -pthread -c test_main.cpp
+	g++ $(FLAGS) -pthread test_main.o tests/*.cpp -o test_main $(INCLUDES)
 
 test_main_parallel: test_main.cpp
 	g++ $(FLAGS) -fopenmp -c test_main.cpp -o test_main_parallel.o