template <class T> class VVector;
template <class T> class VFlatVector;

/* ****************** T Vector ******************** */

/**
   The T_Vector specifies the type of the vector element
 */
template <typename T = double>
class T_Vector : public S_BaseVector<typename mat_traits<T>::TSCAL>
{
public:
  typedef T TELEM;
  typedef typename mat_traits<T>::TSCAL TSCAL;

  T_Vector ();
  virtual ~T_Vector();

  virtual FlatVector<T> FV () const = 0;
  virtual FlatVector<T> FV ()  = 0;


  FlatVector<T> FV (const BaseVector & v2) const 
  {
    return dynamic_cast<const T_Vector &> (v2).FV(); 
  }

  FlatVector<T> FV (BaseVector & v2)  
  {
    return dynamic_cast<T_Vector &> (v2).FV(); 
  }

  virtual int Size() const { return FV().Size(); }

  virtual TempVector Range (int begin, int end) 
  {
    return TempVector (new VFlatVector<T> (end-begin, &FV()[0]+begin));
  }
  
  virtual TempVector Range (int begin, int end) const
  {
    return TempVector (new VFlatVector<T> (end-begin, &FV()[0]+begin));
  }


  T_Vector & operator= (const T_Vector & v)
  {
    FV() = v.FV();
    return *this;
  }

  template<typename TB>
  T_Vector & operator= (const Expr<TB> & v)
  {
    FV() = v;
    return *this;
  }

  T_Vector & operator= (TSCAL s)
  {
    FV() = s;
    return *this;
  }

  T_Vector & operator= (const BaseVector & m)
  {
    FV() = FV(m);
    return *this;
  }

  template <typename T2>
  T_Vector & operator= (const VVecExpr<T2> & v)
  {
    v.AssignTo (1.0, *this);
    return *this;
  }



  T_Vector & operator+= (const T_Vector & v)
  {
    FV() += FV(v);
    return *this;
  }

  template <typename T2>
  T_Vector & operator+= (const VVecExpr<T2> & v)
  {
    v.AddTo (1.0, *this);
    return *this;
  }


  template<typename TB>
  T_Vector & operator+= (const Expr<TB> & v)
  {
    FV() += v;
    return *this;
  }

  T_Vector & operator+= (const BaseVector & v)
  {
    FV() += FV(v);
    return *this;
  }

  T_Vector & operator-= (const T_Vector & v)
  {
    FV() -= FV(v);
    return *this;
  }

  template <typename T2>
  T_Vector & operator-= (const VVecExpr<T2> & v)
  {
    v.AddTo (-1.0, *this);
    return *this;
  }


  template<typename TB>
  T_Vector & operator-= (const Expr<TB> & v)
  {
    FV() -= v;
    return *this;
  }

  T_Vector & operator-= (const BaseVector & v)
  {
    FV() -= FV(v);
    return *this;
  }







  T_Vector & operator*= (TSCAL scal)
  {
    FV() *= scal;
    return *this;
  }


  /*
  T_Vector & Assign (const BaseVector & v)
  {
    s = v.Size();
    data = static_cast<T*> (const_cast<void*> (v.Data()));
    return *this;
  }
  

  virtual const void * Data() const { return data; }
  virtual  void * Data() { return data; }
  virtual int Size() const { return s; }
  virtual int ElementSize() const { return sizeof(T); }
  */
  
  virtual BaseVector & SetScalar (TSCAL s)
  {
    FV() = s;
    return *this;
  }

  virtual TSCAL InnerProduct (const BaseVector & v2) const 
  {
    return ngbla::InnerProduct (FV(), FV(v2));
  }

  virtual double L2Norm () const
  {
    return ngbla::L2Norm (FV());
  }

  virtual BaseVector & Scale (TSCAL scal)
  {
    FV() *= scal;
    return *this;
  }

  virtual BaseVector & Set (TSCAL s, const BaseVector & v)
  {
    FV() = s * FV(v);
    return *this;
  }

  virtual BaseVector & Add (TSCAL s, const BaseVector & v)
  {
    FV() += s * FV(v);
    return *this;
  }

  virtual BaseVector * CreateVector () const
  {
    return new VVector<T> (Size());
  }


  void GetIndirect (const ARRAY<int> & ind, 
		    FlatVector<TSCAL> & v) const
  {
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	static_cast<TELEM*> (v.Data())[i] = FV()(ind[i]);
      else
	static_cast<TELEM*> (v.Data())[i] = 0;
  }

  virtual void SetIndirect (const ARRAY<int> & ind, 
			    const FlatVector<TSCAL> & v) 
  {
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	FV()(ind[i]) = static_cast<const TELEM*> (v.Data())[i];
      else	
	FV()(ind[i]) = 0;
  }

  virtual void AddIndirect (const ARRAY<int> & ind, 
			    const FlatVector<TSCAL> & v) 
  {
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	FV()(ind[i]) += static_cast<const TELEM*> (v.Data())[i];
  }





  virtual void SetRandom () 
  {
    FlatVector<T> fv = FV();
    for (int i = 0; i < Size(); i++)
      fv(i) = double (rand()) / RAND_MAX;
  }

  virtual ostream & Print (ostream & ost) const
  {
    // return (ost << RefMatExpr<FlatVector<T> > (*this) << endl);
    return (ost << FV() << endl);
  }

  virtual void Save(ostream & ost) const;


  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    mu.Append (new MemoryUsageStruct ("Vector: ", Size()*sizeof(T), 1));
  }
  
};