diff --git a/benchmarks/linear_programming/run_mps_files.sh b/benchmarks/linear_programming/run_mps_files.sh
index 20eb3af4e2..8be19a3c39 100755
--- a/benchmarks/linear_programming/run_mps_files.sh
+++ b/benchmarks/linear_programming/run_mps_files.sh
@@ -360,7 +360,7 @@ else
         mapfile -t mps_files < <(find "$MPS_DIR" -type f \( -name "*.mps" -o -name "*.MPS" -o -name "*.SIF" \) | sort)
     else
         # Gather .mps/.MPS and .SIF files in the directory
-        mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps "$MPS_DIR"/*.MPS "$MPS_DIR"/*.SIF 2>/dev/null)
+        mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps "$MPS_DIR"/*.MPS "$MPS_DIR"/*.SIF "$MPS_DIR"/*.mps.gz 2>/dev/null)
     fi
 
     echo "Found ${#mps_files[@]} .mps and .SIF files in $MPS_DIR"
diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu
index 66fa4c2c5f..5db791f7e9 100644
--- a/cpp/src/barrier/barrier.cu
+++ b/cpp/src/barrier/barrier.cu
@@ -1715,6 +1715,26 @@ class iteration_data_t {
     cusparse_view.spmv(alpha, cusparse_u, beta, cusparse_v);
   }
 
+  // v = alpha * A * Dinv * A^T * y + beta * v. Simple interface (plain device vectors,
+  // no pre-built cusparse descriptors) so it can be used as the `a_multiply` callback of
+  // the generic iterative-refinement operator for the ADAT (non-augmented) solve path.
+  void gpu_adat_multiply_simple(f_t alpha,
+                                const rmm::device_uvector<f_t>& y,
+                                f_t beta,
+                                rmm::device_uvector<f_t>& v)
+  {
+    const i_t n = A.n;
+    rmm::device_uvector<f_t> u(n, stream_view_);
+    cusparse_view_.transpose_spmv(1.0, y, 0.0, u);
+    cub::DeviceTransform::Transform(cuda::std::make_tuple(u.data(), d_inv_diag.data()),
+                                    u.data(),
+                                    u.size(),
+                                    cuda::std::multiplies<>{},
+                                    stream_view_.value());
+    RAFT_CHECK_CUDA(stream_view_);
+    cusparse_view_.spmv(alpha, u, beta, v);
+  }
+
   // v = alpha * A * Dinv * A^T * y + beta * v
   void adat_multiply(f_t alpha,
                      const dense_vector_t<i_t, f_t>& y,
@@ -2947,6 +2967,37 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
         settings.log.printf("Linear solve failed\n");
         return -1;
       }
+
+      // Iterative refinement on the ADAT (Schur-complement) system using GMRES.
+      // The direct Cholesky solve can degrade in accuracy on ill-conditioned D near
+      // convergence, as the diagonal D can span many orders of magnitude with small
+      // barrier parameter. In this case, we launch a GMRES-based iterative refinement
+      // loop for added robustness in the Schur-complement (ADAT) approach.
+      // GMRES can handle large, potentially ill-conditioned systems better than simple Richardson
+      // or classical iterative refinement, at the potential cost of higher computational work and
+      // memory. This is only used on the pure Schur-complement (n_dense_columns == 0).
+      if (settings.barrier_iterative_refinement && data.n_dense_columns == 0) {
+        struct adat_op_t {
+          adat_op_t(iteration_data_t<i_t, f_t>& data) : data_(data) {}
+          iteration_data_t<i_t, f_t>& data_;
+          void a_multiply(f_t alpha,
+                          const rmm::device_uvector<f_t>& x,
+                          f_t beta,
+                          rmm::device_uvector<f_t>& y) const
+          {
+            data_.gpu_adat_multiply_simple(alpha, x, beta, y);
+          }
+          void solve(rmm::device_uvector<f_t>& b, rmm::device_uvector<f_t>& x) const
+          {
+            data_.gpu_solve_adat(b, x);
+          }
+        } adat_op(data);
+        const f_t adat_solve_err =
+          iterative_refinement<i_t, f_t, adat_op_t>(adat_op, data.d_h_, data.d_dy_);
+        if (adat_solve_err > 1e-1) {
+          settings.log.printf("||ADAT*dy - h|| %e after IR\n", adat_solve_err);
+        }
+      }
     }  // Close NVTX range
 
     // y_residual <- ADAT*dy - h