some triangle limit set experiments

Makefile

@@ -7,19 +7,31 @@ SPECIAL_OPTIONS=-O0 -g -D_DEBUG
 
 OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
 
-all: triangle_group
+all: triangle_group element_path limit_set
 
 triangle_group: triangle.o linalg.o main.o
 	gcc $(OPTIONS) -o triangle_group triangle.o linalg.o main.o -lm -lgsl -lcblas
 
+element_path: element_path.o linalg.o
+	gcc $(OPTIONS) -o element_path element_path.o linalg.o -lm -lgsl -lcblas
+
+limit_set: limit_set.o linalg.o triangle.o
+	gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o -lm -lgsl -lcblas
+
 main.o: main.c $(HEADERS)
 	gcc $(OPTIONS) -c main.c
 
+element_path.o: element_path.c $(HEADERS)
+	gcc $(OPTIONS) -c element_path.c
+
 linalg.o: linalg.c $(HEADERS)
 	gcc $(OPTIONS) -c linalg.c
 
 triangle.o: triangle.c $(HEADERS)
 	gcc $(OPTIONS) -c triangle.c
 
+limit_set.o: limit_set.c $(HEADERS)
+	gcc $(OPTIONS) -c limit_set.c
+
 clean:
-	rm -f triangle_group triangle.o linalg.o main.o
+	rm -f triangle_group element_path limit_set triangle.o linalg.o main.o element_path.o limit_set.o

linalg.c

@@ -136,3 +136,27 @@ void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
     mu[i] = log(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i+1)));
   }
 }
+
+void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
+{
+	gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
+	int r = gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
+	ERROR(r, "gsl_eigen_nonsymmv failed!\n");
+
+	gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
+
+	int real = 1;
+	for(int i = 0; i < 3; i++)
+		if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
+			real = 0;
+
+	if(!real) {
+		printf("We have non-real eigenvalues!\n");
+		exit(1);
+	}
+
+	for(int i = 0; i < 3; i++)
+	  for(int j = 0; j < 3; j++)
+		  gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
+
+}

linalg.h

@@ -36,5 +36,6 @@ void rotation_matrix(gsl_matrix *g, double *vector);
 void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
 double trace(gsl_matrix *g);
 double determinant(gsl_matrix *g, workspace_t *ws);
+void eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
 
 #endif

main.c

@@ -1,7 +1,7 @@
 #include "triangle.h"
 #include "linalg.h"
 
-#define MAX_ELEMENTS 20000
+#define MAX_ELEMENTS 10000
 
 void initialize_triangle_generators(gsl_matrix **gen, double a1, double a2, double a3, double s)
 {
@@ -9,15 +9,15 @@ void initialize_triangle_generators(gsl_matrix **gen, double a1, double a2, doub
     gsl_matrix_set_identity(gen[i]);
 
   gsl_matrix_set(gen[0], 0, 0, -1);
-  gsl_matrix_set(gen[0], 1, 0, 2*s*cos(M_PI*a1));
-  gsl_matrix_set(gen[0], 2, 0, 2/s*cos(M_PI*a2));
+  gsl_matrix_set(gen[0], 0, 1, 2*s*cos(M_PI*a1));
+  gsl_matrix_set(gen[0], 0, 2, 2/s*cos(M_PI*a2));
 
-  gsl_matrix_set(gen[1], 0, 1, 2/s*cos(M_PI*a1));
+  gsl_matrix_set(gen[1], 1, 0, 2/s*cos(M_PI*a1));
   gsl_matrix_set(gen[1], 1, 1, -1);
-  gsl_matrix_set(gen[1], 2, 1, 2*s*cos(M_PI*a3));
+  gsl_matrix_set(gen[1], 1, 2, 2*s*cos(M_PI*a3));
 
-  gsl_matrix_set(gen[2], 0, 2, 2*s*cos(M_PI*a2));
-  gsl_matrix_set(gen[2], 1, 2, 2/s*cos(M_PI*a3));
+  gsl_matrix_set(gen[2], 2, 0, 2*s*cos(M_PI*a2));
+  gsl_matrix_set(gen[2], 2, 1, 2/s*cos(M_PI*a3));
   gsl_matrix_set(gen[2], 2, 2, -1);
 }
 
@@ -64,7 +64,7 @@ int main(int argc, char *argv[])
   // the real action
   generate_triangle_group(group, MAX_ELEMENTS, 5, 9, 7);
   //  initialize_triangle_generators(gen, 1.0/5, 1.0/5, 1.0/5, 1.0); // Fuchsian
-  initialize_triangle_generators(gen, 3.0/5.0, 3.0/5.0, 3.0/5.0, atof(argv[1]));
+  initialize_triangle_generators(gen, 2.0/5.0, 2.0/5.0, 2.0/5.0, atof(argv[1]));
 
   gsl_matrix_set_identity(matrices[0]);
   for(int i = 1; i < MAX_ELEMENTS; i++)

triangle_group.plt

@@ -1,5 +1,5 @@
 #if(!exists("param")) param = log(0.3810364354550512)
-if(!exists("param")) param = log(2.624420939708161)
+if(!exists("param")) param = log(2.764)
 
 print sprintf("param = %f", exp(param))
 file = sprintf("< ./triangle_group %f", exp(param))
@@ -11,10 +11,10 @@ set yrange [0:20]
 set grid
 set parametric
 
-plot file using 2:3 w p pt 7 ps 1 lc 1 t ""
+plot file using 2:3 w p pt 7 ps 0.5 lc 1 t ""
 
 pause mouse keypress
 if(MOUSE_KEY == 60) param=param-0.02
 if(MOUSE_KEY == 62) param=param+0.02
-if(MOUSE_KEY == 13) param=0
+if(MOUSE_KEY == 13) param=log(2.764)
 if(MOUSE_KEY != 113) reread