Merge branch 'better_boxes'
This commit is contained in:
Florian Stecker
commit
89cbec2860
14
Makefile
14
Makefile
@ -1,13 +1,13 @@
|
||||
HEADERS=triangle.h linalg.h queue.h
|
||||
|
||||
#SPECIAL_OPTIONS=-O0 -g -D_DEBUG
|
||||
SPECIAL_OPTIONS=-O0 -g -D_DEBUG
|
||||
#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
|
||||
SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
|
||||
#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 element_path limit_set
|
||||
all: singular_values element_path limit_set hyperbolic
|
||||
|
||||
singular_values: singular_values.o triangle.o linalg.o
|
||||
gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas
|
||||
@ -18,6 +18,9 @@ element_path: element_path.o linalg.o
|
||||
limit_set: limit_set.o linalg.o triangle.o
|
||||
gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o -lm -lgsl -lcblas
|
||||
|
||||
hyperbolic: hyperbolic.o triangle.o linalg.o
|
||||
gcc $(OPTIONS) -o hyperbolic hyperbolic.o triangle.o linalg.o -lm -lgsl -lcblas
|
||||
|
||||
singular_values.o: singular_values.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c singular_values.c
|
||||
|
||||
@ -33,5 +36,8 @@ triangle.o: triangle.c $(HEADERS)
|
||||
limit_set.o: limit_set.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c limit_set.c
|
||||
|
||||
hyperbolic.o: hyperbolic.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c hyperbolic.c
|
||||
|
||||
clean:
|
||||
rm -f singular_values element_path limit_set triangle.o linalg.o singular_values.o element_path.o limit_set.o
|
||||
rm -f singular_values element_path limit_set hyperbolic triangle.o linalg.o singular_values.o element_path.o limit_set.o hyperbolic.o
|
||||
|
||||
265
hyperbolic.c
Normal file
265
hyperbolic.c
Normal file
@ -0,0 +1,265 @@
|
||||
#include <complex.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include "triangle.h"
|
||||
#include "linalg.h"
|
||||
|
||||
#define POINCARE 1
|
||||
#define LOOP(i) for(int i = 0; i < 3; i++)
|
||||
|
||||
void cartan_matrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
|
||||
{
|
||||
gsl_matrix_set(cartan, 0, 0, -2);
|
||||
gsl_matrix_set(cartan, 0, 1, 2*s*cos(a3));
|
||||
gsl_matrix_set(cartan, 0, 2, 2/s*cos(a2));
|
||||
|
||||
gsl_matrix_set(cartan, 1, 0, 2/s*cos(a3));
|
||||
gsl_matrix_set(cartan, 1, 1, -2);
|
||||
gsl_matrix_set(cartan, 1, 2, 2*s*cos(a1));
|
||||
|
||||
gsl_matrix_set(cartan, 2, 0, 2*s*cos(a2));
|
||||
gsl_matrix_set(cartan, 2, 1, 2/s*cos(a1));
|
||||
gsl_matrix_set(cartan, 2, 2, -2);
|
||||
}
|
||||
|
||||
void initialize_triangle_generators(gsl_matrix **gen, gsl_matrix *cartan)
|
||||
{
|
||||
LOOP(i) {
|
||||
gsl_matrix_set_identity(gen[i]);
|
||||
LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j);
|
||||
}
|
||||
}
|
||||
|
||||
typedef struct { double x[3]; } point;
|
||||
|
||||
point apply(gsl_matrix *g, point x)
|
||||
{
|
||||
point out;
|
||||
LOOP(i) out.x[i] = 0.0;
|
||||
LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(g, i, j) * x.x[j];
|
||||
return out;
|
||||
}
|
||||
|
||||
point incidence(point a, point b)
|
||||
{
|
||||
point c;
|
||||
LOOP(i) c.x[i] = a.x[(i+1)%3]*b.x[(i+2)%3] - a.x[(i+2)%3]*b.x[(i+1)%3];
|
||||
return c;
|
||||
}
|
||||
|
||||
double applyForm(point x, point y, gsl_matrix *form)
|
||||
{
|
||||
double out = 0.0;
|
||||
LOOP(i) LOOP(j) out += x.x[i] * gsl_matrix_get(form, i, j) * y.x[j];
|
||||
return out;
|
||||
}
|
||||
|
||||
void boundary(point x, point y, point *a, point *b, gsl_matrix *form)
|
||||
{
|
||||
double formX = applyForm(x, x, form);
|
||||
double formY = applyForm(y, y, form);
|
||||
double formXY = applyForm(x, y, form);
|
||||
double t1 = (- formXY + sqrt(formXY * formXY - formX * formY)) / formX;
|
||||
double t2 = (- formXY - sqrt(formXY * formXY - formX * formY)) / formX;
|
||||
LOOP(i) a->x[i] = t1 * x.x[i] + y.x[i];
|
||||
LOOP(i) b->x[i] = t2 * x.x[i] + y.x[i];
|
||||
}
|
||||
|
||||
point column(gsl_matrix *m, int j)
|
||||
{
|
||||
point out;
|
||||
LOOP(i) out.x[i] = gsl_matrix_get(m, i, j);
|
||||
return out;
|
||||
}
|
||||
|
||||
point row(gsl_matrix *m, int i)
|
||||
{
|
||||
point out;
|
||||
LOOP(j) out.x[j] = gsl_matrix_get(m, i, j);
|
||||
return out;
|
||||
}
|
||||
|
||||
double coord(point p, int k, gsl_matrix *frame)
|
||||
{
|
||||
double tmp[3] = {0, 0, 0};
|
||||
double factor = 1.0;
|
||||
|
||||
LOOP(i) LOOP(j) tmp[i] += gsl_matrix_get(frame, i, j)*p.x[j];
|
||||
|
||||
#ifdef POINCARE
|
||||
double norm2 = (tmp[0]*tmp[0] + tmp[1]*tmp[1]) / (tmp[2]*tmp[2]);
|
||||
if(fabs(norm2) < 0.5) // just to avoid dividing by 0
|
||||
factor = 1.0 / (1.0 + sqrt(1.0 - norm2));
|
||||
else if(fabs(norm2) < 1.0)
|
||||
factor = (1.0 - sqrt(1.0 - norm2)) / norm2;
|
||||
else
|
||||
factor = 1.0;
|
||||
#endif
|
||||
|
||||
return factor*tmp[k]/tmp[2];
|
||||
}
|
||||
|
||||
void print_tex_header() {
|
||||
char *header =
|
||||
"\\documentclass{standalone}\n"
|
||||
"\\usepackage[utf8]{inputenc}\n"
|
||||
"\\usepackage{tikz}\n"
|
||||
"\\begin{document}\n"
|
||||
"\\begin{tikzpicture}[scale=5]\n";
|
||||
|
||||
printf("%s", header);
|
||||
}
|
||||
|
||||
void print_tex_footer() {
|
||||
char *footer =
|
||||
"\\draw[thick] (0,0) circle (1.0);\n"
|
||||
"\\end{tikzpicture}\n"
|
||||
"\\end{document}\n";
|
||||
|
||||
printf("%s", footer);
|
||||
}
|
||||
|
||||
char *poincare_arc(double x1, double x2, double y1, double y2, char *buffer)
|
||||
{
|
||||
double t = (1 - x1*y1 - x2*y2)/(x1*y2 - x2*y1)/2;
|
||||
double m1 = x1/2 + y1/2 + t*(y2 - x2); // center of circle
|
||||
double m2 = x2/2 + y2/2 + t*(x1 - y1);
|
||||
|
||||
double r = sqrt((x1-m1)*(x1-m1) + (x2-m2)*(x2-m2));
|
||||
double phix = atan2(x2-m2,x1-m1);
|
||||
double phiy = atan2(y2-m2,y1-m1);
|
||||
|
||||
if(phix - phiy > M_PI)
|
||||
phiy += 2*M_PI;
|
||||
else if(phiy - phix > M_PI)
|
||||
phix += 2*M_PI;
|
||||
|
||||
sprintf(buffer, "%f:%f:%f", phix/M_PI*180, phiy/M_PI*180, r);
|
||||
return buffer;
|
||||
}
|
||||
|
||||
void draw_triangle(point *p, gsl_matrix *frame, const char *arguments)
|
||||
{
|
||||
#ifdef POINCARE
|
||||
char buffer1[100], buffer2[100], buffer3[100];
|
||||
|
||||
double x1 = coord(p[0], 0, frame);
|
||||
double y1 = coord(p[0], 1, frame);
|
||||
double x2 = coord(p[1], 0, frame);
|
||||
double y2 = coord(p[1], 1, frame);
|
||||
double x3 = coord(p[2], 0, frame);
|
||||
double y3 = coord(p[2], 1, frame);
|
||||
|
||||
printf("\\draw[%s] (%f, %f) arc (%s) arc (%s) arc (%s);\n", arguments, x1, y1,
|
||||
poincare_arc(x1, y1, x2, y2, buffer1),
|
||||
poincare_arc(x2, y2, x3, y3, buffer2),
|
||||
poincare_arc(x3, y3, x1, y1, buffer3));
|
||||
|
||||
#else
|
||||
printf("\\draw[%s] (%f, %f) -- (%f, %f) -- (%f, %f) -- cycle;\n", arguments,
|
||||
coord(p[0], 0, frame), coord(p[0], 1, frame),
|
||||
coord(p[1], 0, frame), coord(p[1], 1, frame),
|
||||
coord(p[2], 0, frame), coord(p[2], 1, frame));
|
||||
#endif
|
||||
}
|
||||
|
||||
void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments)
|
||||
{
|
||||
char buffer[100];
|
||||
|
||||
double x1 = coord(p1, 0, frame);
|
||||
double y1 = coord(p1, 1, frame);
|
||||
double x2 = coord(p2, 0, frame);
|
||||
double y2 = coord(p2, 1, frame);
|
||||
|
||||
#ifdef POINCARE
|
||||
printf("\\draw[%s] (%f, %f) arc (%s);\n", arguments, x1, y1,
|
||||
poincare_arc(x1, y1, x2, y2, buffer));
|
||||
#else
|
||||
printf("\\draw[%s] (%f, %f) -- (%f, %f);\n", arguments, x1, y1, x2, y2);
|
||||
#endif
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
groupelement_t *group;
|
||||
gsl_matrix **matrices;
|
||||
gsl_matrix *cartan;
|
||||
gsl_matrix *gen[3];
|
||||
gsl_matrix *coxeter[3];
|
||||
gsl_matrix *coxeter_eigenvectors[3];
|
||||
gsl_matrix *frame;
|
||||
workspace_t *ws;
|
||||
int elements = 5000;
|
||||
int p = 3, q = 5, r = 7;
|
||||
|
||||
group = malloc(elements*sizeof(groupelement_t));
|
||||
matrices = malloc(elements*sizeof(gsl_matrix*));
|
||||
for(int i = 0; i < elements; i++)
|
||||
matrices[i] = gsl_matrix_alloc(3, 3);
|
||||
cartan = gsl_matrix_alloc(3, 3);
|
||||
frame = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) gen[i] = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) coxeter[i] = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) coxeter_eigenvectors[i] = gsl_matrix_alloc(3, 3);
|
||||
ws = workspace_alloc(3);
|
||||
|
||||
generate_triangle_group(group, elements, p, q, r);
|
||||
cartan_matrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0);
|
||||
initialize_triangle_generators(gen, cartan);
|
||||
int pos = diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form
|
||||
|
||||
gsl_matrix_set_identity(matrices[0]);
|
||||
for(int i = 1; i < elements; i++)
|
||||
multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
|
||||
|
||||
LOOP(i) multiply_many(ws, coxeter[i], 3, gen[i%3], gen[(i+1)%3], gen[(i+2)%3]);
|
||||
LOOP(i) eigenvectors(coxeter[i], coxeter_eigenvectors[i], ws);
|
||||
|
||||
/*
|
||||
for(int i = 0; i < elements; i++) {
|
||||
printf("%4d: ", i);
|
||||
for(groupelement_t *cur = &group[i]; cur->parent; cur = cur->parent)
|
||||
fputc(cur->letter+'a', stdout);
|
||||
fputc('\n', stdout);
|
||||
}
|
||||
|
||||
return 0;*/
|
||||
|
||||
point coxeter_attracting[3];
|
||||
point coxeter_repelling[3];
|
||||
point reflection_lines[3];
|
||||
point triangle_points[3];
|
||||
point transformed[3];
|
||||
point transformed2[3];
|
||||
|
||||
LOOP(i) coxeter_attracting[i] = column(coxeter_eigenvectors[i], 0);
|
||||
LOOP(i) coxeter_repelling[i] = column(coxeter_eigenvectors[i], 2);
|
||||
LOOP(i) reflection_lines[i] = row(cartan, i);
|
||||
LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]);
|
||||
|
||||
print_tex_header();
|
||||
|
||||
// int indices[10] = {0, 1, 6, 10, 30, 46, 124, 185, 484, 717};
|
||||
// int indices[10] = {0, 1, 4, 10, 22};
|
||||
|
||||
for(int k = 0; k < elements; k++) {
|
||||
if(group[k].length % 2)
|
||||
continue;
|
||||
|
||||
LOOP(i) transformed[i] = apply(matrices[k], triangle_points[i]);
|
||||
draw_triangle(transformed, frame, "black,fill=black!10,line width=0pt");
|
||||
}
|
||||
|
||||
for(int k = 0; k < elements; k++) {
|
||||
if(group[k].length % 2)
|
||||
continue;
|
||||
|
||||
LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]);
|
||||
LOOP(i) transformed2[i] = apply(matrices[k], coxeter_repelling[i]);
|
||||
LOOP(i) draw_line(transformed[i], transformed2[i], frame, "red");
|
||||
}
|
||||
|
||||
print_tex_footer();
|
||||
}
|
||||
75
linalg.c
75
linalg.c
@ -18,6 +18,7 @@ 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);
|
||||
@ -34,6 +35,7 @@ workspace_t *workspace_alloc(int n)
|
||||
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);
|
||||
@ -67,6 +69,33 @@ 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);
|
||||
@ -139,14 +168,15 @@ void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
|
||||
|
||||
void 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(g, ws->eval_complex, ws->evec_complex, 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 = 1;
|
||||
for(int i = 0; i < 3; i++)
|
||||
for(int i = 0; i < ws->n; i++)
|
||||
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
|
||||
real = 0;
|
||||
|
||||
@ -155,8 +185,45 @@ void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
|
||||
exit(1);
|
||||
}
|
||||
|
||||
for(int i = 0; i < 3; i++)
|
||||
for(int j = 0; j < 3; j++)
|
||||
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)));
|
||||
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
27
linalg.h
27
linalg.h
@ -3,6 +3,7 @@
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdarg.h>
|
||||
#include <gsl/gsl_math.h>
|
||||
#include <gsl/gsl_eigen.h>
|
||||
#include <gsl/gsl_blas.h>
|
||||
@ -13,16 +14,17 @@
|
||||
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
|
||||
|
||||
typedef struct _workspace {
|
||||
int n;
|
||||
gsl_eigen_nonsymmv_workspace *work_nonsymmv;
|
||||
gsl_vector *work_sv;
|
||||
gsl_vector_complex *eval_complex;
|
||||
gsl_matrix_complex *evec_complex;
|
||||
gsl_vector *eval_real;
|
||||
gsl_matrix *evec_real;
|
||||
gsl_matrix *tmp;
|
||||
gsl_permutation *permutation;
|
||||
gsl_matrix *stack[20];
|
||||
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_matrix *tmp;
|
||||
gsl_permutation *permutation;
|
||||
gsl_matrix *stack[20];
|
||||
} workspace_t;
|
||||
|
||||
workspace_t *workspace_alloc(int n);
|
||||
@ -30,6 +32,9 @@ void workspace_free(workspace_t *workspace);
|
||||
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);
|
||||
@ -37,5 +42,7 @@ void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||
double trace(gsl_matrix *g);
|
||||
double determinant(gsl_matrix *g, workspace_t *ws);
|
||||
void 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);
|
||||
|
||||
#endif
|
||||
|
||||
@ -28,6 +28,7 @@ int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int
|
||||
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);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user