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

new matrix allocation mechanism

Browse Source
This commit is contained in:
Florian Stecker
2019-02-08 13:35:02 +01:00
parent 3d8378aa16
commit 870ae7d2d2
6 changed files with 152 additions and 154 deletions
Download Patch File Download Diff File Expand all files Collapse all files

116
linalg.c
View File

@@ -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;
}