class BaseVector;
typedef auto_ptr<BaseVector> TempVector;
template <class SCAL> class S_BaseVector;


/**
   Base class to linalg expression templates
 */
template <class T>
class VVecExpr 
{
  const T data;
public:
  ///
  VVecExpr (const T & d) : data(d) { ; }

  /// assign s * vector-expression data to v
  template <class TS>
  void AssignTo (TS s , BaseVector & v) const { data.AssignTo(s, v); }

  /// add s * vector-expression data to v
  template <class TS>
  void AddTo (TS s, BaseVector & v) const { data.AddTo(s, v); }
};


/**
   Base vector for linalg
 */
class BaseVector
{
protected:
  /// size of vector
  int size;
  /// number of doubles per entry
  int entrysize;

public:
  ///
  BaseVector () throw () { ; }
  ///
  virtual ~BaseVector () throw () { ; }

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

  ///
  BaseVector & operator= (const BaseVector & v);
  ///
  BaseVector & operator= (double s);
  ///
  BaseVector & operator= (Complex s);

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

  ///
  BaseVector & operator+= (const BaseVector & v)
  {
    Add (1.0, v);
    return *this;
  }

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

  ///
  BaseVector & operator-= (const BaseVector & v)
  {
    Add (-1.0, v);
    return *this;
  }

  ///
  BaseVector & operator*= (double s)
  {
    return Scale (s);
  }

  ///
  BaseVector & operator*= (Complex s)
  {
    return Scale (s);
  }

  ///
  BaseVector & operator/= (double s)
  {
    if (s == 0)
      throw Exception ("BaseVector::operator/=: division by zero");
    return Scale (1/s);
  }

  ///
  BaseVector & operator/= (Complex s)
  {
    if (s == 0.0)
      throw Exception ("BaseVector::operator/=: division by zero");
    return Scale (1.0/s);
  }

  template <class SCAL>
  S_BaseVector<SCAL> & Spec()
  {
    return dynamic_cast<S_BaseVector<SCAL>&> (*this);
  }
  
  template <class SCAL>
  const S_BaseVector<SCAL> & Spec() const
  {
    return dynamic_cast<const S_BaseVector<SCAL>&> (*this);
  }

  int Size() const throw () { return size; }
  int EntrySize() const throw () { return entrysize; }
  virtual void * Memory () const throw () = 0;

  virtual FlatVector<double> FVDouble () const throw() = 0;
  virtual FlatVector<Complex> FVComplex () const throw() = 0;

  template <class TSCAL>
  TSCAL InnerProduct (const BaseVector & v2) const 
  {
    return dynamic_cast<const S_BaseVector<TSCAL>&> (*this) . 
      InnerProduct (v2);
  }

  double L2Norm () const
  {
    return ngbla::L2Norm (FVDouble());
  }

  BaseVector & Scale (double scal);
  BaseVector & Scale (Complex scal);

  BaseVector & SetScalar (double scal);
  BaseVector & SetScalar (Complex scal);

  BaseVector & Set (double scal, const BaseVector & v);
  BaseVector & Set (Complex scal, const BaseVector & v);

  BaseVector & Add (double scal, const BaseVector & v);
  BaseVector & Add (Complex scal, const BaseVector & v);

  virtual ostream & Print (ostream & ost) const;
  virtual void Save(ostream & ost) const;
  virtual void Load(istream & ist);

  virtual void MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const;
  virtual BaseVector * CreateVector () const;

  virtual void SetRandom ();

  virtual TempVector Range (int begin, int end);
  virtual TempVector Range (int begin, int end) const;


  void GetIndirect (const ARRAY<int> & ind, 
		    FlatVector<double> & v) const;
  void GetIndirect (const ARRAY<int> & ind, 
		    FlatVector<Complex> & v) const;
  void SetIndirect (const ARRAY<int> & ind, 
		    const FlatVector<double> & v);
  void SetIndirect (const ARRAY<int> & ind, 
		    const FlatVector<Complex> & v);
  void AddIndirect (const ARRAY<int> & ind, 
		    const FlatVector<double> & v);
  void AddIndirect (const ARRAY<int> & ind, 
		    const FlatVector<Complex> & v);
};



/**
   Decision between double or Complex
 */
