#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 600000
#define MAX_TEMP_VECTORS 100

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);
int real_eigenvalues(gsl_matrix *g, gsl_vector *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);
int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, 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