From 1f8f3c186377800d349d6f32b4cf2c38069d1baf Mon Sep 17 00:00:00 2001
From: Luc Grosheintz <luc.grosheintz@gmail.com>
Date: Tue, 16 Jul 2024 14:44:11 +0200
Subject: [PATCH] Use Eigen on CPUs.

---
 src/solver/newton/newton.hpp | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/solver/newton/newton.hpp b/src/solver/newton/newton.hpp
index 5d4e05182..522750f1c 100644
--- a/src/solver/newton/newton.hpp
+++ b/src/solver/newton/newton.hpp
@@ -95,6 +95,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
         if (is_converged(X, J, F, eps)) {
             return iter;
         }
+        Eigen::Matrix<double, N, 1> X_solve;
+
+#if defined(CORENEURON_ENABLE_GPU)
         // In Eigen the default storage order is ColMajor.
         // Crout's implementation requires matrices stored in RowMajor order (C-style arrays).
         // Therefore, the transposeInPlace is critical such that the data() method to give the rows
@@ -107,9 +110,10 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
         // Check if J is singular
         if (nmodl::crout::Crout<double>(N, J.data(), pivot.data(), rowmax.data()) < 0) {
             return -1;
-        }
-        Eigen::Matrix<double, N, 1> X_solve;
         nmodl::crout::solveCrout<double>(N, J.data(), F.data(), X_solve.data(), pivot.data());
+#else
+        X_solve = J.partialPivLu().solve(F);
+#endif
         X -= X_solve;
     }
     // If we fail to converge after max_iter iterations, return -1