#ifndef EMATRIX_H
#define EMATRIX_H

#include "eutils.h"

#include <iostream>
#include "efile.h"

#ifdef EUTILS_HAVE_LIBGSL
  #include <gsl/gsl_linalg.h>
#endif

using namespace std;

class ematrix;

class evector
{
 public:
  int w;
  double *vector;
#ifdef EUTILS_HAVE_LIBGSL
  gsl_vector_view gsl_vecview;
#endif

  evector();
  evector(const evector& v);
  evector(int w,double v);
  evector(int w,double *values=0x00);
  ~evector();

  void init(int w,double v);
  void init(int w);

  inline void create(int w) { init(w); }
  void clear();

  double length() const;

  double& operator()(int x);
  const double& operator()(int x) const;
  evector& operator=(const evector& v);
  evector& operator+=(const evector& v);
  evector& operator-=(const evector& v);

  inline double& operator[](int x) { return(vector[x]); }
  inline double& operator[](int x) const { return(vector[x]); }

  inline int size() const { return(w); }

  evector operator*(const ematrix& m) const;
  double  operator*(const evector& v) const;
  evector operator*(double f) const;
};

ostream& operator<<(ostream& stream,const evector& v);

class ematrix
{

 public:
  int w;
  int h;
  double *matrix;
#ifdef EUTILS_HAVE_LIBGSL
  gsl_matrix_view gsl_matview;
#endif

  ematrix();
  ematrix(const ematrix& m);
  ematrix(int w,int h,double *values=0x00);
  ~ematrix();

  void create(const ematrix& m);
  void create(int w,int h);
  void clear();

  void load(efile file);
  void save(efile file);

  double& operator()(int x,int y);
  const double& operator()(int x,int y) const;

  evector row(int i) const;
  evector col(int j) const;

  void swap(int i,int i2);
  void swapcols(int j,int j2);
  void mulrow(int i,double f);
  void addmulrow(int i,double f,int i2);

  ematrix& operator=(const ematrix& m);
  void copytranspose(const ematrix& mat);

  evector operator*(const evector& v) const;
  ematrix operator*(const ematrix& m) const;
};

void nullspace(ematrix& A,ematrix& n);
void nullspacer(ematrix& A,ematrix& n);
int rref(ematrix& A);
void svd(ematrix& A,ematrix& V,evector& S);

evector eigenvalues(ematrix& A);

void pcoa(ematrix& A,evector& eval,ematrix& evec);


ematrix divmat(const ematrix& farr,const ematrix& farr2);
ematrix submat(const ematrix& farr,const ematrix& farr2);
ematrix summat(const ematrix& farr,const ematrix& farr2);
ematrix mulmat(const ematrix& farr,const ematrix& farr2);
ematrix tanmat(const ematrix& farr);
ematrix logmat(const ematrix& farr);
ematrix maxmat(const ematrix& farr,double mval);
ematrix minmat(const ematrix& farr,double mval);
ematrix normmat(const ematrix& farr);
ematrix mulmat(const ematrix& farr,double mval);
ematrix invmat(const ematrix& farr);
ematrix summat(const ematrix& farr,double mval);
ematrix submat(const ematrix& farr,double mval);
ematrix divmat(const ematrix& farr,double mval);

ostream& operator<<(ostream& stream,const ematrix& m);

#endif