lammps · akohlmey · Sep 17, 2020 · Sep 4, 2020 · Sep 5, 2020 · Sep 5, 2020
diff --git a/doc/doxygen/Doxyfile.in b/doc/doxygen/Doxyfile.in
@@ -431,6 +431,7 @@ INPUT                  = @LAMMPS_SOURCE_DIR@/utils.cpp                 \
                          @LAMMPS_SOURCE_DIR@/my_page.h                 \
                          @LAMMPS_SOURCE_DIR@/my_pool_chunk.cpp         \
                          @LAMMPS_SOURCE_DIR@/my_pool_chunk.h           \
+                         @LAMMPS_SOURCE_DIR@/math_eigen.h              \
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded

diff --git a/doc/src/pg_dev_utils.rst b/doc/src/pg_dev_utils.rst
@@ -10,7 +10,7 @@ strings into specific types of numbers with checking for validity.  This
 reduces redundant implementations and encourages consistent behavior.
 
 I/O with status check
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+^^^^^^^^^^^^^^^^^^^^^
 
 These are wrappers around the corresponding C library calls like
 ``fgets()`` or ``fread()``.  They will check if there were errors
@@ -26,6 +26,8 @@ indicating the name of the problematic file, if possible.
 .. doxygenfunction:: sfread
    :project: progguide
 
+----------
+
 String to number conversions with validity check
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -281,6 +283,8 @@ This code example should produce the following output:
    :project: progguide
    :members: what
 
+----------
+
 File reader classes
 ====================
 
@@ -333,6 +337,8 @@ convert numbers, so that LAMMPS will be aborted.
 
 A file that would be parsed by the reader code fragment looks like this:
 
+.. parsed-literal::
+
    # DATE: 2015-02-19 UNITS: metal CONTRIBUTOR: Ray Shan CITATION: Streitz and Mintmire, Phys Rev B, 50, 11996-12003 (1994)
    #
    # X (eV)                J (eV)          gamma (1/\AA)   zeta (1/\AA)    Z (e)
@@ -351,7 +357,6 @@ A file that would be parsed by the reader code fragment looks like this:
    :project: progguide
    :members:
 
-
 ----------
 
 Memory pool classes
@@ -415,3 +420,43 @@ its size is registered later with :cpp:func:`vgot()
 .. doxygenclass:: LAMMPS_NS::MyPoolChunk
    :project: progguide
    :members:
+
+----------
+
+Eigensolver functions
+=====================
+
+The ``MathEigen`` sub-namespace of the ``LAMMPS_NS`` namespace contains
+functions and classes for eigensolvers. Currently only the
+:cpp:func:`jacobi3 function <MathEigen::jacobi3>` is used in various
+places in LAMMPS.  That function is built on top of a group of more
+generic eigensolvers that are maintained in the ``math_eigen_impl.h``
+header file.  This header contains the implementation of three template
+classes:
+
+#. "Jacobi" calculates all of the eigenvalues and eigenvectors
+   of a dense, symmetric, real matrix.
+
+#. The "PEigenDense" class only calculates the principal eigenvalue
+   (ie. the largest or smallest eigenvalue), and its corresponding
+   eigenvector.  However it is much more efficient than "Jacobi" when
+   applied to large matrices (larger than 13x13).  PEigenDense also can
+   understand complex-valued Hermitian matrices.
+
+#. The "LambdaLanczos" class is a generalization of "PEigenDense" which can be
+   applied to arbitrary sparse matrices.
+
+The "math_eigen_impl.h" code is an amalgamation of `jacobi_pd
+<https://github.com/jewettaij/jacobi_pd>`_ by Andrew Jewett at Scripps
+Research (under CC0-1.0 license) and `Lambda Lanczos
+<https://github.com/mrcdr/lambda-lanczos>`_ by Yuya Kurebayashi at
+Tohoku University (under MIT license)
+
+----------
+
+.. doxygenfunction:: MathEigen::jacobi3(double const *const *mat, double *eval, double **evec)
+   :project: progguide
+
+.. doxygenfunction:: MathEigen::jacobi3(double const mat[3][3], double *eval, double evec[3][3])
+   :project: progguide
+
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
@@ -325,6 +325,7 @@ Buyl
 Bybee
 bz
 cadetblue
+calc
 calibre
 caltech
 Caltech
@@ -397,6 +398,7 @@ ChiralIDs
 chiralIDs
 chirality
 Cho
+ChooseOffset
 chris
 Christoph
 Chu
@@ -469,6 +471,7 @@ config
 configfile
 configurational
 conformational
+ConstMatrix
 Contrib
 cooperativity
 coord
@@ -639,7 +642,10 @@ dhex
 dia
 diag
 diagonalization
+diagonalize
 diagonalized
+diagonalizers
+diagonalizing
 Diallo
 diel
 differentiable
@@ -778,8 +784,15 @@ Eggebrecht
 ehex
 eHEX
 Ei
-Eigen
-Eigensolve
+eigen
+eigensolve
+eigensolver
+eigensolvers
+eigendecomposition
+eigenvalue
+eigenvalues
+eigenvector
+eigenvectors
 eij
 Eij
 Eijnden