template <class SCAL>
class S_BaseVector : public BaseVector
{
public:
  S_BaseVector () throw () { ; }
  virtual ~S_BaseVector() throw() { ; }

  S_BaseVector & operator= (double s);

  SCAL InnerProduct (const BaseVector & v2) const
  {
    return ngbla::InnerProduct (FVScal(), 
				dynamic_cast<const S_BaseVector&>(v2).FVScal());
  }

  virtual FlatVector<double> FVDouble () const throw();
  virtual FlatVector<Complex> FVComplex () const throw();
  virtual FlatVector<SCAL> FVScal () const throw() 
  {
    return FlatVector<SCAL> (size * entrysize, Memory());
  }
};




template <>
class S_BaseVector<Complex> : public BaseVector
{
public:
  S_BaseVector () throw() { ; }
  ~S_BaseVector () throw() { ; }

  Complex InnerProduct (const BaseVector & v2) const
  {
    return ngbla::InnerProduct (FVScal(), 
				dynamic_cast<const S_BaseVector&>(v2).FVScal());
  }

  virtual FlatVector<double> FVDouble () const throw();
  virtual FlatVector<Complex> FVComplex () const throw();
  virtual FlatVector<Complex> FVScal () const throw() 
  {
    return FlatVector<Complex> (size * entrysize/2, Memory());
  }
};









/* ********************* Expression templates ***************** */



template <> class VVecExpr<BaseVector>
{
  const BaseVector & v;
public:
  VVecExpr (const BaseVector & av) : v(av) { ; }

  template <class TS>
  void AssignTo (TS s, BaseVector & v2) const { v2.Set (s, v); }
  template <class TS>
  void AddTo (TS s, BaseVector & v2) const { v2.Add (s,  v); }
};



/* ***************************** VSumExpr ************************** */

///
template <class TA, class TB>
class VSumExpr
{
  const TA a;
  const TB b;
  
public:
  VSumExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; }

  template <class TS>
  void AssignTo (TS s, BaseVector & v) const
  { 
    a.AssignTo (s, v);
    b.AddTo (s, v);
  }
  template <class TS>
  void AddTo (TS s, BaseVector & v) const
  { 
    a.AddTo (s, v);
    b.AddTo (s, v);
  }
};



inline VVecExpr<VSumExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > >
operator+ (const BaseVector & a, const BaseVector & b)
{
  typedef VSumExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > TRES;
  return TRES (a, b);
}

template <class TA>
inline VVecExpr<VSumExpr<VVecExpr<TA>, VVecExpr<BaseVector> > >
operator+ (const VVecExpr<TA> & a, const BaseVector & b)
{
  typedef VSumExpr<VVecExpr<TA>, VVecExpr<BaseVector> > TRES;
  return TRES (a, b);
}

template <class TB>
inline VVecExpr<VSumExpr<VVecExpr<BaseVector>, VVecExpr<TB> > >
operator+ (const BaseVector & a, const VVecExpr<TB> & b)
{
  typedef VSumExpr<VVecExpr<BaseVector>, VVecExpr<TB> > TRES;
  return TRES (a, b);
}

template <class TA, class TB>
inline VVecExpr<VSumExpr<VVecExpr<TA>, VVecExpr<TB> > >
operator+ (const VVecExpr<TA> & a, const VVecExpr<TB> & b)
{
  typedef VSumExpr<VVecExpr<TA>, VVecExpr<TB> > TRES;
  return TRES (a, b);
}








/* ***************************** VSubExpr ************************** */

///
template <class TA, class TB>
class VSubExpr
{
  const TA a;
  const TB b;
  
public:
  VSubExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; }


  template <class TS>
  void AssignTo (TS s, BaseVector & v) const
  { 
    a.AssignTo (s, v);
    b.AddTo (-s, v);
  }
  template <class TS>
  void AddTo (TS s, BaseVector & v) const
  { 
    a.AddTo (s, v);
    b.AddTo (-s, v);
  }
};



inline VVecExpr<VSubExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > >
operator- (const BaseVector & a, const BaseVector & b)
{
  typedef VSubExpr<VVecExpr<BaseVector>, VVecExpr<BaseVector> > TRES;
  return TRES (a, b);
}

