add nondiscreteness tests

Makefile

@@ -1,17 +1,23 @@
 HEADERS=triangle.h linalg.h queue.h
 
-SPECIAL_OPTIONS=-O0 -g -D_DEBUG
+#SPECIAL_OPTIONS=-O0 -g -D_DEBUG
 #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
-#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
+SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 #SPECIAL_OPTIONS=
 
 OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
 
-all: discreteness singular_values element_path limit_set hyperbolic ellipticity
+all: discreteness singular_values nondiscrete traces element_path limit_set hyperbolic ellipticity
 
 singular_values: singular_values.o triangle.o linalg.o
 	gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas
 
+nondiscrete: nondiscrete.o triangle.o linalg.o
+	gcc $(OPTIONS) -o nondiscrete triangle.o linalg.o nondiscrete.o -lm -lgsl -lcblas
+
+traces: traces.o triangle.o linalg.o
+	gcc $(OPTIONS) -o traces triangle.o linalg.o traces.o -lm -lgsl -lcblas
+
 discreteness: discreteness.o triangle.o linalg.o
 	gcc $(OPTIONS) -o discreteness triangle.o linalg.o discreteness.o -lm -lgsl -lcblas
 
@@ -30,6 +36,12 @@ ellipticity: ellipticity.o triangle.o linalg.o
 singular_values.o: singular_values.c $(HEADERS)
 	gcc $(OPTIONS) -c singular_values.c
 
+nondiscrete.o: nondiscrete.c $(HEADERS)
+	gcc $(OPTIONS) -c nondiscrete.c
+
+traces.o: traces.c $(HEADERS)
+	gcc $(OPTIONS) -c traces.c
+
 discreteness.o: discreteness.c $(HEADERS)
 	gcc $(OPTIONS) -c discreteness.c
 
@@ -52,4 +64,4 @@ ellipticity.o: ellipticity.c $(HEADERS)
 	gcc $(OPTIONS) -c ellipticity.c
 
 clean:
-	rm -f singular_values discreteness element_path limit_set hyperbolic ellipticity triangle.o linalg.o singular_values.o discreteness.o element_path.o limit_set.o hyperbolic.o ellipticity.o
+	rm -f singular_values nondiscrete traces discreteness element_path limit_set hyperbolic ellipticity triangle.o linalg.o singular_values.o nondiscrete.o traces.o discreteness.o element_path.o limit_set.o hyperbolic.o ellipticity.o

linalg.c

@@ -145,24 +145,25 @@ double determinant(gsl_matrix *g, workspace_t *ws)
   return gsl_linalg_LU_det(ws->tmp, s);
 }
 
-void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
+void jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws)
 {
   gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv);
   gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
-  gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
+  gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_VAL_DESC);
 
   int real = 1;
-  for(int i = 0; i < 4; i++)
+  for(int i = 0; i < ws->n; i++)
     if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
       real = 0;
 
-  if(!real) {
+/*  if(!real) {
     printf("We have non-real eigenvalues!\n");
     exit(1);
-  }
+    }*/
 
-  for(int i = 0; i < ws->n - 1; i++) {
-    mu[i] = log(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i+1)));
+  for(int i = 0; i < ws->n; i++) {
+	  evs[2*i] = GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i));
+	  evs[2*i+1] = GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i));
   }
 }