diff --git a/cpp/src/mip/presolve/third_party_presolve.cpp b/cpp/src/mip/presolve/third_party_presolve.cpp
index 75bad369f8..7c9867d0e7 100644
--- a/cpp/src/mip/presolve/third_party_presolve.cpp
+++ b/cpp/src/mip/presolve/third_party_presolve.cpp
@@ -260,7 +260,11 @@ void check_postsolve_status(const papilo::PostsolveStatus& status)
 {
   switch (status) {
     case papilo::PostsolveStatus::kOk: CUOPT_LOG_INFO("Post-solve status:: succeeded"); break;
-    case papilo::PostsolveStatus::kFailed: CUOPT_LOG_INFO("Post-solve status:: failed"); break;
+    case papilo::PostsolveStatus::kFailed:
+      CUOPT_LOG_INFO(
+        "Post-solve status:: Post solved solution violates constraints. This is most likely due to "
+        "different tolerances.");
+      break;
   }
 }
 
@@ -279,8 +283,8 @@ void set_presolve_methods(papilo::Presolve<f_t>& presolver, problem_category_t c
   presolver.addPresolveMethod(uptr(new papilo::SimpleProbing<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::ParallelRowDetection<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::ParallelColDetection<f_t>()));
-  // FIXME: Postsolve fails with this method
-  // presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing<f_t>()));
+
+  presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::DualFix<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::SimplifyInequalities<f_t>()));
 
@@ -303,16 +307,6 @@ void set_presolve_options(papilo::Presolve<f_t>& presolver,
                           double time_limit)
 {
   presolver.getPresolveOptions().tlim = time_limit;
-
-  if (category == problem_category_t::LP) {
-    presolver.getPresolveOptions().useabsfeas = false;
-    presolver.getPresolveOptions().feastol    = relative_tolerance;
-    presolver.getPresolveOptions().epsilon    = absolute_tolerance;
-  } else if (category == problem_category_t::MIP || category == problem_category_t::IP) {
-    presolver.getPresolveOptions().useabsfeas = true;
-    presolver.getPresolveOptions().feastol    = absolute_tolerance;
-    presolver.getPresolveOptions().epsilon    = absolute_tolerance;
-  }
 }
 
 template <typename i_t, typename f_t>
diff --git a/cpp/tests/linear_programming/pdlp_test.cu b/cpp/tests/linear_programming/pdlp_test.cu
index 47fae97c1e..2598561211 100644
--- a/cpp/tests/linear_programming/pdlp_test.cu
+++ b/cpp/tests/linear_programming/pdlp_test.cu
@@ -219,17 +219,21 @@ TEST(pdlp_class, run_sub_mittleman)
     for (auto solver_mode : solver_mode_list) {
       auto settings             = pdlp_solver_settings_t<int, double>{};
       settings.pdlp_solver_mode = solver_mode;
-      const raft::handle_t handle_{};
-      optimization_problem_solution_t<int, double> solution =
-        solve_lp(&handle_, op_problem, settings);
-      EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL);
-      EXPECT_FALSE(
-        is_incorrect_objective(expected_objective_value,
-                               solution.get_additional_termination_information().primal_objective));
-      test_objective_sanity(op_problem,
-                            solution.get_primal_solution(),
-                            solution.get_additional_termination_information().primal_objective);
-      test_constraint_sanity(op_problem, solution);
+      for (auto [presolve, epsilon] : {std::pair{true, 1e-1}, std::pair{false, 1e-6}}) {
+        settings.presolve = presolve;
+        const raft::handle_t handle_{};
+        optimization_problem_solution_t<int, double> solution =
+          solve_lp(&handle_, op_problem, settings);
+        EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL);
+        EXPECT_FALSE(is_incorrect_objective(
+          expected_objective_value,
+          solution.get_additional_termination_information().primal_objective));
+        test_objective_sanity(op_problem,
+                              solution.get_primal_solution(),
+                              solution.get_additional_termination_information().primal_objective,
+                              epsilon);
+        test_constraint_sanity(op_problem, solution, epsilon, presolve);
+      }
     }
   }
 }
diff --git a/cpp/tests/linear_programming/utilities/pdlp_test_utilities.cuh b/cpp/tests/linear_programming/utilities/pdlp_test_utilities.cuh
index 7a624e0232..e0144b75ba 100644
--- a/cpp/tests/linear_programming/utilities/pdlp_test_utilities.cuh
+++ b/cpp/tests/linear_programming/utilities/pdlp_test_utilities.cuh
@@ -70,7 +70,7 @@ static void test_constraint_sanity(
   const cuopt::mps_parser::mps_data_model_t<int, double>& op_problem,
   const optimization_problem_solution_t<int, double>& solution,
   double epsilon        = tolerance,
-  bool presolve_enabled = true)
+  bool presolve_enabled = false)
 {
   const std::vector<double> primal_vars              = host_copy(solution.get_primal_solution());
   const std::vector<double>& values                  = op_problem.get_constraint_matrix_values();
@@ -83,6 +83,7 @@ static void test_constraint_sanity(
   std::vector<double> residual(solution.get_dual_solution().size(), 0.0);
   std::vector<double> viol(solution.get_dual_solution().size(), 0.0);
 
+  // No dual solution and residual for presolved problems
   if (!presolve_enabled) {
     // CSR SpMV
     for (size_t i = 0; i < offsets.size() - 1; ++i) {
@@ -137,8 +138,14 @@ static void test_constraint_sanity(
   for (size_t i = 0; i < primal_vars.size(); ++i) {
     // Not always stricly true because we apply variable bound clamping on the scaled problem
     // After unscaling it, the variables might not respect exactly (this adding an epsilon)
-    EXPECT_TRUE(primal_vars[i] >= variable_lower_bounds[i] - epsilon &&
-                primal_vars[i] <= variable_upper_bounds[i] + epsilon);
+    auto condition = primal_vars[i] >= variable_lower_bounds[i] - epsilon &&
+                     primal_vars[i] <= variable_upper_bounds[i] + epsilon;
+    if (!condition) {
+      std::cout << "Variable " << i << " is " << primal_vars[i] << " but should be between "
+                << variable_lower_bounds[i] - epsilon << " and "
+                << variable_upper_bounds[i] + epsilon << std::endl;
+    }
+    EXPECT_TRUE(condition);
   }
 }