From 79347f3fb612e203e68e19a708ae91f1d679c2d9 Mon Sep 17 00:00:00 2001
From: Victor Zverovich
Date: Mon, 18 Aug 2014 14:23:16 -0700
Subject: [PATCH] Sort eigenvalues by real part first and then by imaginary
part to make eigen_test more robust.
---
eigen/sort.c | 13 +++++++++++++
eigen/test.c | 4 ++--
2 files changed, 15 insertions(+), 2 deletions(-)
diff --git a/eigen/sort.c b/eigen/sort.c
index 54d1316..0eb5b40 100644
--- a/eigen/sort.c
+++ b/eigen/sort.c
@@ -26,6 +26,15 @@
#include
#include
+/* Compares real parts of a and b and, if they are not approximately equal,
+ * returns Re(a) < Re(b); otherwise returns Im(a) < Im(b). */
+static INLINE_DECL int
+complex_less(gsl_complex a, gsl_complex b)
+{
+ return gsl_fcmp(GSL_REAL(a), GSL_REAL(b), GSL_DBL_EPSILON) == 0 ?
+ GSL_IMAG(a) < GSL_IMAG(b) : GSL_REAL(a) < GSL_REAL(b);
+}
+
/* The eigen_sort below is not very good, but it is simple and
* self-contained. We can always implement an improved sort later. */
@@ -208,7 +217,11 @@ gsl_eigen_nonsymmv_sort (gsl_vector_complex * eval,
test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
break;
case GSL_EIGEN_SORT_VAL_ASC:
+ test = complex_less(ej, ek);
+ break;
case GSL_EIGEN_SORT_VAL_DESC:
+ test = complex_less(ek, ej);
+ break;
default:
GSL_ERROR ("invalid sort type", GSL_EINVAL);
}
diff --git a/eigen/test.c b/eigen/test.c
index 20b9f8a..0cc291a 100644
--- a/eigen/test.c
+++ b/eigen/test.c
@@ -1267,8 +1267,8 @@ test_eigen_gen_pencil(const gsl_matrix * A, const gsl_matrix * B,
}
/* sort eval and evalv and test them */
- gsl_eigen_nonsymmv_sort(w->eval, NULL, GSL_EIGEN_SORT_ABS_ASC);
- gsl_eigen_nonsymmv_sort(w->evalv, NULL, GSL_EIGEN_SORT_ABS_ASC);
+ gsl_eigen_nonsymmv_sort(w->eval, NULL, GSL_EIGEN_SORT_VAL_ASC);
+ gsl_eigen_nonsymmv_sort(w->evalv, NULL, GSL_EIGEN_SORT_VAL_ASC);
test_eigenvalues_complex(w->evalv, w->eval, "gen", desc);
gsl_eigen_genv_sort(w->alphav, w->betav, w->evec, GSL_EIGEN_SORT_ABS_ASC);
--
1.9.1