From 5eb23c1985c212f517cd2409df2a0df963a5e6de Mon Sep 17 00:00:00 2001
From: Cory <cory@balaton.dev>
Date: Wed, 20 Sep 2023 20:00:07 +0200
Subject: [PATCH] Sort eigenvalues and eigenvectors

---
 src/jacobi.cpp | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/jacobi.cpp b/src/jacobi.cpp
index 00fd9eb..0542c18 100644
--- a/src/jacobi.cpp
+++ b/src/jacobi.cpp
@@ -6,7 +6,9 @@
  *  @bug The eigenvalues fail the test. When comparing them to arma::eigsym
  *  there is one value that is way off when testing with a 6x6 matrix.
  */
+#include <algorithm>
 #include <cmath>
+#include <utility>
 
 #include "utils.hpp"
 #include "jacobi.hpp"
@@ -75,9 +77,16 @@ void jacobi_eigensolver(const arma::mat& A,
 
     } while (max_offdiag >= eps && ++iterations < maxiter);
 
+
+    arma::vec res = arma::diagvec(A_m);
+    arma::uvec sorted_indeces = arma::sort_index(res);
+    eigenvalues.resize(res.n_elem);
+    eigenvectors.resize(res.n_elem, res.n_elem);
     // The diagonal elements of A_m are the eigenvalues
-    eigenvalues = arma::diagvec(A_m);
-    eigenvectors = R;
+    for (int i=0; i < res.n_elem; i++) {
+        eigenvalues[i] = res[sorted_indeces[i]]; 
+        eigenvectors.insert_cols(i, R.col(sorted_indeces[i]));
+    }
 
     converged = max_offdiag < eps;
     A_m.print();