Avoid transposing J.

nmodl/ode.py

@@ -315,6 +315,7 @@ def solve_non_lin_system(eq_strings, vars, constants, function_calls):
     vecJcode = []
     for i, j in itertools.product(range(jacobian.rows), range(jacobian.cols)):
         flat_index = i + jacobian.rows * j
+        flat_index = i*jacobian.cols + j
 
         rhs = sp.ccode(jacobian[i,j].simplify().subs(X_vec_map).evalf(), user_functions=custom_fcts)
         vecJcode.append(f"J[{flat_index}] = {rhs}")

src/solver/newton/newton.hpp

@@ -60,7 +60,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
     // Vector to store result of function F(X):
     Eigen::Matrix<double, N, 1> F;
     // Matrix to store jacobian of F(X):
-    Eigen::Matrix<double, N, N> J;
+    Eigen::Matrix<double, N, N, Eigen::RowMajor> J;
     // Solver iteration count:
     int iter = -1;
     while (++iter < max_iter) {
@@ -72,12 +72,6 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
             // we have converged: return iteration count
             return iter;
         }
-        // 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
-        // instead of the columns.
-        if (!J.IsRowMajor)
-            J.transposeInPlace();
         Eigen::Matrix<int, N, 1> pivot;
         Eigen::Matrix<double, N, 1> rowmax;
         // Check if J is singular