Replace eigensolver code in LAMMPS with math_eigen.h and updated docs #2347

jewettaij · 2020-09-06T04:31:09Z

Summary

This pull-request replaces all Numerical Recipes eigensolver code in LAMMPS with the code in "math_eigen.h", which contains a variety of eigensolvers for matrices of varying size and sparsity. It also moves "math_eigen.h" into the main "src" directory, and adds documentation (in pg_utils.rst and doxygen).

(The code in lib/colvars directory was not modified. However I have made a separate pull request to replace the Numerical Recipes code located there as well. If that pull request is accepted, hopefully that code will eventually migrate its way into LAMMPS later on.)

Related Issues

fixes #2321

Author(s)

Andrew Jewett, Scripps Research, https://moltemplate.org

Yuya Kurebayashi, Tohoku University, https://github.com/mrcdr

Licensing

By submitting this pull request, I agree, that my contribution will be included in LAMMPS and redistributed under either the GNU General Public License version 2 (GPL v2) or the GNU Lesser General Public License version 2.1 (LGPL v2.1). The code in this pull request is also compatible with the MIT license.

Backward Compatibility

Effort was made to make minimal modifications to the existing LAMMPS code. No header files were changed (other than "math_extra.h"). The "jacobi3()" functions are stand-in replacements for the "jacobi()" functions found in "math_extra.h/cpp" and "fix_poems.cpp".

As for the "math_eigen.h" file, the eigensolvers there have been tested, and their unit-tests are available at their respective github pages:

Based on that code, Yuya and I have written a (somewhat messy) rough draft of the unit test code for "math_eigen.h". But it is not included in this pull-request.

Implementation Notes

The largest bug risk comes from possible differences in the order in which the eigenvalues and eigenvectors are returned. In the new "jacobi3()" function, the eigenvalues (and eigenvectors) are sorted according to eigenvalue in decreasing order. In the original jacobi() functions, it seems that no attempt was made to sort the eigenvalues. I don't think that this discrepancy would lead to behavioral problems, but I point it out. (EDIT: Sorting the eigenvalues does indeed effect the behavior rigid-body simulations, but does not seem to cause any harm. See below.)

This is a template class library. All implementations are fully separated from their declarations, however both currently appear in the header file ("math_eigen.h"). To avoid bloat in the compiled binary, the most common template class instances for float and double are defined explicitly in the "math_eigen.cpp" file. When you wish to use one of these instances, you can avoid recompiling the template class again by declaring the instance you need using the "extern" keyword:

// --- hypothetical .cpp file ----
#include"math_eigen.h"
extern template class MathEigen::Jacobi<double, double*, double**, double const*const*>;

