From 78a3a7b9c24ef099297fc038b7acff281aef7ff8 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 2 Apr 2024 23:28:14 -0400
Subject: [PATCH] add optional argument to enable overriding sort order or
 eigenvalues/vectors

---
 src/math_eigen.cpp | 24 ++++++++++++++++++++----
 src/math_eigen.h   |  5 +++--
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/src/math_eigen.cpp b/src/math_eigen.cpp
index 65c3fa806a..c79a61f84a 100644
--- a/src/math_eigen.cpp
+++ b/src/math_eigen.cpp
@@ -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
 
diff --git a/src/math_eigen.h b/src/math_eigen.h
index 9f2431e14a..c7a3853767 100644
--- a/src/math_eigen.h
+++ b/src/math_eigen.h
@@ -22,13 +22,14 @@ namespace MathEigen {
  * \param  mat   the 3x3 matrix you wish to diagonalize
  * \param  eval  store the eigenvalues here
  * \param  evec  store the eigenvectors here...
+ * \param  sort  order eigenvalues and -vectors (-1 decreasing (default), 1 increasing, 0 unsorted)
  * \return       0 if eigenvalue calculation converged, 1 if it failed */
 
-int jacobi3(double const *const *mat, double *eval, double **evec);
+int jacobi3(double const *const *mat, double *eval, double **evec, int sort = -1);
 
 /** \overload */
 
-int jacobi3(double const mat[3][3], double *eval, double evec[3][3]);
+int jacobi3(double const mat[3][3], double *eval, double evec[3][3], int sort = -1);
 
 }    // namespace MathEigen