#ifndef __MYMATRIX__
#define __MYMATRIX__
#include <sys/types.h>
#include <assert.h>
#include <stdlib.h>
#include <math.h>
#include "Vector.hh"
using namespace std;
extern void GetSpace(size_t n);
extern void FreeSpace();


class Matrix{
private:
  size_t size1;
  size_t size2;
  size_t lda;
  double * data;
  int owner;

  
public:
  friend class Vector;
  friend void eig(const Matrix& M, Matrix& V, Vector& er, Vector& ei);
  friend void eig(const Matrix& M, Vector& er, Vector& ei);
  friend void svd(const Matrix& M, Vector& S);
  friend void svd(const Matrix& M, Vector& S, Matrix& U, Matrix& V);
  friend void solve(const Matrix& M, const Vector& v, Vector& soln);
  friend void solve(const Matrix& M, const Vector& v);
  friend void mult(const Matrix& M1, const Matrix& M2, Matrix& prod);
  friend void mult(const Matrix& M1, Matrix& M2);
  friend void mult(const Matrix& M, const Vector& v, Vector& prod);
  friend void mult(const Matrix& M, Vector& v);
  friend void transpose(const Matrix& M1, Matrix& M2);
  friend void invert(const Matrix& M1, Matrix& M2);
  friend double det(const Matrix& M);
  friend void fft(const Vector& v, Vector& vf);
  friend void fft_destroy_plan();//deallocate all fftw plans
  friend void ifft(const Vector& vf, Vector& v);
  friend void diff(const Vector& v, Vector& vd);
  friend void diff(const Matrix& M, Matrix & N);//differentiates row-by-row
  friend void shift(const Vector& v, double h, Vector& w);
  friend void shift(const Matrix& M, double h, Matrix& N); //shifts row-wise
  friend void filter(const Vector& v, double alpha, Vector& w);
  friend void expm(const Matrix& M1, Matrix& M2);
  friend void spline(const Vector& t, const Vector& v,
		     double leftD, double rightD,
		     const Vector& tt, Vector& vv);
  friend void spline(const Vector& t, const Matrix& M,
		     const Vector& leftD, const Vector& rightD,
		     const Vector& tt, Matrix& MM);
  friend void cspline(const Vector& t, const Vector& v, 
		      const Vector& tt, Vector& vv);
  friend void cspline(const Vector& t, const Matrix& M, 
		      const Vector& tt, Matrix& MM);
  friend void matrixsmooth(const Matrix& X, double alpha, Matrix& XX);
public:
  //empty constructor
  Matrix(){
    size1 = 0;
    size2 = 0;
    lda = 0;
    data = NULL;
    owner = 0;
  }
  
  //only constructor to actually allocate space for data
  Matrix(size_t m, size_t n){
    size1 = m;
    size2 = n;
    lda = m;
    data = new double[m*n];
    owner = 1;
  }

    
  //copy constructor, makes *this a shadow of M
  Matrix(const Matrix& M){
    size1 = M.size1;
    size2 = M.size2;
    lda = M.lda;
    data = M.data;
    owner = 0;
  }

  //makes *this a shadow of M(i1:i2, j1:j2)
  Matrix(const Matrix& M, size_t i1, size_t i2, size_t j1, size_t j2){

    assert((i1<=i2)&&(j1<=j2)&&(i2<M.size1)&&(j2<M.size2));
    size1 = (i2-i1)+1;
    size2 = (j2-j1)+1;
    lda = M.lda;
    data = M.data + j1*lda+i1;
    owner = 0;
  }

  //makes *this a shadow of M
  void shadow(const Matrix & M){
    assert(!owner);
    size1 = M.size1;
    size2 = M.size2;
    lda = M.lda;
    data = M.data;
    owner = 0;
  }

  //makes *this a shadow of M(i1:i2,j1:j2)
  void shadow(const Matrix & M, size_t i1, size_t i2, size_t j1, size_t j2){
    assert(!owner);
    assert((i1<=i2)&&(j1<=j2)&&(i2<M.size1)&&(j2<M.size2));

    size1 = (i2-i1)+1;
    size2 = (j2-j1)+1;
    lda = M.lda;
    data = M.data + j1*lda+i1;
    owner = 0;
  }
  
 
  //makes (*this) a m by n matrix shadowing v(i1:i2)
  void shadow(const Vector& v, size_t i1, size_t i2, size_t m, size_t n){
    assert((i2-i1+1)==m*n);
    assert(i1 <= i2);
    assert(i2 < v.getSize());
    assert(!owner);
    #ifndef DEBUG
    i2 = 0; //takes out compile time warning
    #endif
    size1 = m;
    size2 = n;
    lda = m;
    data = v.data + i1;
    owner = 0;
  }

  //makes *this a m by n matrix shadowing v(i1:i1+m*n-1)
  void shadow(const Vector& v, size_t i1, size_t m, size_t n){
    this->shadow(v,i1,i1+m*n-1,m,n);
  }

  //makes *this a m by n matrix shadowing v(0:m*n-1)
  void shadow(const Vector& v, size_t m, size_t n){
    this->shadow(v,0,m*n-1,m,n);
  }

