Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use Eigen on CPUs. #1349

Open
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
Copy link
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).

Copy link
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
Loading