@@ -893,6 +906,11 @@ eV
 evalue
 Evanseck
 evdwl
+evector
+evec
+evecs
+eval
+evals
 Everaers
 Evgeny
 evirials
@@ -1069,6 +1087,7 @@ Germann
 Germano
 gerolf
 Gerolf
+Gershgorin
 gettimeofday
 gewald
 Gezelter
@@ -1184,6 +1203,7 @@ Henkelman
 Henkes
 henrich
 Henrich
+Hermitian
 Herrmann
 Hertizian
 hertzian
@@ -1257,6 +1277,7 @@ icosahedral
 idealgas
 IDR
 idx
+ie
 ielement
 ieni
 ifdefs
@@ -1279,6 +1300,7 @@ Imageint
 Imagemagick
 imd
 Impey
+impl
 impropers
 Impropers
 includelink
@@ -1348,6 +1370,8 @@ isothermal
 isotropically
 isovolume
 Isralewitz
+iter
+iters
 iteratively
 Ith
 Itsets
@@ -1524,6 +1548,7 @@ Kub
 Kubo
 Kumagai
 Kumar
+Kurebayashi
 Kuronen
 Kusters
 Kutta
@@ -1534,12 +1559,14 @@ Ladd
 lagrangian
 lambdai
 lamda
+LambdaLanczos
 lammps
 Lammps
 LAMMPS
 lammpsplot
 Lampis
 Lamoureux
+Lanczos
 Lande
 Landron
 langevin
@@ -1847,6 +1874,7 @@ Microscale
 midnightblue
 mie
 Mie
+Mij
 Mikami
 Militzer
 Minary
@@ -1945,6 +1973,7 @@ Muccioli
 mui
 Mukherjee
 Mulders
+mult
 multi
 multibody
 Multibody
@@ -2331,6 +2360,7 @@ peachpuff
 Pearlman
 Pedersen
 peID
+PEigenDense
 Peng
 peptide
 peratom
@@ -2559,6 +2589,8 @@ rdf
 RDideal
 rdx
 reacter
+realTypeMap
+real_t
 README
 realtime
 reamin
@@ -2741,6 +2773,7 @@ Schuring
 Schwen
 screenshot
 screenshots
+Scripps
 Scripta
 sdk
 sdpd
@@ -3069,6 +3102,7 @@ Tmin
 tmp
 tN
 Tobias
+Tohoku
 tokenizer
 tokyo
 tol
@@ -3132,6 +3166,8 @@ tu
 Tuckerman
 tue
 tunable
+tuple
+tuples
 Turkand
 Tutein
 tweakable
@@ -3425,6 +3461,7 @@ Yuh
 yukawa
 Yukawa
 Yusof
+Yuya
 yx
 yy
 yz

diff --git a/src/.gitignore b/src/.gitignore
@@ -33,7 +33,6 @@
 /pair_kim.h
 
 /superpose3d.h
-/math_eigen.h
 
 /kokkos.cpp
 /kokkos.h

diff --git a/src/BODY/body_nparticle.cpp b/src/BODY/body_nparticle.cpp
@@ -17,6 +17,7 @@
 #include "atom_vec_body.h"
 #include "error.h"
 #include "math_extra.h"
+#include "math_eigen.h"
 #include "memory.h"
 #include "my_pool_chunk.h"
 
@@ -136,7 +137,7 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
 
   double *inertia = bonus->inertia;
   double evectors[3][3];
-  int ierror = MathExtra::jacobi(tensor,inertia,evectors);
+  int ierror = MathEigen::jacobi3(tensor,inertia,evectors);
   if (ierror) error->one(FLERR,
                          "Insufficient Jacobi rotations for body nparticle");
 

diff --git a/src/BODY/body_rounded_polygon.cpp b/src/BODY/body_rounded_polygon.cpp
@@ -22,6 +22,7 @@
 #include "domain.h"
 #include "error.h"
 #include "math_extra.h"
+#include "math_eigen.h"
 #include "memory.h"
 #include "my_pool_chunk.h"
 
@@ -198,7 +199,7 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
 
   double *inertia = bonus->inertia;
   double evectors[3][3];
-  int ierror = MathExtra::jacobi(tensor,inertia,evectors);
+  int ierror = MathEigen::jacobi3(tensor,inertia,evectors);
   if (ierror) error->one(FLERR,
                          "Insufficient Jacobi rotations for body nparticle");
 

diff --git a/src/BODY/body_rounded_polyhedron.cpp b/src/BODY/body_rounded_polyhedron.cpp
@@ -21,6 +21,7 @@
 #include "atom_vec_body.h"
 #include "error.h"
 #include "math_extra.h"
+#include "math_eigen.h"
 #include "memory.h"
 #include "my_pool_chunk.h"
 
@@ -242,7 +243,7 @@ void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
 
   double *inertia = bonus->inertia;
   double evectors[3][3];
-  int ierror = MathExtra::jacobi(tensor,inertia,evectors);
+  int ierror = MathEigen::jacobi3(tensor,inertia,evectors);
   if (ierror) error->one(FLERR,
                          "Insufficient Jacobi rotations for body nparticle");