ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add optional argument to enable overriding sort order or eigenvalues/vectors

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-04-02 23:28:14 -04:00
parent 21745538a7
commit 78a3a7b9c2
2 changed files with 23 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/math_eigen.cpp
View File

@ -31,7 +31,7 @@ using namespace MathEigen;
typedef Jacobi<double, double *, double (*)[3], double const (*)[3]> Jacobi_v1;
typedef Jacobi<double, double *, double **, double const *const *> Jacobi_v2;
int MathEigen::jacobi3(double const mat[3][3], double *eval, double evec[3][3])
int MathEigen::jacobi3(double const mat[3][3], double *eval, double evec[3][3], int sort)
{
// make copy of const matrix
@ -44,7 +44,15 @@ int MathEigen::jacobi3(double const mat[3][3], double *eval, double evec[3][3])
// create instance of generic Jacobi class and get eigenvalues and -vectors
Jacobi_v1 ecalc3(3, M, midx);
int ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v1::SORT_DECREASING_EVALS);
int ierror = 1;
if (sort == -1)
ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v1::SORT_DECREASING_EVALS);
else if (sort == 0)
ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v1::DO_NOT_SORT);
else if (sort == 1)
ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v1::SORT_INCREASING_EVALS);
if (ierror) return ierror;
// transpose the evec matrix
@ -54,7 +62,7 @@ int MathEigen::jacobi3(double const mat[3][3], double *eval, double evec[3][3])
return ierror;
}
int MathEigen::jacobi3(double const *const *mat, double *eval, double **evec)
int MathEigen::jacobi3(double const *const *mat, double *eval, double **evec, int sort)
{
// make copy of const matrix
@ -67,7 +75,15 @@ int MathEigen::jacobi3(double const *const *mat, double *eval, double **evec)
// create instance of generic Jacobi class and get eigenvalues and -vectors
Jacobi_v2 ecalc3(3, M, midx);
int ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v2::SORT_DECREASING_EVALS);
int ierror = 1;
if (sort == -1)
ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v2::SORT_DECREASING_EVALS);
else if (sort == 0)
ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v2::DO_NOT_SORT);
else if (sort == 1)
ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v2::SORT_INCREASING_EVALS);
if (ierror) return ierror;
// transpose the evec matrix