template <int S, class T> class Vec;
template <class T> class SysVector;
template <class T> class FlatVector;
template <class T> class Vector;


extern void CheckVecRange(int s, int i);
extern void CheckVecRange(int s, int i, int j);


/*
   Static base class for \Ref{FlatVector}.
   Provides size, range check methods.
 */
/*
class BaseFlatVector
{
protected:
  unsigned int s;
  BaseFlatVector (int as) : s(as) { ; }

  void CheckRange(int i) const;
  void CheckRange(int i, int j) const;
};
*/



/**
   A simple vector.
   Has size and generic data-pointer. 
   No memory allocation/deallocation. User must provide memory.
 */
template <typename T = double>
class FlatVector : public MatExpr<FlatVector<T> > 
{
protected:
  /// vector size
  int s;
  /// the data
  T * data;
public:
  typedef T TELEM;
  typedef typename mat_traits<T>::TSCAL TSCAL;

  /// 
  FlatVector () { ; }
  /// set size and mem
  FlatVector (unsigned int as, T * adata) : s(as), data(adata) { ; }

  /// set size and mem
  FlatVector (unsigned int as, void * adata) : s(as), data(static_cast<TELEM*> (adata)) { ; }

  /// useful to put FlatVector over Vec
  template <typename T2>
  explicit FlatVector (const MatExpr<T2> & v)
    : s(v.Height()), data(const_cast<T*>(&v.Spec()(0)))
  { ; }

  template <int S>
  FlatVector (Vec<S,TSCAL> & v)
    : s(v.Size()), 
      data(const_cast<T*>(&v(0)))
  { ; }

  
  FlatVector (int as, LocalHeap & lh) 
    : s(as), data((T*)lh.Alloc(s*sizeof(T))) { ; }

  /*
  template<class T2>
  FlatVector (ExprAndMemory<T2> e)
    : s(e.Exp().Height()), data(e.Memory())
  {
    for (int i = 0; i < s; i++)
      data[i] = e.Exp()(i);
  }
  */

  FlatVector (const SysVector<TSCAL> & sv)
    : s(sv.Size()*sv.BlockSize() / mat_traits<T>::VDIM), 
      data (sv(0))
  {
    ;
  }


  void AssignMemory (int as, LocalHeap & lh) 
  {
    s = as;
    data = (T*)lh.Alloc(s*sizeof(T));
  }

  void AssignMemory (int as, T * mem) 
  {
    s = as;
    data = mem;
  }


  
  FlatVector & operator= (const FlatVector & v)
  {
    for (int i = 0; i < s; i++)
      data[i] = v(i);
    return *this;
  }

  template<typename TB>
  FlatVector & operator= (const Expr<TB> & v)
  {
    return MatExpr<FlatVector>::operator= (v);
  }

  FlatVector & operator= (TSCAL scal)
  {
    for (int i = 0; i < s; i++)
      data[i] = scal; 
    return *this;
  }



  TELEM & operator() (int i) 
  {
#ifdef CHECK_RANGE
    CheckVecRange(s,i);
#endif
    return data[i]; 
  }

  const TELEM & operator() (int i) const
  {
#ifdef CHECK_RANGE
    CheckVecRange(s,i);
#endif
    return data[i]; 
  }

  const TELEM & operator() (int i, int j) const
  {
#ifdef CHECK_RANGE
    CheckVecRange(s,i);
#endif
    return data[i]; 
  }

  TELEM & operator[] (int i) 
  {
#ifdef CHECK_RANGE
    CheckVecRange(s,i);
#endif
    return data[i]; 
  }

  const TELEM & operator[] (int i) const
  {
#ifdef CHECK_RANGE
    CheckVecRange(s,i);
#endif
    return data[i]; 
  }


  FlatVector<T> Range(int first, int last)
  { return FlatVector<T> (last-first+1, data+first); }
  const FlatVector<T> Range(int first, int last) const
  { return FlatVector<T> (last-first+1, data+first); }

  int Size () const { return s; }
  int Height () const { return s; }
  int Width () const { return 1; }
  const void * Data () const { return static_cast<const void*>(data); }
  void * Data () { return static_cast<void*>(data); }
};



/**
   A Vector class with memory allocation/deallocation
 */
template <typename T = double>
class Vector : public FlatVector<T>
{
public:
  typedef typename mat_traits<T>::TSCAL TSCAL;

  explicit Vector (int as) : FlatVector<T> (as, new T[as]) { ; }
  ~Vector() { delete [] this->data; }

  Vector & operator= (TSCAL scal)
  {
    FlatVector<T>::operator= (scal);
    return *this;
  }

  void SetSize(int as)
  {
    if (this->s == as) return;
    delete [] this->data;
    this->s = as;
    this->data = new T[this->s];
  }

  template<typename TB>
  Vector & operator= (const Expr<TB> & v)
  {
    MatExpr<FlatVector<T> >::operator= (v);
    return *this;
  }
};











// A system vector (non completely functional)
template <typename T = double>
class FlatSysVector : public MatExpr<FlatSysVector<T> >
{
protected:
  int s;
  int blocksize;
  T * data;
public:
  typedef FlatVector<T> TELEM;
  typedef typename mat_traits<T>::TSCAL TSCAL;
  
