Skip to content
This repository was archived by the owner on Mar 25, 2025. It is now read-only.

Use Eigen on CPUs. #1349

Open
1uc wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Use Eigen on CPUs. #1349

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions src/solver/newton/newton.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

@pramodk pramodk Jul 17, 2024

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add a brief description/comment if this is done for performance reasons or something else?

I don't remember this well and hence asking for clarity: do we expect small differences in the results between the two? If there are two code paths, just think from CPU vs GPU debugging or results comparison perspectives. (obviously, not thinking fp differences as they exist on cpu vs gpu. But more from solver/accuracy perspective).

@1uc 1uc Jul 17, 2024

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, without using NMODL is a 25% regression, compared to NOCMODL on the Masoli circuit we've been running recently.

That's a tricky question. The variation we implement ourselves is also "partial pivoting". There's two ways of doing LU decompositions: Dolittle and Crout. They result in almost the same matrices. Mostly I think they just iterate over the matrix differently. For a well-conditioned matrix, I doubt it makes any difference.

X -= X_solve;
}
// If we fail to converge after max_iter iterations, return -1
Expand Down