From aea8496d6364d05b00801241e388eb0dc863735d Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Thu, 31 May 2018 00:13:53 +0200 Subject: [PATCH] these boxes could work! --- Makefile | 4 +-- linalg.c | 75 +++++++++++++++++++++++++++++++++++++++++++++++++++--- linalg.h | 27 ++++++++++++-------- triangle.c | 1 + 4 files changed, 91 insertions(+), 16 deletions(-) diff --git a/Makefile b/Makefile index ecd48d8..f651f9b 100644 --- a/Makefile +++ b/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) diff --git a/linalg.c b/linalg.c index 66ef82d..69b857f 100644 --- a/linalg.c +++ b/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; +} diff --git a/linalg.h b/linalg.h index b6cdb1f..11cbe80 100644 --- a/linalg.h +++ b/linalg.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -13,16 +14,17 @@ #define FCMP(x, y) gsl_fcmp(x, y, 1e-10) typedef struct _workspace { - int n; - gsl_eigen_nonsymmv_workspace *work_nonsymmv; - gsl_vector *work_sv; - gsl_vector_complex *eval_complex; - gsl_matrix_complex *evec_complex; - gsl_vector *eval_real; - gsl_matrix *evec_real; - gsl_matrix *tmp; - gsl_permutation *permutation; - gsl_matrix *stack[20]; + 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; + gsl_vector *eval_real; + gsl_matrix *evec_real; + gsl_matrix *tmp; + gsl_permutation *permutation; + gsl_matrix *stack[20]; } workspace_t; workspace_t *workspace_alloc(int n); @@ -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 diff --git a/triangle.c b/triangle.c index 5b23d10..770c7fe 100644 --- a/triangle.c +++ b/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);