  //returns shadow of (*this)(i1:i2,j1:j2)
  Matrix getShadow(size_t i1, size_t i2, size_t j1, size_t j2)const{
      Matrix temp;
      assert((i1<=i2)&&(j1<=j2)&&(i2<size1)&&(j2<size2));
      temp.size1 = (i2-i1)+1;
      temp.size2 = (j2-j1)+1;
      temp.lda = lda;
      temp.data = data + j1*lda + i1;
      temp.owner = 0;
      return temp;
  }
  
  //destructor, checks if *this is a shadow or owns space
  ~Matrix(){
    if(owner!=0)
      delete[] data;
  }

  //returns number of rows of *this
  size_t getm() const{
    return size1;
  }
  
  //returns number of columns of *this
  size_t getn() const{
    return size2;
  }

  //returns reference to (i,j)th entry of *this
  double& operator()(size_t i, size_t j){
    assert(i<size1);
    assert(j<size2);
    return data[j*lda+i];
  }
	
  //returns reference ot (i,j)th entry of *this
  const double& operator()(size_t i, size_t j)const{
    assert(i<size1);
    assert(j<size2);
    return data[j*lda+i];

  }

  //returns reference to (i)th entry for rows and columns
  double& operator()(size_t i){
    assert((size1==1)||(size2==1));
    if(size1==1)
      return (*this)(0,i);
    else if(size2==1)
      return (*this)(i,0);
  }
      

  //*this(i,j) is set to M(i,j) for all (i,j)
  Matrix& operator=(const Matrix& M){
    assert((size1==M.size1)&&(size2==M.size2));
    for(size_t j=0; j < size2; j++)
      memcpy((char *)(data+lda*j), (char *)(M.data+M.lda*j), 
	     size1*sizeof(double));
    return(*this);
  }

  
  // returns a vector that shadows column j
  Vector getColumn(size_t j) const{
    Vector v;
    assert(j<size2);
    v.size = size1;
    v.data = data + j*lda;
    v.owner = 0;
    return v;
  }

  //transpose square matrix
  void transpose(){
    assert(getm()==getn());
    for(size_t i=0; i < getm(); i++)
      for(size_t j = i+1; j < getn(); j++){
	double temp = (*this)(i,j);
	(*this)(i,j) = (*this)(j,i);
	(*this)(j,i) = temp;
      }
  }
	
     
  //*this = *this+M
  void add(const Matrix& M){
    assert((size1 == M.size1)&&(size2==M.size2));
    for(size_t j=0; j < size2; j++)
      for(size_t i=0; i < size1; i++)
	data[j*lda+i] += (M.data)[j*M.lda+i];
  }
  
  //*this = *this-M
  void sub(const Matrix& M){
    assert((size1 == M.size1)&&(size2==M.size2));
    for(size_t j=0; j < size2; j++)
      for(size_t i=0; i < size1; i++)
	data[j*lda+i] -= (M.data)[j*M.lda+i];
  }

  //*this = (*this).*M
  void mul(const Matrix& M){
    assert((size1 == M.size1)&&(size2==M.size2));
    for(size_t j=0; j < size2; j++)
      for(size_t i=0; i < size1; i++)
	data[j*lda+i] *= (M.data)[j*M.lda+i];
  }
  
  //*this = (*this)./M
  void div(const Matrix& M){
    assert((size1 == M.size1)&&(size2==M.size2));
    for(size_t j=0; j < size2; j++)
      for(size_t i=0; i < size1; i++)
	data[j*lda+i] /= (M.data)[j*M.lda+i]; 
  }
  
  //*this = x * (*this)
  void scale(const double x){
    for(size_t j=0; j < size2; j++)
      for(size_t i=0; i < size1; i++)
	data[j*lda+i] *= x; 
  }

  //*this = (*this) + x
  void add_constant(const double x){
    for(size_t j=0; j < size2; j++)
      for(size_t i=0; i < size1; i++)
	data[j*lda+i] += x; 
  }


  
  //return 2-norm of *this
  double norm(){
      size_t n;
      n = (getm()<=getn())? getm(): getn();
      Vector S(n);
      svd(*this, S);
      return S[0];
  }
      

  //return frobenius norm of *this
  double fronorm() const{
      double ans = 0;
      for(size_t i = 0; i < getn(); i++){
	  double temp = (this->getColumn(i)).norm();
	  ans += temp*temp;
      }

      return sqrt(ans);
  }

  //output matrix to file fname
  void output(const char *fname)const{
    ofstream ofile(fname);
    ofile.setf(ios::scientific, ios::floatfield);
    ofile.precision(16);
    size_t i,j;
    for(i=0; i < size1; i++){
      for(j=0; j < size2; j++)
	ofile<<(*this)(i,j)<<" ";
      ofile<<endl;
    }
    ofile.close();
  }

  //input matrix from file fname
  void input(const char *fname){
    ifstream ifile(fname);
    assert(!ifile.fail());
    size_t i,j;
    for(i=0; i < size1; i++)
      for(j=0; j < size2; j++){
	assert(!ifile.eof());
	ifile>>(*this)(i,j);
      }
    ifile.close();
  }
      
  //use with care!
  // double * getRawData()const{
  // return data;
  //}
  
};






#endif