From 8146f0ee6005f8b19f3906642ab4c7181aaf4b0d Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Fri, 8 Feb 2019 13:35:46 +0100 Subject: [PATCH] fix indentation --- linalg.c | 180 +++++++++++++++++++++++++++---------------------------- 1 file changed, 90 insertions(+), 90 deletions(-) diff --git a/linalg.c b/linalg.c index 97a7be8..704ca7c 100644 --- a/linalg.c +++ b/linalg.c @@ -15,75 +15,75 @@ 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); - result->eval_real = gsl_vector_alloc(n); - result->evec_real = gsl_matrix_alloc(n, n); - result->permutation = gsl_permutation_alloc(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); + result->eval_real = gsl_vector_alloc(n); + result->evec_real = gsl_matrix_alloc(n, n); + result->permutation = gsl_permutation_alloc(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; + 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; } 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); - gsl_vector_free(workspace->eval_real); - gsl_matrix_free(workspace->evec_real); - gsl_permutation_free(workspace->permutation); + 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); + gsl_vector_free(workspace->eval_real); + gsl_matrix_free(workspace->evec_real); + gsl_permutation_free(workspace->permutation); - 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); + 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 ********************************************************/ void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws) { - int s; - gsl_matrix *tmp = getTempMatrix(ws); + int s; + 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); + gsl_matrix_memcpy(tmp, in); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + 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) { - gsl_matrix *tmp = getTempMatrix(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, tmp); // in * conjugator^{-1} - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, tmp, 0.0, out); + 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, conjugator, tmp, 0.0, out); - releaseTempMatrices(ws, 1); + releaseTempMatrices(ws, 1); } 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) @@ -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, ...) { va_list args; - va_start(args, n); + va_start(args, n); 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); } - va_end(args); + va_end(args); } 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) { - gsl_matrix_view view = gsl_matrix_view_array(data, x, y); - gsl_matrix_memcpy(g, &view.matrix); + gsl_matrix_view view = gsl_matrix_view_array(data, x, y); + gsl_matrix_memcpy(g, &view.matrix); } void rotation_matrix(gsl_matrix *g, double *vector) { - double normalized[3]; - double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]); - for(int i = 0; i < 3; i++) - normalized[i] = vector[i] / norm; - gsl_matrix_set_identity(g); - gsl_matrix_set(g, 0, 0, cos(norm)); - gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]); - 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, 1, cos(norm)); - 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, 1, +sin(norm) * normalized[0]); - gsl_matrix_set(g, 2, 2, cos(norm)); - for(int i = 0; i < 3; i++) - for(int j = 0; j < 3; j++) - g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j]; + double normalized[3]; + double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]); + for(int i = 0; i < 3; i++) + normalized[i] = vector[i] / norm; + gsl_matrix_set_identity(g); + gsl_matrix_set(g, 0, 0, cos(norm)); + gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]); + 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, 1, cos(norm)); + 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, 1, +sin(norm) * normalized[0]); + gsl_matrix_set(g, 2, 2, cos(norm)); + for(int i = 0; i < 3; i++) + for(int j = 0; j < 3; j++) + g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j]; } 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) { - int s; - double result; - gsl_matrix *tmp = getTempMatrix(ws); + int 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); + 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; + releaseTempMatrices(ws, 1); + return result; } int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws) { - gsl_eigen_nonsymmv_params(1, 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_params(1, 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); - int real = 1; - for(int i = 0; i < 4; i++) - if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) - real = 0; + int real = 1; + for(int i = 0; i < 4; i++) + if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) + real = 0; - if(!real) - return 0; // non-real eigenvalues! + if(!real) + return 0; // non-real eigenvalues! - 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))); - } + 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))); + } - return 1; + return 1; } 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; 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))); + 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))); success = 1;