diff --git a/draw.c b/draw.c index 248858e..c5eff01 100644 --- a/draw.c +++ b/draw.c @@ -13,22 +13,22 @@ vector_t cross(vector_t a, vector_t b) vector_t apply(DrawingContext *ctx, gsl_matrix *m, vector_t x) { - gsl_vector *tmp = getTempVector(ctx); - gsl_vector *tmp2 = getTempVector(ctx); + gsl_vector *tmp = getTempVector(ctx->ws); + gsl_vector *tmp2 = getTempVector(ctx->ws); vector_t out; LOOP(i) gsl_vector_set(tmp, i, x.x[i]); gsl_blas_dgemv(CblasNoTrans, 1.0, m, tmp, 0.0, tmp2); LOOP(i) out.x[i] = gsl_vector_get(tmp2, i); - releaseTempVectors(ctx, 2); + releaseTempVectors(ctx->ws, 2); } int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) { - gsl_matrix *tmp = getTempMatrix(ctx); - gsl_matrix *ev = getTempMatrix(ctx); - gsl_matrix **gen = getTempMatrices(ctx, 3); + gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_matrix *ev = getTempMatrix(ctx->ws); + gsl_matrix **gen = getTempMatrices(ctx->ws, 3); initializeTriangleGenerators(gen, ctx->cartan); @@ -42,16 +42,16 @@ int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) LOOP(i) LOOP(j) out[i].x[j] = gsl_matrix_get(ev, j, i); - releaseTempMatrices(ctx, 5); + releaseTempMatrices(ctx->ws, 5); return count; } void transformFrameStd(DrawingContext *ctx, vector_t *x, gsl_matrix *out) { - gsl_matrix *tmp = getTempMatrix(ctx); - gsl_vector *fourth = getTempVector(ctx); - gsl_vector *lambda = getTempVector(ctx); + gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_vector *fourth = getTempVector(ctx->ws); + gsl_vector *lambda = getTempVector(ctx->ws); int s; LOOP(i) LOOP(j) gsl_matrix_set(out, j, i, x[i].x[j]); @@ -63,8 +63,8 @@ void transformFrameStd(DrawingContext *ctx, vector_t *x, gsl_matrix *out) gsl_matrix_fprintf(stdout, out, "%f"); - releaseTempMatrices(ctx, 1); - releaseTempVectors(ctx, 2); + releaseTempMatrices(ctx->ws, 1); + releaseTempVectors(ctx->ws, 2); } // level 1: the elementary drawing functions, drawPoint, drawSegment2d diff --git a/limit_set.c b/limit_set.c index 03770f5..270098c 100644 --- a/limit_set.c +++ b/limit_set.c @@ -29,16 +29,16 @@ void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan) int computeLimitCurve(DrawingContext *ctx) { workspace_t *ws = ctx->ws; - gsl_matrix *cartan_pos = getTempMatrix(ctx); - gsl_matrix *cob_pos = getTempMatrix(ctx); - gsl_matrix *coxeter_pos = getTempMatrix(ctx); - gsl_matrix *coxeter_fixedpoints_pos = getTempMatrix(ctx); - gsl_matrix *fixedpoints_pos = getTempMatrix(ctx); - gsl_matrix *coxeter = getTempMatrix(ctx); - gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx); - gsl_matrix *fixedpoints = getTempMatrix(ctx); - gsl_matrix **gen = getTempMatrices(ctx, 3); - gsl_matrix **elements = getTempMatrices(ctx, ctx->n_group_elements); + gsl_matrix *cartan_pos = getTempMatrix(ctx->ws); + gsl_matrix *cob_pos = getTempMatrix(ctx->ws); + gsl_matrix *coxeter_pos = getTempMatrix(ctx->ws); + gsl_matrix *coxeter_fixedpoints_pos = getTempMatrix(ctx->ws); + gsl_matrix *fixedpoints_pos = getTempMatrix(ctx->ws); + gsl_matrix *coxeter = getTempMatrix(ctx->ws); + gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws); + gsl_matrix *fixedpoints = getTempMatrix(ctx->ws); + gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements); groupelement_t *group = ctx->group; int success = 0; @@ -92,7 +92,7 @@ int computeLimitCurve(DrawingContext *ctx) success = 1; error_out: - releaseTempMatrices(ctx, 11+ctx->n_group_elements); + releaseTempMatrices(ctx->ws, 11+ctx->n_group_elements); return success; } diff --git a/linalg.c b/linalg.c index 998496f..97a7be8 100644 --- a/linalg.c +++ b/linalg.c @@ -24,11 +24,16 @@ workspace_t *workspace_alloc(int n) result->evec_complex = gsl_matrix_complex_alloc(n, n); result->eval_real = gsl_vector_alloc(n); result->evec_real = gsl_matrix_alloc(n, n); - result->tmp = gsl_matrix_alloc(n, n); result->permutation = gsl_permutation_alloc(n); - for(int i = 0; i < 20; i++) - result->stack[i] = gsl_matrix_alloc(n, n); + result->tmp_mat = malloc(MAX_TEMP_MATRICES*sizeof(gsl_matrix*)); + for(int i = 0; i < MAX_TEMP_MATRICES; i++) + result->tmp_mat[i] = gsl_matrix_alloc(3, 3); + result->tmp_mat_used = 0; + result->tmp_vec = malloc(MAX_TEMP_MATRICES*sizeof(gsl_vector*)); + for(int i = 0; i < MAX_TEMP_MATRICES; i++) + result->tmp_vec[i] = gsl_vector_alloc(3); + result->tmp_vec_used = 0; return result; } @@ -41,10 +46,14 @@ void workspace_free(workspace_t *workspace) gsl_matrix_complex_free(workspace->evec_complex); gsl_vector_free(workspace->eval_real); gsl_matrix_free(workspace->evec_real); - gsl_matrix_free(workspace->tmp); gsl_permutation_free(workspace->permutation); - for(int i = 0; i < 20; i++) - gsl_matrix_free(workspace->stack[i]); + + for(int i = 0; i < MAX_TEMP_MATRICES; i++) + gsl_matrix_free(workspace->tmp_mat[i]); + free(workspace->tmp_mat); + for(int i = 0; i < MAX_TEMP_VECTORS; i++) + gsl_vector_free(workspace->tmp_vec[i]); + free(workspace->tmp_vec); } /************************************************** basic operations ********************************************************/ @@ -52,16 +61,24 @@ void workspace_free(workspace_t *workspace) void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws) { int s; - gsl_matrix_memcpy(ws->tmp, in); - gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s); - gsl_linalg_LU_invert(ws->tmp, ws->permutation, out); + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, in); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + gsl_linalg_LU_invert(tmp, ws->permutation, out); + + releaseTempMatrices(ws, 1); } void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws) { + gsl_matrix *tmp = getTempMatrix(ws); + invert(conjugator, out, ws); // use out to temporarily store inverse conjugator - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, ws->tmp); // in * conjugator^{-1} - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, ws->tmp, 0.0, out); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, tmp); // in * conjugator^{-1} + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, tmp, 0.0, out); + + releaseTempMatrices(ws, 1); } void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out) @@ -71,14 +88,18 @@ void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *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]); + gsl_matrix *tmp = getTempMatrix(ws); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, tmp); + gsl_matrix_memcpy(a, tmp); + releaseTempMatrices(ws, 1); } 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]); + gsl_matrix *tmp = getTempMatrix(ws); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, tmp); + gsl_matrix_memcpy(b, tmp); + releaseTempMatrices(ws, 1); } void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...) @@ -98,11 +119,15 @@ void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...) void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws) { - gsl_matrix_memcpy(ws->tmp, g); - gsl_linalg_SV_decomp(ws->tmp, ws->evec_real, ws->eval_real, ws->work_sv); + gsl_matrix *tmp = getTempMatrix(ws); - for(int i = 0; i < ws->n - 1; i++) - mu[i] = log(gsl_vector_get(ws->eval_real, i) / gsl_vector_get(ws->eval_real, i+1)); + gsl_matrix_memcpy(tmp, g); + gsl_linalg_SV_decomp(tmp, ws->evec_real, ws->eval_real, ws->work_sv); + + for(int i = 0; i < ws->n - 1; i++) + mu[i] = log(gsl_vector_get(ws->eval_real, i) / gsl_vector_get(ws->eval_real, i+1)); + + releaseTempMatrices(ws, 1); } void initialize(gsl_matrix *g, double *data, int x, int y) @@ -140,9 +165,15 @@ double trace(gsl_matrix *g) double determinant(gsl_matrix *g, workspace_t *ws) { int s; - gsl_matrix_memcpy(ws->tmp, g); - gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s); - return gsl_linalg_LU_det(ws->tmp, s); + double result; + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, g); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + result = gsl_linalg_LU_det(tmp, s); + + releaseTempMatrices(ws, 1); + return result; } int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws) @@ -168,9 +199,12 @@ int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws) int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) { - gsl_matrix_memcpy(ws->stack[0], g); + gsl_matrix *g_ = getTempMatrix(ws); + int success = 0; + + gsl_matrix_memcpy(g_, g); gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); - int r = gsl_eigen_nonsymmv(ws->stack[0], ws->eval_complex, ws->evec_complex, 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); @@ -181,21 +215,27 @@ int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) real = 0; if(!real) - return 0; // non-real eigenvalues! + goto eigenvectors_out; 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))); - return 1; + success = 1; + +eigenvectors_out: + releaseTempMatrices(ws, 1); + return success; } // only fills in the real eigenvectors and returns their count int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) { - gsl_matrix_memcpy(ws->stack[0], g); + gsl_matrix *g_ = getTempMatrix(ws); + + gsl_matrix_memcpy(g_, g); gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); - int r = gsl_eigen_nonsymmv(ws->stack[0], ws->eval_complex, ws->evec_complex, 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); @@ -210,23 +250,30 @@ int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) } } + releaseTempMatrices(ws, 1); return real; } 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); + gsl_matrix *g_ = getTempMatrix(ws); + + gsl_matrix_memcpy(g_, g); + int r = gsl_eigen_symmv (g_, eval, evec, ws->work_symmv); ERROR(r, "gsl_eigen_symmv failed!\n"); gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC); + + releaseTempMatrices(ws, 1); } // 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); + gsl_matrix *A_ = getTempMatrix(ws); + + gsl_matrix_memcpy(A_, A); + int r = gsl_eigen_symmv (A_, 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); @@ -243,9 +290,6 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws) *gsl_matrix_ptr(cob, i, j) *= sqrt(fabs(gsl_vector_get(ws->eval_real, i))); } + releaseTempMatrices(ws, 1); 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 5e5a946..59931fa 100644 --- a/linalg.h +++ b/linalg.h @@ -13,6 +13,9 @@ #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);} #define FCMP(x, y) gsl_fcmp(x, y, 1e-10) +#define MAX_TEMP_MATRICES 60000 +#define MAX_TEMP_VECTORS 100 + typedef struct _workspace { int n; gsl_eigen_nonsymmv_workspace *work_nonsymmv; @@ -22,9 +25,12 @@ typedef struct _workspace { gsl_matrix_complex *evec_complex; gsl_vector *eval_real; gsl_matrix *evec_real; - gsl_matrix *tmp; gsl_permutation *permutation; - gsl_matrix *stack[20]; + + gsl_matrix **tmp_mat; + int tmp_mat_used; + gsl_vector **tmp_vec; + int tmp_vec_used; } workspace_t; workspace_t *workspace_alloc(int n); @@ -46,4 +52,45 @@ int real_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); + +// matrix allocation stuff + +static gsl_matrix **getTempMatrices(workspace_t *ws, int n) +{ + ERROR(ws->tmp_mat_used + n > MAX_TEMP_MATRICES, "Ran out of temporary matrices. Consider increasing MAX_TEMP_MATRICES\n"); + int index = ws->tmp_mat_used; + ws->tmp_mat_used += n; + return ws->tmp_mat + index; +} + +static gsl_matrix *getTempMatrix(workspace_t *ws) +{ + return *getTempMatrices(ws, 1); +} + +static void releaseTempMatrices(workspace_t *ws, int n) +{ + ERROR(ws->tmp_mat_used - n < 0, "Released more matrices then in use\n"); + ws->tmp_mat_used -= n; +} + +static gsl_vector **getTempVectors(workspace_t *ws, int n) +{ + ERROR(ws->tmp_vec_used + n > MAX_TEMP_VECTORS, "Ran out of temporary vectors. Consider increasing MAX_TEMP_VECTORS\n"); + int index = ws->tmp_vec_used; + ws->tmp_vec_used += n; + return ws->tmp_vec + index; +} + +static gsl_vector *getTempVector(workspace_t *ws) +{ + return *getTempVectors(ws, 1); +} + +static void releaseTempVectors(workspace_t *ws, int n) +{ + ERROR(ws->tmp_vec_used - n < 0, "Released more vectors then in use\n"); + ws->tmp_vec_used -= n; +} + #endif diff --git a/main.c b/main.c index 6937919..72d6d55 100644 --- a/main.c +++ b/main.c @@ -43,15 +43,6 @@ void setupContext(DrawingContext *ctx) ctx->cartan = gsl_matrix_alloc(3, 3); ctx->cob = gsl_matrix_alloc(3, 3); ctx->ws = workspace_alloc(3); - ctx->tmp = malloc(MAX_TEMP_MATRICES*sizeof(gsl_matrix*)); - for(int i = 0; i < MAX_TEMP_MATRICES; i++) - ctx->tmp[i] = gsl_matrix_alloc(3, 3); - ctx->tmp_used = 0; - ctx->tmp_vec = malloc(MAX_TEMP_MATRICES*sizeof(gsl_vector*)); - for(int i = 0; i < MAX_TEMP_MATRICES; i++) - ctx->tmp_vec[i] = gsl_vector_alloc(3); - ctx->tmp_vec_used = 0; - } void destroyContext(DrawingContext *ctx) @@ -63,13 +54,6 @@ void destroyContext(DrawingContext *ctx) gsl_matrix_free(ctx->cob); workspace_free(ctx->ws); - - for(int i = 0; i < MAX_TEMP_MATRICES; i++) - gsl_matrix_free(ctx->tmp[i]); - free(ctx->tmp); - for(int i = 0; i < MAX_TEMP_VECTORS; i++) - gsl_vector_free(ctx->tmp_vec[i]); - free(ctx->tmp_vec); } void updateMatrices(DrawingContext *ctx) @@ -212,39 +196,6 @@ int main() updateMatrices(screen_context); computeLimitCurve(screen_context); - // do stuff -/* - DrawingContext *ctx = screen_context; - gsl_matrix *tmp = getTempMatrix(ctx); - gsl_matrix *ev = getTempMatrix(ctx); - gsl_matrix *evinv = getTempMatrix(ctx); - gsl_matrix **gen = getTempMatrices(ctx, 3); - gsl_matrix *transform = getTempMatrix(ctx); - vector_t eigenvectors[3][3]; - char words[3][4] = {"abc", "bca", "cab"}; - - initializeTriangleGenerators(gen, ctx->cartan); - - gsl_matrix_set_identity(tmp); - for(int j = 0; j < strlen(words[0]); j++) - multiply_right(tmp, gen[words[0][j]-'a'], ctx->ws); - real_eigenvectors(tmp, ev, ctx->ws); - invert(ev, evinv, ctx->ws); - - - LOOP(i) { - gsl_matrix_set_identity(tmp); - for(int j = 0; j < strlen(words[i]); j++) - multiply_right(tmp, gen[words[i][j]-'a'], ctx->ws); - real_eigenvectors(tmp, ev, ctx->ws); - LOOP(j) LOOP(k) eigenvectors[i][j].x[k] = gsl_matrix_get(ev, k, j); - } - - // solve the following: (x + y + z) (l) = lambda w - - return 0; -*/ - info = initCairo(0, KeyPressMask, 200, 200, "Triangle group"); if(!info) return 1; diff --git a/main.h b/main.h index 251a0c5..4fc4a3f 100644 --- a/main.h +++ b/main.h @@ -10,8 +10,6 @@ #define LOOP(i) for(int i = 0; i < 3; i++) #define NUM_GROUP_ELEMENTS 50000 -#define MAX_TEMP_MATRICES 60000 -#define MAX_TEMP_VECTORS 100 typedef struct { double x[3]; @@ -49,50 +47,8 @@ typedef struct { // temporary; matrices can only be freed from the top, but that's enough for us workspace_t *ws; - gsl_matrix **tmp; - int tmp_used; - gsl_vector **tmp_vec; - int tmp_vec_used; } DrawingContext; -static gsl_matrix **getTempMatrices(DrawingContext *ctx, int n) -{ - ERROR(ctx->tmp_used + n > MAX_TEMP_MATRICES, "Ran out of temporary matrices. Consider increasing MAX_TEMP_MATRICES\n"); - int index = ctx->tmp_used; - ctx->tmp_used += n; - return ctx->tmp + index; -} - -static gsl_matrix *getTempMatrix(DrawingContext *ctx) -{ - return *getTempMatrices(ctx, 1); -} - -static void releaseTempMatrices(DrawingContext *ctx, int n) -{ - ERROR(ctx->tmp_used - n < 0, "Released more matrices then in use\n"); - ctx->tmp_used -= n; -} - -static gsl_vector **getTempVectors(DrawingContext *ctx, int n) -{ - ERROR(ctx->tmp_vec_used + n > MAX_TEMP_VECTORS, "Ran out of temporary vectors. Consider increasing MAX_TEMP_VECTORS\n"); - int index = ctx->tmp_vec_used; - ctx->tmp_vec_used += n; - return ctx->tmp_vec + index; -} - -static gsl_vector *getTempVector(DrawingContext *ctx) -{ - return *getTempVectors(ctx, 1); -} - -static void releaseTempVectors(DrawingContext *ctx, int n) -{ - ERROR(ctx->tmp_vec_used - n < 0, "Released more vectors then in use\n"); - ctx->tmp_vec_used -= n; -} - // implemented in limit_set.c void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s); void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan);