diff --git a/resolve/LinSolverDirectCuSolverGLU.cpp b/resolve/LinSolverDirectCuSolverGLU.cpp
index 977cc5a2c..5608fd7e6 100644
--- a/resolve/LinSolverDirectCuSolverGLU.cpp
+++ b/resolve/LinSolverDirectCuSolverGLU.cpp
@@ -1,8 +1,6 @@
 #include "LinSolverDirectCuSolverGLU.hpp"
 
-#include <algorithm>
 #include <cstring> // includes memcpy
-#include <vector>
 
 #include <resolve/Profiling.hpp>
 #include <resolve/matrix/Csr.hpp>
@@ -236,10 +234,6 @@ namespace ReSolve
         M_col[count++] = U_col[j];
       }
     }
-    for (index_type i = 0; i < n; ++i) // this is crucial, turns out somehow the indices are not sorted
-    {
-      std::sort(M_col + M_row[i], M_col + M_row[i + 1]);
-    }
   }
 
   void LinSolverDirectCuSolverGLU::combineFactors(matrix::Sparse* L, matrix::Sparse* U)
diff --git a/resolve/LinSolverDirectCuSolverRf.cpp b/resolve/LinSolverDirectCuSolverRf.cpp
index 83d49b8ef..5f9f267f2 100644
--- a/resolve/LinSolverDirectCuSolverRf.cpp
+++ b/resolve/LinSolverDirectCuSolverRf.cpp
@@ -1,9 +1,7 @@
 #include "LinSolverDirectCuSolverRf.hpp"
 
-#include <algorithm>
 #include <cassert>
 #include <cstring> // includes memcpy
-#include <vector>
 
 #include <resolve/matrix/Csc.hpp>
 #include <resolve/matrix/Csr.hpp>
@@ -108,16 +106,6 @@ namespace ReSolve
 
     status_cusolverrf_ = cusolverRfSetResetValuesFastMode(handle_cusolverrf_, CUSOLVERRF_RESET_VALUES_FAST_MODE_ON);
     error_sum += status_cusolverrf_;
-    // sort L and U columns
-    for (index_type i = 0; i < n; ++i)
-    {
-      std::sort(L->getColData(memory::HOST) + L->getRowData(memory::HOST)[i],
-                L->getColData(memory::HOST) + L->getRowData(memory::HOST)[i + 1]);
-      std::sort(U->getColData(memory::HOST) + U->getRowData(memory::HOST)[i],
-                U->getColData(memory::HOST) + U->getRowData(memory::HOST)[i + 1]);
-    }
-    L->setUpdated(memory::HOST);
-    U->setUpdated(memory::HOST);
     L->syncData(memory::DEVICE);
     U->syncData(memory::DEVICE);
     status_cusolverrf_ = cusolverRfSetupDevice(n,
diff --git a/resolve/LinSolverDirectKLU.cpp b/resolve/LinSolverDirectKLU.cpp
index 9764740a8..c9396a881 100644
--- a/resolve/LinSolverDirectKLU.cpp
+++ b/resolve/LinSolverDirectKLU.cpp
@@ -1,7 +1,9 @@
 #include "LinSolverDirectKLU.hpp"
 
+#include <algorithm>
 #include <cassert>
 #include <cstring>
+#include <vector>
 
 #include <resolve/matrix/Csc.hpp>
 #include <resolve/matrix/Csr.hpp>
@@ -295,6 +297,15 @@ namespace ReSolve
                            nullptr,
                            nullptr,
                            &Common_);
+      // Sort the column indices in L and U to ensure they are in ordered CSR format
+      // WARNING: Values are not sorted. We currently don't use values from KLU across solvers. If we ever decide to, we will need to change this.
+      for (index_type i = 0; i < A_->getNumRows(); ++i)
+      {
+        std::sort(L_->getColData(memory::HOST) + L_->getRowData(memory::HOST)[i], // Sort L's column indices from the start of Row i to the start of Row i+1 (non inclusive).
+                  L_->getColData(memory::HOST) + L_->getRowData(memory::HOST)[i + 1]);
+        std::sort(U_->getColData(memory::HOST) + U_->getRowData(memory::HOST)[i],
+                  U_->getColData(memory::HOST) + U_->getRowData(memory::HOST)[i + 1]);
+      }
 
       L_->setUpdated(memory::HOST);
       U_->setUpdated(memory::HOST);
diff --git a/resolve/LinSolverDirectRocSolverRf.cpp b/resolve/LinSolverDirectRocSolverRf.cpp
index 40d9a2aba..b737a7c3a 100644
--- a/resolve/LinSolverDirectRocSolverRf.cpp
+++ b/resolve/LinSolverDirectRocSolverRf.cpp
@@ -842,10 +842,6 @@ namespace ReSolve
         M_col[count++] = U_col[j];
       }
     }
-    for (index_type i = 0; i < n; ++i) // this is crucial, turns out somehow the indices are not sorted
-    {
-      std::sort(M_col + M_row[i], M_col + M_row[i + 1]);
-    }
   }
 
   /**