From f310173eaa124ec40f0e97f089666aceff46ff6a Mon Sep 17 00:00:00 2001 From: GitHub Actions Bot <> Date: Tue, 16 Jul 2024 13:22:33 +0000 Subject: [PATCH] References for 'https://github.com/BlueBrain/nmodl/pull/1349'. --- global/neuron/thread_newton.cpp | 7 ++++++- kinetic/coreneuron/X2Y.cpp | 7 ++++++- kinetic/neuron/X2Y.cpp | 7 ++++++- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/global/neuron/thread_newton.cpp b/global/neuron/thread_newton.cpp index 4f08bda..1d2d594 100644 --- a/global/neuron/thread_newton.cpp +++ b/global/neuron/thread_newton.cpp @@ -279,6 +279,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix& X, // we have converged: return iteration count return iter; } + Eigen::Matrix 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& X, // Check if J is singular if (nmodl::crout::Crout(N, J.data(), pivot.data(), rowmax.data()) < 0) return -1; - Eigen::Matrix X_solve; nmodl::crout::solveCrout(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 index fd20fb7..a0b8f78 100644 --- a/kinetic/coreneuron/X2Y.cpp +++ b/kinetic/coreneuron/X2Y.cpp @@ -288,6 +288,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix& X, // we have converged: return iteration count return iter; } + Eigen::Matrix 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& X, // Check if J is singular if (nmodl::crout::Crout(N, J.data(), pivot.data(), rowmax.data()) < 0) return -1; - Eigen::Matrix X_solve; nmodl::crout::solveCrout(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 index b69c155..b94e9a5 100644 --- a/kinetic/neuron/X2Y.cpp +++ b/kinetic/neuron/X2Y.cpp @@ -279,6 +279,9 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix& X, // we have converged: return iteration count return iter; } + Eigen::Matrix 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& X, // Check if J is singular if (nmodl::crout::Crout(N, J.data(), pivot.data(), rowmax.data()) < 0) return -1; - Eigen::Matrix X_solve; nmodl::crout::solveCrout(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