compute jordan projection for triangle rotation group
This commit is contained in:
Florian Stecker
commit
0c07494110
25
Makefile
Normal file
25
Makefile
Normal file
@ -0,0 +1,25 @@
|
|||||||
|
HEADERS=triangle.h linalg.h queue.h
|
||||||
|
|
||||||
|
#SPECIAL_OPTIONS=-O0 -g -D_DEBUG
|
||||||
|
#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
|
||||||
|
SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
|
||||||
|
#SPECIAL_OPTIONS=
|
||||||
|
|
||||||
|
OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
|
||||||
|
|
||||||
|
all: singular_values
|
||||||
|
|
||||||
|
singular_values: singular_values.o triangle.o linalg.o
|
||||||
|
gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas
|
||||||
|
|
||||||
|
singular_values.o: singular_values.c $(HEADERS)
|
||||||
|
gcc $(OPTIONS) -c singular_values.c
|
||||||
|
|
||||||
|
linalg.o: linalg.c $(HEADERS)
|
||||||
|
gcc $(OPTIONS) -c linalg.c
|
||||||
|
|
||||||
|
triangle.o: triangle.c $(HEADERS)
|
||||||
|
gcc $(OPTIONS) -c triangle.c
|
||||||
|
|
||||||
|
clean:
|
||||||
|
rm -f singular_values triangle.o linalg.o singular_values.o
|
||||||
353
linalg.c
Normal file
353
linalg.c
Normal file
@ -0,0 +1,353 @@
|
|||||||
|
#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 <gsl/gsl_complex_math.h>
|
||||||
|
#include <memory.h>
|
||||||
|
#include <math.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->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;
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
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 solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws)
|
||||||
|
{
|
||||||
|
int s;
|
||||||
|
gsl_matrix *tmp = getTempMatrix(ws);
|
||||||
|
|
||||||
|
gsl_matrix_memcpy(tmp, A);
|
||||||
|
gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
|
||||||
|
gsl_linalg_LU_solve(tmp, ws->permutation, b, result);
|
||||||
|
|
||||||
|
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, 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)
|
||||||
|
{
|
||||||
|
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_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_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, ...)
|
||||||
|
{
|
||||||
|
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 *tmp = getTempMatrix(ws);
|
||||||
|
|
||||||
|
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)
|
||||||
|
{
|
||||||
|
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;
|
||||||
|
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 eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
|
||||||
|
{
|
||||||
|
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(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);
|
||||||
|
|
||||||
|
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)
|
||||||
|
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)));
|
||||||
|
|
||||||
|
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 *g_ = getTempMatrix(ws);
|
||||||
|
|
||||||
|
gsl_matrix_memcpy(g_, g);
|
||||||
|
gsl_eigen_nonsymmv_params(0, 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);
|
||||||
|
|
||||||
|
int real = 0;
|
||||||
|
|
||||||
|
for(int i = 0; i < ws->n; i++) {
|
||||||
|
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) == 0) {// real
|
||||||
|
if(evec_real) {
|
||||||
|
for(int j = 0; j < ws->n; j++)
|
||||||
|
gsl_matrix_set(evec_real, j, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, j, i)));
|
||||||
|
}
|
||||||
|
real++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
releaseTempMatrices(ws, 1);
|
||||||
|
return real;
|
||||||
|
}
|
||||||
|
|
||||||
|
void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)
|
||||||
|
{
|
||||||
|
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 *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);
|
||||||
|
|
||||||
|
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)));
|
||||||
|
}
|
||||||
|
|
||||||
|
releaseTempMatrices(ws, 1);
|
||||||
|
return positive;
|
||||||
|
}
|
||||||
|
|
||||||
|
// computes a matrix in SL(3, R) which projectively transforms (e1, e2, e3, e1+e2+e3) to the 4 given vectors
|
||||||
|
void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws)
|
||||||
|
{
|
||||||
|
gsl_matrix *tmp = getTempMatrix(ws);
|
||||||
|
gsl_vector *coeff = getTempVector(ws);
|
||||||
|
int s;
|
||||||
|
double det, scale;
|
||||||
|
|
||||||
|
for(int i = 0; i < 3; i++)
|
||||||
|
for(int j = 0; j < 3; j++)
|
||||||
|
gsl_matrix_set(tmp, i, j, gsl_vector_get(vertices[j], i));
|
||||||
|
|
||||||
|
gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
|
||||||
|
gsl_linalg_LU_solve(tmp, ws->permutation, vertices[3], coeff);
|
||||||
|
det = gsl_linalg_LU_det(tmp, s);
|
||||||
|
|
||||||
|
for(int i = 0; i < 3; i++)
|
||||||
|
det *= gsl_vector_get(coeff, i);
|
||||||
|
scale = 1/cbrt(det);
|
||||||
|
|
||||||
|
for(int i = 0; i < 3; i++)
|
||||||
|
for(int j = 0; j < 3; j++)
|
||||||
|
gsl_matrix_set(result, i, j, scale*gsl_vector_get(vertices[j], i)*gsl_vector_get(coeff, j));
|
||||||
|
|
||||||
|
releaseTempMatrices(ws, 1);
|
||||||
|
releaseTempVectors(ws, 1);
|
||||||
|
}
|
||||||
|
|
||||||
|
void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws)
|
||||||
|
{
|
||||||
|
gsl_matrix *tmp = getTempMatrix(ws);
|
||||||
|
gsl_matrix *rot_basis = getTempMatrix(ws);
|
||||||
|
|
||||||
|
gsl_matrix_memcpy(tmp, rotation);
|
||||||
|
gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
|
||||||
|
int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
|
||||||
|
ERROR(r, "gsl_eigen_nonsymmv failed!\n");
|
||||||
|
|
||||||
|
double arg, minarg = 5; // greater than pi
|
||||||
|
int minidx;
|
||||||
|
for(int i = 0; i < 3; i++) {
|
||||||
|
arg = gsl_complex_arg(gsl_vector_complex_get(ws->eval_complex, i));
|
||||||
|
if(abs(arg) < minarg)
|
||||||
|
{
|
||||||
|
minidx = i;
|
||||||
|
minarg = abs(arg);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
ERROR(FCMP(minarg, 0.0) != 0, "rotation_frame() failed! No eigenvalue was 1.\n");
|
||||||
|
|
||||||
|
for(int i = 0; i < 3; i++) {
|
||||||
|
gsl_complex x = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+1)%3);
|
||||||
|
gsl_complex y = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+2)%3);
|
||||||
|
gsl_complex z = gsl_matrix_complex_get(ws->evec_complex, i, minidx);
|
||||||
|
gsl_matrix_set(result, i, 0, GSL_REAL(x)+GSL_REAL(y));
|
||||||
|
gsl_matrix_set(result, i, 1, GSL_IMAG(x)-GSL_IMAG(y));
|
||||||
|
gsl_matrix_set(result, i, 2, GSL_REAL(z));
|
||||||
|
}
|
||||||
|
|
||||||
|
releaseTempMatrices(ws, 2);
|
||||||
|
}
|
||||||
98
linalg.h
Normal file
98
linalg.h
Normal file
@ -0,0 +1,98 @@
|
|||||||
|
#ifndef LINALG_H
|
||||||
|
#define LINALG_H
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <stdarg.h>
|
||||||
|
#include <gsl/gsl_math.h>
|
||||||
|
#include <gsl/gsl_eigen.h>
|
||||||
|
#include <gsl/gsl_blas.h>
|
||||||
|
#include <gsl/gsl_linalg.h>
|
||||||
|
#include <memory.h>
|
||||||
|
|
||||||
|
#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 10
|
||||||
|
#define MAX_TEMP_VECTORS 10
|
||||||
|
|
||||||
|
typedef struct _workspace {
|
||||||
|
int n;
|
||||||
|
gsl_eigen_nonsymmv_workspace *work_nonsymmv;
|
||||||
|
gsl_eigen_symmv_workspace *work_symmv;
|
||||||
|
gsl_vector *work_sv;
|
||||||
|
gsl_vector_complex *eval_complex;
|
||||||
|
gsl_matrix_complex *evec_complex;
|
||||||
|
gsl_vector *eval_real;
|
||||||
|
gsl_matrix *evec_real;
|
||||||
|
gsl_permutation *permutation;
|
||||||
|
|
||||||
|
gsl_matrix **tmp_mat;
|
||||||
|
int tmp_mat_used;
|
||||||
|
gsl_vector **tmp_vec;
|
||||||
|
int tmp_vec_used;
|
||||||
|
} workspace_t;
|
||||||
|
|
||||||
|
workspace_t *workspace_alloc(int n);
|
||||||
|
void workspace_free(workspace_t *workspace);
|
||||||
|
void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws);
|
||||||
|
void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws);
|
||||||
|
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws);
|
||||||
|
void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out);
|
||||||
|
void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
|
||||||
|
void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
|
||||||
|
void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...);
|
||||||
|
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||||
|
void initialize(gsl_matrix *g, double *data, int x, int y);
|
||||||
|
void rotation_matrix(gsl_matrix *g, double *vector);
|
||||||
|
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||||
|
double trace(gsl_matrix *g);
|
||||||
|
double determinant(gsl_matrix *g, workspace_t *ws);
|
||||||
|
int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
|
||||||
|
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);
|
||||||
|
void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws);
|
||||||
|
void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, 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
|
||||||
46
queue.h
Normal file
46
queue.h
Normal file
@ -0,0 +1,46 @@
|
|||||||
|
#ifndef QUEUE_H
|
||||||
|
#define QUEUE_H
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
#define QUEUE_SIZE 1500000
|
||||||
|
|
||||||
|
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
|
||||||
|
#ifdef _DEBUG
|
||||||
|
#define LOG(msg, ...) fprintf(stderr, msg, ##__VA_ARGS__)
|
||||||
|
#else
|
||||||
|
#define LOG(msg, ...)
|
||||||
|
#endif
|
||||||
|
|
||||||
|
typedef struct {
|
||||||
|
unsigned int start;
|
||||||
|
unsigned int end;
|
||||||
|
int data[QUEUE_SIZE];
|
||||||
|
} queue_t;
|
||||||
|
|
||||||
|
static void queue_init(queue_t *q)
|
||||||
|
{
|
||||||
|
q->start = q->end = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
static void queue_put(queue_t *q, int x)
|
||||||
|
{
|
||||||
|
q->data[q->end++] = x;
|
||||||
|
q->end %= QUEUE_SIZE;
|
||||||
|
|
||||||
|
ERROR(q->start == q->end, "The queue is full! Increase QUEUE_SIZE\n");
|
||||||
|
}
|
||||||
|
|
||||||
|
static int queue_get(queue_t *q)
|
||||||
|
{
|
||||||
|
if(q->start == q->end)
|
||||||
|
return -1;
|
||||||
|
|
||||||
|
int result = q->data[q->start++];
|
||||||
|
q->start %= QUEUE_SIZE;
|
||||||
|
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
210
singular_values.c
Normal file
210
singular_values.c
Normal file
@ -0,0 +1,210 @@
|
|||||||
|
#include "triangle.h"
|
||||||
|
#include "linalg.h"
|
||||||
|
|
||||||
|
#define MAX_ELEMENTS 14000
|
||||||
|
//#define DRAW_PICTURE 1
|
||||||
|
|
||||||
|
void initialize_triangle_generators(workspace_t *ws, gsl_matrix **gen, double a1, double a2, double a3, double s, double t)
|
||||||
|
{
|
||||||
|
gsl_matrix **tmp;
|
||||||
|
tmp = getTempMatrices(ws, 6);
|
||||||
|
|
||||||
|
double rho1sqrt = sqrt(s*(s+2*cos(2*M_PI*a1))+1);
|
||||||
|
double rho2sqrt = sqrt(s*(s+2*cos(2*M_PI*a2))+1);
|
||||||
|
double rho3sqrt = sqrt(s*(s+2*cos(2*M_PI*a3))+1);
|
||||||
|
|
||||||
|
for(int i = 0; i < 6; i++)
|
||||||
|
gsl_matrix_set_zero(tmp[i]);
|
||||||
|
|
||||||
|
gsl_matrix_set(tmp[0], 0, 0, 1);
|
||||||
|
gsl_matrix_set(tmp[0], 1, 0, -rho3sqrt/t);
|
||||||
|
gsl_matrix_set(tmp[0], 2, 0, -rho2sqrt*t);
|
||||||
|
gsl_matrix_set(tmp[0], 1, 1, -1);
|
||||||
|
gsl_matrix_set(tmp[0], 2, 2, -1);
|
||||||
|
|
||||||
|
gsl_matrix_set(tmp[1], 0, 0, -1);
|
||||||
|
gsl_matrix_set(tmp[1], 0, 1, -rho3sqrt*t);
|
||||||
|
gsl_matrix_set(tmp[1], 1, 1, 1);
|
||||||
|
gsl_matrix_set(tmp[1], 2, 1, -rho1sqrt/t);
|
||||||
|
gsl_matrix_set(tmp[1], 2, 2, -1);
|
||||||
|
|
||||||
|
gsl_matrix_set(tmp[2], 0, 0, -1);
|
||||||
|
gsl_matrix_set(tmp[2], 1, 1, -1);
|
||||||
|
gsl_matrix_set(tmp[2], 0, 2, -rho2sqrt/t);
|
||||||
|
gsl_matrix_set(tmp[2], 1, 2, -rho1sqrt*t);
|
||||||
|
gsl_matrix_set(tmp[2], 2, 2, 1);
|
||||||
|
|
||||||
|
gsl_matrix_set(tmp[3], 0, 0, 1);
|
||||||
|
gsl_matrix_set(tmp[3], 1, 1, s);
|
||||||
|
gsl_matrix_set(tmp[3], 2, 2, 1/s);
|
||||||
|
|
||||||
|
gsl_matrix_set(tmp[4], 0, 0, 1/s);
|
||||||
|
gsl_matrix_set(tmp[4], 1, 1, 1);
|
||||||
|
gsl_matrix_set(tmp[4], 2, 2, s);
|
||||||
|
|
||||||
|
gsl_matrix_set(tmp[5], 0, 0, s);
|
||||||
|
gsl_matrix_set(tmp[5], 1, 1, 1/s);
|
||||||
|
gsl_matrix_set(tmp[5], 2, 2, 1);
|
||||||
|
|
||||||
|
multiply_many(ws, gen[0], 3, tmp[2], tmp[3], tmp[1]);
|
||||||
|
multiply_many(ws, gen[1], 3, tmp[0], tmp[4], tmp[2]);
|
||||||
|
multiply_many(ws, gen[2], 3, tmp[1], tmp[5], tmp[0]);
|
||||||
|
|
||||||
|
invert(gen[0], gen[3], ws);
|
||||||
|
invert(gen[1], gen[4], ws);
|
||||||
|
invert(gen[2], gen[5], ws);
|
||||||
|
|
||||||
|
releaseTempMatrices(ws, 6);
|
||||||
|
}
|
||||||
|
|
||||||
|
char *print_word(groupelement_t *g, char *str)
|
||||||
|
{
|
||||||
|
int i = g->length - 1;
|
||||||
|
|
||||||
|
str[g->length] = 0;
|
||||||
|
while(g->parent) {
|
||||||
|
str[i--] = 'a' + g->letter;
|
||||||
|
g = g->parent;
|
||||||
|
}
|
||||||
|
|
||||||
|
return str;
|
||||||
|
}
|
||||||
|
|
||||||
|
int main(int argc, char *argv[])
|
||||||
|
{
|
||||||
|
groupelement_t *group;
|
||||||
|
workspace_t *ws;
|
||||||
|
gsl_matrix **matrices;
|
||||||
|
gsl_matrix *gen[6];
|
||||||
|
gsl_matrix *tmp, *tmp2, *coxeter;
|
||||||
|
double mu[6], tr, trinv;
|
||||||
|
char buf[100];
|
||||||
|
|
||||||
|
/*
|
||||||
|
if(argc < 2) {
|
||||||
|
fprintf(stderr, "Too few arguments!\n");
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
// allocate stuff
|
||||||
|
group = malloc(MAX_ELEMENTS*sizeof(groupelement_t));
|
||||||
|
ws = workspace_alloc(3);
|
||||||
|
matrices = malloc(MAX_ELEMENTS*sizeof(gsl_matrix*));
|
||||||
|
for(int i = 0; i < MAX_ELEMENTS; i++)
|
||||||
|
matrices[i] = gsl_matrix_alloc(3, 3);
|
||||||
|
for(int i = 0; i < 6; i++)
|
||||||
|
gen[i] = gsl_matrix_alloc(3, 3);
|
||||||
|
tmp = gsl_matrix_alloc(3, 3);
|
||||||
|
tmp2 = gsl_matrix_alloc(3, 3);
|
||||||
|
coxeter = gsl_matrix_alloc(3, 3);
|
||||||
|
|
||||||
|
#ifdef DRAW_PICTURE
|
||||||
|
int width = 800;
|
||||||
|
int height = 800;
|
||||||
|
printf("P3\n%d %d\n255\n", width, height);
|
||||||
|
|
||||||
|
for(int y = 0; y < height; y++) {
|
||||||
|
for(int x = 0; x < width; x++) {
|
||||||
|
double s = exp(((double)x/width-0.5)*4);
|
||||||
|
double t = exp(((double)y/height-0.5)*4);
|
||||||
|
#else
|
||||||
|
double s = atof(argv[1]);
|
||||||
|
double t = atof(argv[2]);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
// the real action
|
||||||
|
generate_triangle_group(group, MAX_ELEMENTS, 5, 6, 7);
|
||||||
|
initialize_triangle_generators(ws, gen, 1.0/5.0, 1.0/6.0, 1.0/7.0, s, t);
|
||||||
|
|
||||||
|
gsl_matrix_set_identity(matrices[0]);
|
||||||
|
for(int i = 1; i < MAX_ELEMENTS; i++) {
|
||||||
|
if(group[i].length % 2 != 0)
|
||||||
|
continue;
|
||||||
|
|
||||||
|
int parent = group[i].parent->id;
|
||||||
|
int grandparent = group[i].parent->parent->id;
|
||||||
|
int letter;
|
||||||
|
|
||||||
|
if(group[parent].letter == 0 && group[i].letter == 1)
|
||||||
|
letter = 3; // p = ab
|
||||||
|
else if(group[parent].letter == 1 && group[i].letter == 2)
|
||||||
|
letter = 4; // q = bc
|
||||||
|
else if(group[parent].letter == 2 && group[i].letter == 0)
|
||||||
|
letter = 5; // r = ca
|
||||||
|
else if(group[parent].letter == 1 && group[i].letter == 0)
|
||||||
|
letter = 0; // p^{-1} = ba
|
||||||
|
else if(group[parent].letter == 2 && group[i].letter == 1)
|
||||||
|
letter = 1; // q^{-1} = cb
|
||||||
|
else if(group[parent].letter == 0 && group[i].letter == 2)
|
||||||
|
letter = 2; // r^{-1} = ac
|
||||||
|
|
||||||
|
multiply(matrices[grandparent], gen[letter], matrices[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
double excentricity;
|
||||||
|
double max_excentricity = 0;
|
||||||
|
int max_excentricity_index = 0;
|
||||||
|
|
||||||
|
for(int i = 0; i < MAX_ELEMENTS; i++) {
|
||||||
|
if(group[i].length % 2 != 0)
|
||||||
|
continue;
|
||||||
|
|
||||||
|
// jordan projection
|
||||||
|
gsl_matrix_memcpy(tmp, matrices[i]);
|
||||||
|
gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
|
||||||
|
int r = gsl_eigen_nonsymmv(tmp, 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_evs = 0;
|
||||||
|
|
||||||
|
for(int j = 0; j < 3; j++)
|
||||||
|
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0)
|
||||||
|
real_evs++;
|
||||||
|
|
||||||
|
if(real_evs != 3)
|
||||||
|
continue;
|
||||||
|
|
||||||
|
mu[0] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 0)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1)));
|
||||||
|
mu[1] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 2)));
|
||||||
|
|
||||||
|
excentricity = log(mu[1])/log(mu[0]);
|
||||||
|
if(excentricity > max_excentricity + 1e-3) {
|
||||||
|
max_excentricity = excentricity;
|
||||||
|
max_excentricity_index = i;
|
||||||
|
}
|
||||||
|
|
||||||
|
#ifndef DRAW_PICTURE
|
||||||
|
printf("%d %f %f %f %f 0x0000ff %d %s\n", group[i].length, mu[0], mu[1], tr, trinv, i, print_word(&group[i], buf));
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
#ifdef DRAW_PICTURE
|
||||||
|
int color = (int)((max_excentricity-1)/(max_excentricity+4)*255);
|
||||||
|
if(color < 0)
|
||||||
|
color = 0;
|
||||||
|
if(color > 255)
|
||||||
|
color = 255;
|
||||||
|
printf("%d %d %d\n", color, color, color);
|
||||||
|
}
|
||||||
|
fprintf(stderr, "y = %d\n", y);
|
||||||
|
}
|
||||||
|
#else
|
||||||
|
fprintf(stderr, "max = %f %s\n", max_excentricity, print_word(&group[max_excentricity_index], buf));
|
||||||
|
#endif
|
||||||
|
|
||||||
|
// free stuff
|
||||||
|
for(int i = 0; i < MAX_ELEMENTS; i++)
|
||||||
|
gsl_matrix_free(matrices[i]);
|
||||||
|
for(int i = 0; i < 6; i++)
|
||||||
|
gsl_matrix_free(gen[i]);
|
||||||
|
gsl_matrix_free(tmp);
|
||||||
|
gsl_matrix_free(tmp2);
|
||||||
|
gsl_matrix_free(coxeter);
|
||||||
|
free(matrices);
|
||||||
|
free(group);
|
||||||
|
workspace_free(ws);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
23
singular_values.plt
Normal file
23
singular_values.plt
Normal file
@ -0,0 +1,23 @@
|
|||||||
|
if(!exists("logt")) logt = log(1.5)
|
||||||
|
if(!exists("logs")) logs = log(1.0)
|
||||||
|
|
||||||
|
file = sprintf("< ./singular_values %f %f", exp(logs), exp(logt))
|
||||||
|
title = sprintf("s = %f, t = %f", exp(logs), exp(logt))
|
||||||
|
# print title
|
||||||
|
|
||||||
|
set zeroaxis
|
||||||
|
set size square
|
||||||
|
set xrange [0:25]
|
||||||
|
set yrange [0:25]
|
||||||
|
set grid
|
||||||
|
set parametric
|
||||||
|
|
||||||
|
plot file using (log($2)):(log($3)) w p pt 7 ps 0.5 lc 1 t title
|
||||||
|
# plot file using 2:3 w p pt 7 ps 0.5 lc 1 t title
|
||||||
|
|
||||||
|
pause mouse keypress
|
||||||
|
if(MOUSE_KEY == 60) logt=logt-0.02
|
||||||
|
if(MOUSE_KEY == 62) logt=logt+0.02
|
||||||
|
if(MOUSE_KEY == 44) logs=logs-0.02
|
||||||
|
if(MOUSE_KEY == 46) logs=logs+0.02
|
||||||
|
if(MOUSE_KEY != 113) reread
|
||||||
72
triangle.c
Normal file
72
triangle.c
Normal file
@ -0,0 +1,72 @@
|
|||||||
|
#include <stdio.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <memory.h>
|
||||||
|
|
||||||
|
#include "triangle.h"
|
||||||
|
#include "queue.h"
|
||||||
|
|
||||||
|
static groupelement_t *find_adj(groupelement_t *cur, int gen, int other, int order)
|
||||||
|
{
|
||||||
|
groupelement_t *next = cur->adj[other];
|
||||||
|
for(int i = 0; i < order - 1; i++) {
|
||||||
|
if(!next)
|
||||||
|
break;
|
||||||
|
next = next->adj[gen];
|
||||||
|
if(!next)
|
||||||
|
break;
|
||||||
|
next = next->adj[other];
|
||||||
|
}
|
||||||
|
return next;
|
||||||
|
}
|
||||||
|
|
||||||
|
int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3)
|
||||||
|
{
|
||||||
|
queue_t q;
|
||||||
|
int id, n;
|
||||||
|
groupelement_t *cur, *inv;
|
||||||
|
|
||||||
|
int k[3] = {k1, k2, k3};
|
||||||
|
|
||||||
|
memset(group, 0, size*sizeof(groupelement_t));
|
||||||
|
group[0].letter = -1;
|
||||||
|
n = 1;
|
||||||
|
queue_init(&q);
|
||||||
|
queue_put(&q, 0);
|
||||||
|
|
||||||
|
while((id = queue_get(&q)) != -1) {
|
||||||
|
cur = &group[id];
|
||||||
|
|
||||||
|
for(int gen = 0; gen < 3; gen++) {
|
||||||
|
if(!cur->adj[gen])
|
||||||
|
cur->adj[gen] = find_adj(cur, gen, (gen+2)%3, k[(gen+1)%3]);
|
||||||
|
|
||||||
|
if(!cur->adj[gen])
|
||||||
|
cur->adj[gen] = find_adj(cur, gen, (gen+1)%3, k[(gen+2)%3]);
|
||||||
|
|
||||||
|
if(!cur->adj[gen] && n < size) {
|
||||||
|
cur->adj[gen] = &group[n];
|
||||||
|
cur->adj[gen]->id = n;
|
||||||
|
cur->adj[gen]->parent = cur;
|
||||||
|
cur->adj[gen]->letter = gen;
|
||||||
|
cur->adj[gen]->length = cur->length + 1;
|
||||||
|
queue_put(&q, n);
|
||||||
|
n++;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(cur->adj[gen])
|
||||||
|
cur->adj[gen]->adj[gen] = cur;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int i = 0; i < size; i++) {
|
||||||
|
cur = &group[i];
|
||||||
|
inv = &group[0];
|
||||||
|
while(cur->parent && inv) {
|
||||||
|
inv = inv->adj[cur->letter];
|
||||||
|
cur = cur->parent;
|
||||||
|
}
|
||||||
|
group[i].inverse = inv;
|
||||||
|
}
|
||||||
|
|
||||||
|
return n;
|
||||||
|
}
|
||||||
21
triangle.h
Normal file
21
triangle.h
Normal file
@ -0,0 +1,21 @@
|
|||||||
|
#ifndef TRIANGLE_H
|
||||||
|
#define TRIANGLE_H
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <memory.h>
|
||||||
|
|
||||||
|
#include "queue.h"
|
||||||
|
|
||||||
|
typedef struct _groupelement {
|
||||||
|
int id;
|
||||||
|
int length;
|
||||||
|
struct _groupelement *adj[3];
|
||||||
|
struct _groupelement *parent;
|
||||||
|
struct _groupelement *inverse;
|
||||||
|
int letter;
|
||||||
|
} groupelement_t;
|
||||||
|
|
||||||
|
int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);
|
||||||
|
|
||||||
|
#endif
|
||||||
Loading…
Reference in New Issue
Block a user