template <class TA>
inline VVecExpr<VSubExpr<VVecExpr<TA>, VVecExpr<BaseVector> > >
operator- (const VVecExpr<TA> & a, const BaseVector & b)
{
  typedef VSubExpr<VVecExpr<TA>, VVecExpr<BaseVector> > TRES;
  return TRES (a, b);
}

template <class TB>
inline VVecExpr<VSubExpr<VVecExpr<BaseVector>, VVecExpr<TB> > >
operator- (const BaseVector & a, const VVecExpr<TB> & b)
{
  typedef VSubExpr<VVecExpr<BaseVector>, VVecExpr<TB> > TRES;
  return TRES (a, b);
}

template <class TA, class TB>
inline VVecExpr<VSubExpr<VVecExpr<TA>, VVecExpr<TB> > >
operator- (const VVecExpr<TA> & a, const VVecExpr<TB> & b)
{
  typedef VSubExpr<VVecExpr<TA>, VVecExpr<TB> > TRES;
  return TRES (a, b);
}





/* ************************* Scal * Vec ******************** */


///
template <class TA, class TSCAL>
class VScaleExpr
{
  const TA a;
  const TSCAL scal;
  
public:
  VScaleExpr (const TA & aa, const TSCAL & as) : a(aa), scal(as) { ; }


  template <class TS>
  void AssignTo (TS s, BaseVector & v) const
  { 
    a.AssignTo (scal * s, v);
  }
  template <class TS>
  void AddTo (TS s, BaseVector & v) const
  { 
    a.AddTo (scal * s, v);
  }
};



inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, double> >
operator* (const BaseVector & a, const double & b)
{
  typedef VScaleExpr<VVecExpr<BaseVector>, double> TRES;
  return TRES (a, b);
}

template <class TA>
inline VVecExpr<VScaleExpr<VVecExpr<TA>, double> >
operator* (const VVecExpr<TA> & a, const double & b)
{
  typedef VScaleExpr<VVecExpr<TA>, double> TRES;
  return TRES (a, b);
}



inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, Complex> >
operator* (const BaseVector & a, const Complex & b)
{
  typedef VScaleExpr<VVecExpr<BaseVector>, Complex> TRES;
  return TRES (a, b);
}

template <class TA>
inline VVecExpr<VScaleExpr<VVecExpr<TA>, Complex> >
operator* (const VVecExpr<TA> & a, const Complex & b)
{
  typedef VScaleExpr<VVecExpr<TA>, Complex> TRES;
  return TRES (a, b);
}





inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, double> >
operator* (const double & b, const BaseVector & a)
{
  typedef VScaleExpr<VVecExpr<BaseVector>, double> TRES;
  return TRES (a, b);
}

template <class TA>
inline VVecExpr<VScaleExpr<VVecExpr<TA>, double> >
operator* (const double & b, const VVecExpr<TA> & a)
{
  typedef VScaleExpr<VVecExpr<TA>, double> TRES;
  return TRES (a, b);
}


inline VVecExpr<VScaleExpr<VVecExpr<BaseVector>, Complex> >
operator* (const Complex & b, const BaseVector & a)
{
  typedef VScaleExpr<VVecExpr<BaseVector>, Complex> TRES;
  return TRES (a, b);
}

template <class TA>
inline VVecExpr<VScaleExpr<VVecExpr<TA>, Complex> >
operator* (const Complex & b, const VVecExpr<TA> & a)
{
  typedef VScaleExpr<VVecExpr<TA>, Complex> TRES;
  return TRES (a, b);
}



/* *********************** operator<< ********************** */

///
inline ostream & operator<< (ostream & ost, const BaseVector & v)
{
  return v.Print(ost);
}

///
inline double InnerProduct (const BaseVector & v1, const BaseVector & v2)
{
  return dynamic_cast<const S_BaseVector<double>&>(v1).InnerProduct(v2); 
}

///
template <class SCAL>
inline SCAL S_InnerProduct (const BaseVector & v1, const BaseVector & v2)
{
  return dynamic_cast<const S_BaseVector<SCAL>&>(v1).InnerProduct(v2); 
}

///
inline double L2Norm (const BaseVector & v)
{
  return v.L2Norm();
}