#ifndef __MYVECTOR__
#define __MYVECTOR__
#include <sys/types.h>
#include <string.h>
#include <assert.h>
#include <iostream>
#include <fstream>
using namespace std;
const double PI = 3.1415926535897932384e+00;
const double EEXP = 2.718281828459045235360287;
class Matrix;
class SparseMatrix;
extern "C" double dnrm2_(int *, double *, int *);
class Vector{
private:
  size_t size;
  double * data;
  int owner;
public:
  friend class Matrix;
  friend void eig(const Matrix& M, Matrix& V, Vector& er, Vector& ei);//er and ei 
  //are real and imaginary parts of eigenvalues; V as in LAPACK
  friend void eig(const Matrix& M, Vector& er, Vector& ei);//er and ei as above
  friend void svd(const Matrix& M, Vector& S);
  friend void svd(const Matrix& M, Vector& S, Matrix& U, Matrix& V);//S=singular
  //values U=left singular vectors V=right singular vectors
  friend void solve(const Matrix& M, const Vector& v, Vector& soln);
  friend void solve(const Matrix& M, Vector & v);//overwrites soln on v
  friend void mult(const Matrix& M1, const Matrix& M2, Matrix& prod);
  friend void mult(const Matrix& M1, Matrix& M2);//overwrites product on M2
  friend void mult(const Matrix& M, const Vector& v, Vector& prod);
  friend void mult(const Matrix& M, Vector& v);//overwrites product on v
  friend double det(const Matrix& M);//determinant
  friend void fft(const Vector& v, Vector& vf);
  friend void ifft(const Vector& vf, Vector& v);
  friend void fft_destroy_plan();//deallocate all fftw plans
  friend void diff(const Vector& v, Vector& vd);//vd = v'(t)
  friend void shift(const Vector& v, double h, Vector& w);//w(t)=v(t+h)
  friend void filter(const Vector& v, double alpha, Vector& w);
  friend void expm(const Matrix& M1, Matrix& M2);//M2 = expm(M1)
  friend void solve(SparseMatrix& M, const Vector& v, Vector& soln);
  friend void mult(SparseMatrix& M, const Vector& v, Vector& prod);
  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);
public:
  //empty constructor
  Vector(){
    size = 0;
    data = NULL;
    owner = 0;
  }
  
  //only constructor to allocate space for data
  Vector(size_t  n){
    size = n;
    data = new double[n];
    owner = 1;
  }
  
  
  //*this becomes shadow of Vector v (copy constructor)
  Vector(const Vector& v){
    size = v.size;
    data = v.data;
    owner = 0;
  } 

  

  //*this becomes shadow of v[i:k]
  Vector(const Vector& v, size_t i, size_t k){
    size_t n;
    assert(k>=i);
    assert(k < v.size);
    n = (k-i) + 1; 
    
    size = n;
    data = v.data;
    owner = 0;
  }

  //makes *this shadow v
  void shadow(const Vector& v){
    assert(!owner);
    size = v.size;
    data = v.data;
    owner = 0;
  }
	
  //makes *this shadow v(i:k)
  void shadow(const Vector& v,  size_t i,  size_t k){ 
    size_t n;
    assert(!owner);
    assert(k>=i);
    assert(k<v.size);
  
    n = (k-i)+1;
    size = n;
    data = v.data + i;
    owner = 0;
  }
  
  //returns a Vector shadowing (*this)(i:k)
  Vector getShadow(size_t i, size_t k){
      Vector v;
      assert(k>=i);
      assert(k<size);
      v.size = (k-i)+1;
      v.data = data + i;
      v.owner = 0;
      return v;
  }
      

  
  //destructor, finds out if *this is a shadow or owns space
  ~Vector(){
    if(owner!=0)
      delete[] data;
  }
  

  //returns (*this)(i)
  double& operator()(size_t i){
    assert(i < size);
    return(data[i]);
  }
  
  //returns (*this)(i)
  const double& operator()(size_t i) const{
    assert(i<size);
    return(data[i]);
  }
		
  //returns (*this)(i)
  double& operator[](size_t i){
    assert(i < size);
    return(data[i]);
  }
  
  //retuns (*this)(i)
  const double& operator[](size_t i) const{
    assert(i < size);
    return(data[i]);
  }

  //retuns length of (*this)
  size_t getSize() const{
    return size;
  }

  //(*this)(i) = v(i) for all i
  Vector& operator=(const Vector& v){
    assert(size==v.size);
    assert(((data<v.data)&&((data+size-1)<v.data))||
	   ((data>v.data)&& (data>(v.data+size-1)))); //no overlap
    memcpy((char *)data, (char *)v.data, size*sizeof(double));
    return(*this);
  }


  //*this = *this + v
  void add(const Vector& v){
    assert(size == v.size);
    for(size_t i=0; i < size; i++)
      data[i] += (v.data)[i];
  }

  //*this = *this - v
  void sub(const Vector& v){
    assert(size == v.size);
    for(size_t i=0; i < size; i++)
      data[i] -= (v.data)[i];
  }

  //*this = *this .* v
  void mul(const Vector& v){
    assert(size == v.size);
    for(size_t i=0; i < size; i++)
      data[i] *= (v.data)[i];
  }

  //*this = *this ./ v
  void div(const Vector& v){
    assert(size == v.size);
    for(size_t i=0; i < size; i++)
      data[i] /= (v.data)[i];
  }

  //*this = x * (*this)
  void scale(const double x){
    for(size_t i=0; i < size; i++)
      data[i] *= x;
  }

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

  
  
 
  //returns 2-norm of *this
  double norm() const{
    int dim = size;
    double * vec = data;
    int stride = 1;
    return dnrm2_(&dim, vec, &stride); 
  }

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

  //input vector from file fname
  void input(const char* fname){
    ifstream ifile(fname);
    assert(!ifile.fail());
    size_t i;
    for(i=0; i < size; i++){
      assert(!ifile.eof());
      ifile>>(*this)(i);
    }
    ifile.close();
  }

  //use with care!
  // double * getRawData() const{
  //  return data;
  //}
};
	

#endif