References for 'BlueBrain/nmodl#1349'.

BlueBrain · Jul 16, 2024 · f310173 · f310173
1 parent 237859f
commit f310173
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 3 deletions.
diff --git a/global/neuron/thread_newton.cpp b/global/neuron/thread_newton.cpp
@@ -279,6 +279,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
             // we have converged: return iteration count
             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
@@ -290,8 +293,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

diff --git a/kinetic/coreneuron/X2Y.cpp b/kinetic/coreneuron/X2Y.cpp
@@ -288,6 +288,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
             // we have converged: return iteration count
             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
@@ -299,8 +302,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

diff --git a/kinetic/neuron/X2Y.cpp b/kinetic/neuron/X2Y.cpp
@@ -279,6 +279,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
             // we have converged: return iteration count
             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
@@ -290,8 +293,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