diff --git a/Makefile b/Makefile index b5475c4..1e66863 100644 --- a/Makefile +++ b/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 diff --git a/linalg.c b/linalg.c index 69b857f..ba0a30e 100644 --- a/linalg.c +++ b/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)); } }