#ifndef LINALG_H
#define LINALG_H

#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>

#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#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];
} workspace_t;

workspace_t *workspace_alloc(int n);
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 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);
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);

#endif