florian
/
triangle_group_limit_set
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

fix indentation

This commit is contained in:
Florian Stecker 2019-02-08 13:35:46 +01:00
parent 870ae7d2d2
commit 8146f0ee60
1 changed files with 90 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

180
linalg.c
View File

@ -15,75 +15,75 @@
workspace_t *workspace_alloc(int n) workspace_t *workspace_alloc(int n)
{ {
workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t)); workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t));
result->n = n; result->n = n;
result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n); result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n);
result->work_symmv = gsl_eigen_symmv_alloc(n); result->work_symmv = gsl_eigen_symmv_alloc(n);
result->work_sv = gsl_vector_alloc(n); result->work_sv = gsl_vector_alloc(n);
result->eval_complex = gsl_vector_complex_alloc(n); result->eval_complex = gsl_vector_complex_alloc(n);
result->evec_complex = gsl_matrix_complex_alloc(n, n); result->evec_complex = gsl_matrix_complex_alloc(n, n);
result->eval_real = gsl_vector_alloc(n); result->eval_real = gsl_vector_alloc(n);
result->evec_real = gsl_matrix_alloc(n, n); result->evec_real = gsl_matrix_alloc(n, n);
result->permutation = gsl_permutation_alloc(n); result->permutation = gsl_permutation_alloc(n);
result->tmp_mat = malloc(MAX_TEMP_MATRICES*sizeof(gsl_matrix*)); result->tmp_mat = malloc(MAX_TEMP_MATRICES*sizeof(gsl_matrix*));
for(int i = 0; i < MAX_TEMP_MATRICES; i++) for(int i = 0; i < MAX_TEMP_MATRICES; i++)
result->tmp_mat[i] = gsl_matrix_alloc(3, 3); result->tmp_mat[i] = gsl_matrix_alloc(3, 3);
result->tmp_mat_used = 0; result->tmp_mat_used = 0;
result->tmp_vec = malloc(MAX_TEMP_MATRICES*sizeof(gsl_vector*)); result->tmp_vec = malloc(MAX_TEMP_MATRICES*sizeof(gsl_vector*));
for(int i = 0; i < MAX_TEMP_MATRICES; i++) for(int i = 0; i < MAX_TEMP_MATRICES; i++)
result->tmp_vec[i] = gsl_vector_alloc(3); result->tmp_vec[i] = gsl_vector_alloc(3);
result->tmp_vec_used = 0; result->tmp_vec_used = 0;
return result; return result;
} }
void workspace_free(workspace_t *workspace) void workspace_free(workspace_t *workspace)
{ {
gsl_eigen_nonsymmv_free(workspace->work_nonsymmv); gsl_eigen_nonsymmv_free(workspace->work_nonsymmv);
gsl_eigen_symmv_free(workspace->work_symmv); gsl_eigen_symmv_free(workspace->work_symmv);
gsl_vector_free(workspace->work_sv); gsl_vector_free(workspace->work_sv);
gsl_vector_complex_free(workspace->eval_complex); gsl_vector_complex_free(workspace->eval_complex);
gsl_matrix_complex_free(workspace->evec_complex); gsl_matrix_complex_free(workspace->evec_complex);
gsl_vector_free(workspace->eval_real); gsl_vector_free(workspace->eval_real);
gsl_matrix_free(workspace->evec_real); gsl_matrix_free(workspace->evec_real);
gsl_permutation_free(workspace->permutation); gsl_permutation_free(workspace->permutation);
for(int i = 0; i < MAX_TEMP_MATRICES; i++) for(int i = 0; i < MAX_TEMP_MATRICES; i++)
gsl_matrix_free(workspace->tmp_mat[i]); gsl_matrix_free(workspace->tmp_mat[i]);
free(workspace->tmp_mat); free(workspace->tmp_mat);
for(int i = 0; i < MAX_TEMP_VECTORS; i++) for(int i = 0; i < MAX_TEMP_VECTORS; i++)
gsl_vector_free(workspace->tmp_vec[i]); gsl_vector_free(workspace->tmp_vec[i]);
free(workspace->tmp_vec); free(workspace->tmp_vec);
} }
/************************************************** basic operations ********************************************************/ /************************************************** basic operations ********************************************************/
void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws) void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws)
{ {
int s; int s;
gsl_matrix *tmp = getTempMatrix(ws); gsl_matrix *tmp = getTempMatrix(ws);
gsl_matrix_memcpy(tmp, in); gsl_matrix_memcpy(tmp, in);
gsl_linalg_LU_decomp(tmp, ws->permutation, &s); gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
gsl_linalg_LU_invert(tmp, ws->permutation, out); gsl_linalg_LU_invert(tmp, ws->permutation, out);
releaseTempMatrices(ws, 1); releaseTempMatrices(ws, 1);
} }
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws) void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws)
{ {
gsl_matrix *tmp = getTempMatrix(ws); gsl_matrix *tmp = getTempMatrix(ws);
invert(conjugator, out, ws); // use out to temporarily store inverse conjugator invert(conjugator, out, ws); // use out to temporarily store inverse conjugator
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, tmp); // in * conjugator^{-1} 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); gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, tmp, 0.0, out);
releaseTempMatrices(ws, 1); releaseTempMatrices(ws, 1);
} }
void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out) void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out)
{ {
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, 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) void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
@ -105,7 +105,7 @@ void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...) void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...)
{ {
va_list args; va_list args;
va_start(args, n); va_start(args, n);
gsl_matrix_set_identity(out); gsl_matrix_set_identity(out);
@ -114,7 +114,7 @@ void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...)
multiply_right(out, cur, ws); multiply_right(out, cur, ws);
} }
va_end(args); va_end(args);
} }
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws) void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
@ -132,69 +132,69 @@ void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
void initialize(gsl_matrix *g, double *data, int x, int y) void initialize(gsl_matrix *g, double *data, int x, int y)
{ {
gsl_matrix_view view = gsl_matrix_view_array(data, x, y); gsl_matrix_view view = gsl_matrix_view_array(data, x, y);
gsl_matrix_memcpy(g, &view.matrix); gsl_matrix_memcpy(g, &view.matrix);
} }
void rotation_matrix(gsl_matrix *g, double *vector) void rotation_matrix(gsl_matrix *g, double *vector)
{ {
double normalized[3]; double normalized[3];
double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]); double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
for(int i = 0; i < 3; i++) for(int i = 0; i < 3; i++)
normalized[i] = vector[i] / norm; normalized[i] = vector[i] / norm;
gsl_matrix_set_identity(g); gsl_matrix_set_identity(g);
gsl_matrix_set(g, 0, 0, cos(norm)); gsl_matrix_set(g, 0, 0, cos(norm));
gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]); gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]);
gsl_matrix_set(g, 0, 2, +sin(norm) * normalized[1]); gsl_matrix_set(g, 0, 2, +sin(norm) * normalized[1]);
gsl_matrix_set(g, 1, 0, +sin(norm) * normalized[2]); gsl_matrix_set(g, 1, 0, +sin(norm) * normalized[2]);
gsl_matrix_set(g, 1, 1, cos(norm)); gsl_matrix_set(g, 1, 1, cos(norm));
gsl_matrix_set(g, 1, 2, -sin(norm) * normalized[0]); gsl_matrix_set(g, 1, 2, -sin(norm) * normalized[0]);
gsl_matrix_set(g, 2, 0, -sin(norm) * normalized[1]); gsl_matrix_set(g, 2, 0, -sin(norm) * normalized[1]);
gsl_matrix_set(g, 2, 1, +sin(norm) * normalized[0]); gsl_matrix_set(g, 2, 1, +sin(norm) * normalized[0]);
gsl_matrix_set(g, 2, 2, cos(norm)); gsl_matrix_set(g, 2, 2, cos(norm));
for(int i = 0; i < 3; i++) for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++) for(int j = 0; j < 3; j++)
g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j]; g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j];
} }
double trace(gsl_matrix *g) double trace(gsl_matrix *g)
{ {
return gsl_matrix_get(g, 0, 0) + gsl_matrix_get(g, 1, 1) + gsl_matrix_get(g, 2, 2); return gsl_matrix_get(g, 0, 0) + gsl_matrix_get(g, 1, 1) + gsl_matrix_get(g, 2, 2);
} }
double determinant(gsl_matrix *g, workspace_t *ws) double determinant(gsl_matrix *g, workspace_t *ws)
{ {
int s; int s;
double result; double result;
gsl_matrix *tmp = getTempMatrix(ws); gsl_matrix *tmp = getTempMatrix(ws);
gsl_matrix_memcpy(tmp, g); gsl_matrix_memcpy(tmp, g);
gsl_linalg_LU_decomp(tmp, ws->permutation, &s); gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
result = gsl_linalg_LU_det(tmp, s); result = gsl_linalg_LU_det(tmp, s);
releaseTempMatrices(ws, 1); releaseTempMatrices(ws, 1);
return result; return result;
} }
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws) int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
{ {
gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv); gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv);
gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, 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_ABS_DESC);
int real = 1; int real = 1;
for(int i = 0; i < 4; i++) for(int i = 0; i < 4; i++)
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
real = 0; real = 0;
if(!real) if(!real)
return 0; // non-real eigenvalues! return 0; // non-real eigenvalues!
for(int i = 0; i < ws->n - 1; i++) { 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))); mu[i] = log(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i+1)));
} }
return 1; return 1;
} }
int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
@ -218,8 +218,8 @@ int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
goto eigenvectors_out; goto eigenvectors_out;
for(int i = 0; i < ws->n; i++) for(int i = 0; i < ws->n; i++)
for(int j = 0; j < ws->n; j++) 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))); gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
success = 1; success = 1;