these boxes could work!

Makefile

@@ -1,8 +1,8 @@
 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)

linalg.c

@@ -18,6 +18,7 @@ workspace_t *workspace_alloc(int n)
   workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t));
   result->n = n;
   result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n);
+  result->work_symmv = gsl_eigen_symmv_alloc(n);
   result->work_sv = gsl_vector_alloc(n);
   result->eval_complex = gsl_vector_complex_alloc(n);
   result->evec_complex = gsl_matrix_complex_alloc(n, n);
@@ -34,6 +35,7 @@ workspace_t *workspace_alloc(int n)
 void workspace_free(workspace_t *workspace)
 {
   gsl_eigen_nonsymmv_free(workspace->work_nonsymmv);
+  gsl_eigen_symmv_free(workspace->work_symmv);
   gsl_vector_free(workspace->work_sv);
   gsl_vector_complex_free(workspace->eval_complex);
   gsl_matrix_complex_free(workspace->evec_complex);
@@ -67,6 +69,33 @@ void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out)
   gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, out);
 }
 
+void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
+{
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, ws->stack[0]);
+	gsl_matrix_memcpy(a, ws->stack[0]);
+}
+
+void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
+{
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, ws->stack[0]);
+	gsl_matrix_memcpy(b, ws->stack[0]);
+}
+
+void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...)
+{
+	va_list args;
+    va_start(args, n);
+
+	gsl_matrix_set_identity(out);
+
+	for(int i = 0; i < n; i++) {
+		gsl_matrix *cur = va_arg(args, gsl_matrix *);
+		multiply_right(out, cur, ws);
+	}
+
+    va_end(args);
+}
+
 void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
 {
   gsl_matrix_memcpy(ws->tmp, g);
@@ -139,14 +168,15 @@ void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
 
 void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
 {
+	gsl_matrix_memcpy(ws->stack[0], g);
 	gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
-	int r = gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
+	int r = gsl_eigen_nonsymmv(ws->stack[0], 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++)
+	for(int i = 0; i < ws->n; i++)
 		if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
 			real = 0;
 
@@ -155,8 +185,45 @@ void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
 		exit(1);
 	}
 
-	for(int i = 0; i < 3; i++)
-	  for(int j = 0; j < 3; j++)
+	for(int i = 0; i < ws->n; i++)
+	  for(int j = 0; j < ws->n; j++)
 		  gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
 
 }
+
+void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)
+{
+	gsl_matrix_memcpy(ws->stack[0], g);
+	int r = gsl_eigen_symmv (ws->stack[0], eval, evec, ws->work_symmv);
+	ERROR(r, "gsl_eigen_symmv failed!\n");
+
+	gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
+}
+
+// returns number of positive directions and matrix transforming TO diagonal basis
+int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
+{
+	gsl_matrix_memcpy(ws->stack[0], A);
+	int r = gsl_eigen_symmv (ws->stack[0], ws->eval_real, cob, ws->work_symmv);
+	ERROR(r, "gsl_eigen_symmv failed!\n");
+
+	gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC);
+
+	gsl_matrix_transpose(cob);
+
+	int positive = 0;
+
+	for(int i = 0; i < ws->n; i++) {
+		if(gsl_vector_get(ws->eval_real, i) > 0)
+			positive++;
+
+		for(int j = 0; j < ws->n; j++)
+			*gsl_matrix_ptr(cob, i, j) *= sqrt(fabs(gsl_vector_get(ws->eval_real, i)));
+	}
+
+	return positive;
+
+//	printf("Eigenvalues: %.10f, %.10f, %.10f\n", gsl_vector_get(ws->eval_real, 0), gsl_vector_get(ws->eval_real, 1), gsl_vector_get(ws->eval_real, 2));
+
+//	return 0;
+}

linalg.h

@@ -3,6 +3,7 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <stdarg.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_blas.h>
@@ -15,6 +16,7 @@
 typedef struct _workspace {
 	int n;
 	gsl_eigen_nonsymmv_workspace *work_nonsymmv;
+	gsl_eigen_symmv_workspace *work_symmv;
 	gsl_vector *work_sv;
 	gsl_vector_complex *eval_complex;
 	gsl_matrix_complex *evec_complex;
@@ -30,6 +32,9 @@ void workspace_free(workspace_t *workspace);
 void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws);
 void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws);
 void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out);
+void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
+void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
+void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...);
 void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
 void initialize(gsl_matrix *g, double *data, int x, int y);
 void rotation_matrix(gsl_matrix *g, double *vector);
@@ -37,5 +42,7 @@ 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);
+void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws);
+int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws);
 
 #endif

triangle.c

@@ -28,6 +28,7 @@ int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int
 	int k[3] = {k1, k2, k3};
 
 	memset(group, 0, size*sizeof(groupelement_t));
+	group[0].letter = -1;
 	n = 1;
 	queue_init(&q);
 	queue_put(&q, 0);