229 lines
7.1 KiB
C
229 lines
7.1 KiB
C
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <gsl/gsl_math.h>
|
|
#include <gsl/gsl_eigen.h>
|
|
#include <gsl/gsl_blas.h>
|
|
#include <gsl/gsl_linalg.h>
|
|
#include <memory.h>
|
|
|
|
#include "linalg.h"
|
|
|
|
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
|
|
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
|
|
|
|
/*********************************************** temporary storage ********************************************************/
|
|
|
|
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->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);
|
|
|
|
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_matrix_free(workspace->tmp);
|
|
gsl_permutation_free(workspace->permutation);
|
|
for(int i = 0; i < 20; i++)
|
|
gsl_matrix_free(workspace->stack[i]);
|
|
}
|
|
|
|
/************************************************** basic operations ********************************************************/
|
|
|
|
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);
|
|
}
|
|
|
|
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *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);
|
|
}
|
|
|
|
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);
|
|
gsl_linalg_SV_decomp(ws->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));
|
|
}
|
|
|
|
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);
|
|
}
|
|
|
|
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 trace(gsl_matrix *g)
|
|
{
|
|
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;
|
|
gsl_matrix_memcpy(ws->tmp, g);
|
|
gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s);
|
|
return gsl_linalg_LU_det(ws->tmp, s);
|
|
}
|
|
|
|
int jordan_calc(gsl_matrix *g, double *evs, 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);
|
|
|
|
int real = 1;
|
|
for(int i = 0; i < ws->n; i++)
|
|
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
|
|
real = 0;
|
|
|
|
/* if(!real) {
|
|
printf("We have non-real eigenvalues!\n");
|
|
exit(1);
|
|
}*/
|
|
|
|
for(int i = 0; i < ws->n; i++) {
|
|
evs[2*i] = GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i));
|
|
evs[2*i+1] = GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i));
|
|
}
|
|
|
|
return real;
|
|
}
|
|
|
|
int 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(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 = 0;
|
|
for(int j = 0; j < ws->n; j++) {
|
|
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0) {
|
|
for(int i = 0; i < ws->n; i++)
|
|
gsl_matrix_set(evec_real, i, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
|
|
real++;
|
|
}
|
|
}
|
|
|
|
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);
|
|
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;
|
|
}
|