  FlatSysVector (int as, int bs, T * adata) 
    : s(as), blocksize(bs), data(adata) { ; }
  
  FlatSysVector (int as, int bs, LocalHeap & lh) 
    : s(as), blocksize(bs), data (new (lh) T[as*bs])
  { ; }
  
  FlatSysVector (const FlatSysVector<TSCAL> & sv)
    : s(sv.Size()),
      blocksize(sv.BlockSize()),
      data (sv(0))
  {
    ;
  }

  
  FlatSysVector & operator= (const FlatSysVector & v)
  {
    for (int i = 0; i < this->s * this->bs; i++)
      data[i] = v.data[i];
    return *this;
  }
  
  template<typename TB>
  FlatSysVector & operator= (const Expr<TB> & v)
  {
    return MatExpr<FlatSysVector>::operator= (v);
  }
  
  FlatSysVector & operator= (TSCAL s)
  {
    for (int i = 0; i < this->s*this->bs; i++)
      data[i] = s;
    return *this;
  }

  

  TELEM operator() (int i) 
  {
    return FlatVector<T> (blocksize, &data[i*blocksize]); 
  }

  const TELEM operator() (int i) const
  {
    return FlatVector<T> (blocksize, &data[i*blocksize]); 
  }

  const TELEM operator() (int i, int j) const
  {
    return FlatVector<T> (blocksize, &data[i*blocksize]); 
  }


  FlatSysVector<T> Range(int first, int last)
  { return FlatSysVector<T> (last-first+1, blocksize, data+(first*blocksize)); }
  const FlatSysVector<T> Range(int first, int last) const
  { return FlatSysVector<T> (last-first+1, blocksize, data+(first*blocksize)); }

  int Size () const { return s; }
  int Height () const { return s; }
  int Width () const { return 1; }
};










/**
   A vector of fixed size.
   Useful as entry type in system vectors.
 */
template <int S, typename T = double>
class Vec : public MatExpr<Vec<S,T> > // , protected BaseVec
{
  /// the values
  T data[S];
public:
  /// type of the elements
  typedef T TELEM;
  /// is the element double or complex ?
  typedef typename mat_traits<T>::TSCAL TSCAL;
  /// a vec is a S times 1 matrix, the according colume vector
  typedef Vec<S, typename mat_traits<T>::TV_COL> TV_COL;
  /// a vec is a S times 1 matrix, the according row vector
  typedef Vec<1, typename mat_traits<T>::TV_ROW> TV_ROW;

  enum { SIZE = S };
  /// height of matrix
  enum { HEIGHT = S };
  /// with of matrix
  enum { WIDTH  = 1 };

  /// constructor, no initialization
  Vec () { ; }
  /// copy vector
  Vec (const Vec & v) : MatExpr<Vec> ()
  {
    for (int i = 0; i < S; i++)
      data[i] = v.data[i];
  }

  /// initialize with values
  Vec (const TSCAL & scal)
  {
    for (int i = 0; i < S; i++)
      data[i] = scal;
  }

  /// initialize with expression
  template<typename TB>
  Vec (const Expr<TB> & v)
  {
    (*this) = v;
  }

  /// initialize zeroth and first elements
  Vec (const T & s1, const T & s2)
  {
    data[0] = s1;
    data[1] = s2;
  }

  /// initialize zeroth, first, and second elements
  Vec (const T & s1, const T & s2, const T & s3)
  {
    data[0] = s1;
    data[1] = s2;
    data[2] = s3;
  }


  /// copy vector
  Vec & operator= (const Vec & v)
  {
    for (int i = 0; i < S; i++)
      data[i] = v.data[i];
    return *this;
  }

  /// assign scalar value
  Vec & operator= (TSCAL scal)
  {
    for (int i = 0; i < S; i++)
      data[i] = scal;
    return *this;
  }

  /// assign expression
  template<typename TB>
  Vec & operator= (const Expr<TB> & v)
  {
    MatExpr<Vec>::operator= (v);
    return *this;
  }


  /// access vector
  TELEM & operator() (int i) 
  {
#ifdef CHECK_RANGE
    CheckVecRange(S,i);
#endif
    return data[i]; 
  }

  /// access vector
  const TELEM & operator() (int i) const 
  {
#ifdef CHECK_RANGE
    CheckVecRange(S,i);
#endif
    return data[i]; 
  }

  /// access vector
  const TELEM & operator() (int i, int j) const 
  {
#ifdef CHECK_RANGE
    CheckVecRange(S,i);
#endif
    return data[i]; 
  }

  /// vector size
  int Size () const { return S; }
  /// corresponding matrix height
  int Height () const { return S; }
  /// corresponding matrix with
  int Width () const { return 1; }
};


/// cross product
inline Vec<3> Cross (const Vec<3> & a, const Vec<3> & b)
{
  Vec<3> prod;
  prod(0) = a(1) * b(2) - a(2) * b(1);
  prod(1) = a(2) * b(0) - a(0) * b(2);
  prod(2) = a(0) * b(1) - a(1) * b(0);
  return prod;
}

/// output vector
template<int S, typename T>
inline ostream & operator<< (ostream & ost, const Vec<S,T> & v)
{
  for (int i = 0; i < S; i++)
    ost << " " << setw(7) << v(i);
  return ost;
}