(Tiresome detail: This is not necessary if you are invoking the jacobi3() function on 3x3 matrices. Since most of the current LAMMPS files use jacobi3(), I haven't used the "extern" declaration in any of the existing cpp files that I modified so far. But this option exists for people who want to diagonalize larger matrices and avoid bloat.)

More details: Why this approach? (feel free to skip)

If you move the implementation details into "math_eigen.cpp", it's a headache to get the PEigenDense and LambdaLanczos to work with complex valued matrices. (I was not able to get LAMMPS to compile if I simply add this line: "template class LambdaLanczos<std::complex<double>>;" to "math_eigen.cpp". SEE UPDATE BELOW.)
Generally speaking, the ability for these template classes to work with a variety of exotic matrix implementations (eg double**, double [3][3], cvm::matrix2d, as well as complex-valued matrices...) is one of their strengths. For this reason, I actually prefer to keep the implementations in the header file instead of moving them to a cpp file. Keeping implementations in the .h file and using the "extern" keyword seems to be a good compromise that avoids most of the bloat, but keeps the option open in the future for people to use new matrix representations. (For example, fix_poems.cpp uses double**. fix_rigid.cpp uses double [3][3]. The colvars library defines its own "cvm::matrix2d" type. The new code from "math_eigen.h" accepts all of these matrices directly without asking the user to convert them to double** first.)

Incidentally, the "math_eigen.h" file is only included in a few cpp files that need to calculate eigenvalues. So keeping its implementation in the header file should not slow down compilation noticeably. (EDIT: Compiling LAMMPS takes about 3 seconds longer now using this new code.)

UPDATE: Declaring a variable of type LambdaLanczos<std::complex<double>> inside a function does work. I suppose I could define an empty function within "math_eigen.cpp" which contains a local variable of this type. But it's a confusing workaround, and I'm happy keeping everything in "math_eigen.h".

If you still want to move all implementation details into "math_eigen.cpp", just cut lines 423-1368 from "math_eigen.h" and paste them into "math_eigen.cpp".

Post Submission Checklist

Further Information, Files, and Links

As mentioned earlier, the two main eigensolvers in this pull request are from these github repositorries:

https://github.com/jewettaij/jacobi_pd

https://github.com/mrcdr/lambda-lanczos

…ead of "math_extra.h". Code in "lib" now uses its own abbreviated local version of the "math_eigen.h" file (which is named "jacobi_pd.h"), since it is not clear how code from "lib/" can access the code in "src/"

…en.h". moved "math_eigen.h" into the main "src" directory.

… LAMMPS doxygen style.

…"false_positives.txt"

I replaced the 'Σ' character that was causing problems with PDF generation with the word "sum". Hopefully PDF documentation generation works now.

jewettaij · 2020-09-06T04:59:16Z

The unit tests succeed but the compilation tests fail. Interesting. For what it's worth, the code and the documentation (HTML and PDF) build successfully on my computer (running ubuntu 18.04, using conventional make).

akohlmey · 2020-09-06T08:59:42Z

The unit tests succeed but the compilation tests fail. Interesting. For what it's worth, the code and the documentation (HTML and PDF) build successfully on my computer (running ubuntu 18.04, using conventional make).

They fail because of a deficiency in superpose3d.h which uses assertions without including the corresponding header.

https://ci.lammps.org/job/dev/job/pull_requests/job/centos7/job/cmake_mpi_openmp_bigbig_static/1233/console

jewettaij · 2020-09-06T09:34:05Z

Thanks Axel. I did not think to test the code with USER-REACTION enabled. I'm too sleepy to wait for the test results, but hopefully this fixes everything.

akohlmey · 2020-09-06T09:34:42Z

Thanks Axel. I did not think to test the code with USER-REACTION enabled. I'm too sleepy to wait for the test results, but hopefully this fixes everything.

You also have merge conflicts and whitespace/permission issues

…" to make sphinx happy. Fixed some new spelling mistakes.

jewettaij · 2020-09-06T16:41:33Z

Thanks Axel. I did not think to test the code with USER-REACTION enabled. I'm too sleepy to wait for the test results, but hopefully this fixes everything.

You also have merge conflicts and whitespace/permission issues

...and some spelling errors I missed.
(...How did "peachpuff" end up in "false_positives.txt"?... no matter)

Also: When I run "make spelling" I also see that sphinx is complaining about line 830 of "pg_developer.rst" (which is in the "String processing" section). It seems that the real problem is that sphinx doesn't like the default string argument on line 239 of "utils.h". So I got around this problem this by editing "utils.h" and creating two different versions of the "trim_and_count_words()" function.

jewettaij · 2020-09-06T16:52:38Z

Oops. There are some linker errors. Hold on. I'll fix this in a few minutes. (Sorry. I got impatient waiting for LAMMPS to compile before pushing my commit because changing "utils.h" triggers a full recompilation.) Fixing this now.

akohlmey · 2020-09-06T17:13:47Z

Oops. There are some linker errors. Hold on. I'll fix this in a few minutes. (Sorry. I got impatient waiting for LAMMPS to compile before pushing my commit because changing "utils.h" triggers a full recompilation.) Fixing this now.

You are trying to fix things the wrong way around. You are 220 commits behind master and should merge in upstream first and with that resolve the merge conflicts. That would have likely removed the issue with the documentation build as well.

akohlmey · 2020-09-06T17:51:15Z

Oops. There are some linker errors. Hold on. I'll fix this in a few minutes. (Sorry. I got impatient waiting for LAMMPS to compile before pushing my commit because changing "utils.h" triggers a full recompilation.) Fixing this now.

This is not the fix we want. The current master can translate and link this just fine.

jewettaij · 2020-09-06T18:12:38Z

Oops. There are some linker errors. Hold on. I'll fix this in a few minutes. (Sorry. I got impatient waiting for LAMMPS to compile before pushing my commit because changing "utils.h" triggers a full recompilation.) Fixing this now.

This is not the fix we want. The current master can translate and link this just fine.

I will revert "utils.h" and "utils.cpp".

jewettaij · 2020-09-06T18:27:45Z

current with remote, finally.

jewettaij · 2020-09-06T18:40:09Z

Incidentally, after reverting to the previous version of "utils.h" and "utils.cpp", the message in big red letters returns. (But it looks like it's probably a warning message, not an error message, so no need to worry.)

doc/src/pg_developer.rst

jewettaij · 2020-09-15T01:24:26Z

Sounds good.
Let me know if you need anything.

akohlmey · 2020-09-16T16:09:56Z

@jewettaij I had a closer look at the code and the changes it introduces and thought about how to best address the concerns I have about all-header implementations versus the needs of LAMMPS. I think the best approach is to have the "math_eigen.h" header and the "math_eigen.cpp" file simply contain the prototypes for the two variants of the jacobi3 function with the optional arguments removed and applied directly. What is the "math_eigen.h" file will be renamed to "math_eigen_impl.h" and "math_eigen.cpp" (but not "math_eigen.h") will include it, same to "superpose3d.h". That should avoid major changes to the implementation and at the same time satisfy my desire to keep things clean and lean.

Thanks for all the effort with adding documentation, but I think for the LAMMPS purposes, it is sufficient to just document the jacobi3 function instead of the entire set of template classes, since they are currently not used anywhere else. I do want to keep mentioning them and the fact that they can be accessed in "math_eigen_impl.h", but for details I would prefer to point people to some webpage. Do you have a preference? The GitHub page where you maintain the code? or some other location?
Please let me know.

I'm also updating the unit tests and will take care of resetting the affected regression tests.

akohlmey · 2020-09-16T18:19:10Z

@sjplimp @rbberger this is now satisfying my preferences and seems to be a decent replacement for the previous version. Please let me know if there is anything else you would like to discuss or change. Once this is merged we need to reset some of the regression tests. In the long run, this should be for the better, because the eigenvectors are now consistently sorted by magnitude of the eigenvalues. The old implementation would provide them in different order on some platforms when the values were close and thus making cross-platform testing a bit more difficult.

sjplimp

The math_eigen.cpp file needs the standard LAMMPS license info
and contributing authors entry.

akohlmey · 2020-09-16T22:24:27Z

@sjplimp copyright and author attributes are added. I've also added URLs to the two projects on github where the code is maintained and integrated the existing tester tool (for the generic code, the specific functions are implicitly tested by the atom style tester and the fix rigid testers) into our unittest framework.

jewettaij · 2020-09-17T00:47:46Z

@sjplimp @rbberger this is now satisfying my preferences and seems to be a decent replacement for the previous version.

Cool.

Please let me know if there is anything else you would like to discuss or change. Once this is merged we need to reset some of the regression tests. In the long run, this should be for the better, because the eigenvectors are now consistently sorted by magnitude of the eigenvalues.

To clarify: (by default) they are sorted by the eigenvalues. If you prefer to sort them by their absolute value (by default), then edit lines 45 and 69 of "math_eigen.cpp" and 181, 210 of "math_eigen_impl.h" and replace: "SORT_DECREASING_EVALS" with "SORT_DECREASING_ABS_EVALS".
(I'm happy with either choice. EDIT: I just remembered the eigenvalues are positive when calculating moments of inertia, so it doesn't matter.)

The old implementation would provide them in different order on some platforms when the values were close and thus making cross-platform testing a bit more difficult.

jewettaij added 6 commits September 4, 2020 04:10

updated math_eigen.h and moved it into the main "src" directory

7f6fc8a

replaced all numerical recipes jacobi() code with code from "math_eig…

58fa72d

…en.h". moved "math_eigen.h" into the main "src" directory.

fixed a type error in math_eigen.h

4b66414

some corrections before rewriting the doxygen comments to comply with…

7cb774e

… LAMMPS doxygen style.

updated documentation files: "pg_developer.rst", "math_eigen.h", and …

4ad6d99

…"false_positives.txt"

jewettaij requested review from rbberger and sjplimp as code owners September 6, 2020 04:31

purged greek characters from "pd_developer.rst"

c16321e

I replaced the 'Σ' character that was causing problems with PDF generation with the word "sum". Hopefully PDF documentation generation works now.

removed assert() statements from "superpose3d.h"

2ee2009

changed the "trim_and_count_words()" function declaration in "utils.h…

b295222

…" to make sphinx happy. Fixed some new spelling mistakes.

jewettaij requested a review from akohlmey as a code owner September 6, 2020 16:41

fixed linker problems in "utils.h" and "utils.cpp"

d97e8d8

jewettaij added 2 commits September 6, 2020 10:59

Merge branch 'master' into math_eigen

eb6f7f9

reverting to previous version of util.h and util.cpp

fabf762

akohlmey reviewed Sep 6, 2020

View reviewed changes

doc/src/pg_developer.rst Outdated Show resolved Hide resolved

akohlmey added enhancement test_for_regression Enable to trigger regression tests maintenance labels Sep 6, 2020

akohlmey added 4 commits September 15, 2020 21:37

Merge branch 'master' into math_eigen

f599869

only do Zstd tests if it is enabled/found.

48ad860

update rigid tests for new eigensolver

981d60f

update atom_style tester for new jacobi implementation

f7a939d

akohlmey self-assigned this Sep 16, 2020

akohlmey added 4 commits September 16, 2020 14:06

separate the MathEigen implementation into its own header file

10991ee

update docs

fd9e39b

relax epsilon for rigid body integrators

0931cdd

include math_eigen.h, remove it from .gitignore

bd542ef

akohlmey removed the work_in_progress label Sep 16, 2020

akohlmey added this to the Stable Release Summer 2020 milestone Sep 16, 2020

akohlmey self-requested a review September 16, 2020 18:15

sjplimp requested changes Sep 16, 2020

View reviewed changes

akohlmey added 4 commits September 16, 2020 17:35

update copyright headers to match LAMMPS' conventions

942ed3a

add URLs to upstream projects for eigensolvers

96f4178

adapt unit test for Jacobi class

a8a9fb6

remove trailing whitespace

b176cdf

akohlmey changed the title ~~replaced all eigensolver code in LAMMPS with math_eigen.h and updated docs~~ Replace eigensolver code in LAMMPS with math_eigen.h and updated docs Sep 16, 2020

akohlmey requested a review from sjplimp September 16, 2020 22:19

akohlmey added the ready_for_merge label Sep 16, 2020

sjplimp approved these changes Sep 16, 2020

View reviewed changes

akohlmey approved these changes Sep 17, 2020

View reviewed changes

akohlmey merged commit e924fc6 into lammps:master Sep 17, 2020

giacomofiorin mentioned this pull request Sep 18, 2020

Update Colvars to version 2020-09-17 #2376

Merged

7 tasks

mrcdr mentioned this pull request Oct 5, 2020

PR to relate "math_eigen_impl.h" with an author #2409

Merged

11 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Replace eigensolver code in LAMMPS with math_eigen.h and updated docs #2347

Replace eigensolver code in LAMMPS with math_eigen.h and updated docs #2347

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 6, 2020

akohlmey commented Sep 6, 2020

jewettaij commented Sep 6, 2020

akohlmey commented Sep 6, 2020

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 6, 2020

akohlmey commented Sep 6, 2020

akohlmey commented Sep 6, 2020

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 6, 2020

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 15, 2020

akohlmey commented Sep 16, 2020

akohlmey commented Sep 16, 2020

sjplimp left a comment

akohlmey commented Sep 16, 2020

jewettaij commented Sep 17, 2020 •

edited

Loading

Replace eigensolver code in LAMMPS with math_eigen.h and updated docs #2347

Replace eigensolver code in LAMMPS with math_eigen.h and updated docs #2347

Conversation

jewettaij commented Sep 6, 2020 • edited Loading

jewettaij commented Sep 6, 2020

akohlmey commented Sep 6, 2020

jewettaij commented Sep 6, 2020

akohlmey commented Sep 6, 2020

jewettaij commented Sep 6, 2020 • edited Loading

jewettaij commented Sep 6, 2020

akohlmey commented Sep 6, 2020

akohlmey commented Sep 6, 2020

jewettaij commented Sep 6, 2020 • edited Loading

jewettaij commented Sep 6, 2020

jewettaij commented Sep 6, 2020 • edited Loading

jewettaij commented Sep 15, 2020

akohlmey commented Sep 16, 2020

akohlmey commented Sep 16, 2020

sjplimp left a comment

Choose a reason for hiding this comment

akohlmey commented Sep 16, 2020

jewettaij commented Sep 17, 2020 • edited Loading

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 6, 2020 •

edited

Loading

jewettaij commented Sep 17, 2020 •

edited

Loading