/*********************************************************************/
/* File:   basevector.cpp                                            */
/* Author: Joachim Schoeberl                                         */
/* Date:   7. Feb. 2003                                              */
/*********************************************************************/

/* 
   base class in vector hierarchy
*/

#include <la.hpp>
namespace ngla
{
  using namespace ngla;


  /*
  BaseVector :: BaseVector () throw()
  {
    ;
  }

  BaseVector :: ~BaseVector () throw()
  {
    ;
  }
  */

  BaseVector & BaseVector :: operator= (const BaseVector & v)
  {
    Set (1.0, v);
    return *this;
  }

  BaseVector & BaseVector :: operator= (double s)
  {
    SetScalar (s);
    return *this;
  }

  BaseVector & BaseVector :: operator= (Complex s)
  {
    SetScalar (s);
    return *this;
  }


  BaseVector & BaseVector :: Scale (double scal)
  {
    FVDouble() *= scal;
    return *this;
  }
  BaseVector & BaseVector :: Scale (Complex scal)
  {
    FVComplex() *= scal;
    return *this;
  }

  BaseVector & BaseVector :: SetScalar (double scal)
  {
    FVDouble() = scal;
    return *this;
  }

  BaseVector & BaseVector :: SetScalar (Complex scal)
  {
    FVComplex() = scal;
    return *this;
  }

  BaseVector & BaseVector :: Set (double scal, const BaseVector & v)
  {
    FVDouble() = scal * v.FVDouble();
    return *this;
  }

  BaseVector & BaseVector :: Set (Complex scal, const BaseVector & v)
  {
    FVComplex() = scal * v.FVComplex();
    return *this;
  }
    

  BaseVector & BaseVector :: Add (double scal, const BaseVector & v)
  {
    FVDouble() += scal * v.FVDouble();
    return *this;
  }

  BaseVector & BaseVector :: Add (Complex scal, const BaseVector & v)
  {
    FVComplex() += scal * v.FVComplex();
    return *this;
  }






  TempVector BaseVector :: Range (int begin, int end)
  {
    throw Exception ("BaseVector::Range called");
  }

  TempVector BaseVector ::Range (int begin, int end) const
  {
    throw Exception ("BaseVector::Range const called");
  }


  ostream & BaseVector :: Print (ostream & ost) const
  {
    throw Exception ("BaseVector::Print called");
  }
  
  void BaseVector :: Save(ostream & ost) const
  {
    FlatVector<double> fv = FVDouble();
    for (int i = 0; i < fv.Size(); i++)
      SaveBin (ost, fv(i));
  }

  void BaseVector :: Load(istream & ist) 
  {
    FlatVector<double> fv = FVDouble();
    for (int i = 0; i < fv.Size(); i++)
      LoadBin (ist, fv(i));
  }



  void BaseVector :: MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  { 
    ;
  }

  BaseVector * BaseVector :: CreateVector () const
  {
    cout << "Create vec called for base class" << endl;
    return 0;
  }


  void BaseVector :: SetRandom () 
  {
    FlatVector<double> fv = FVDouble();
    for (int i = 0; i < fv.Size(); i++)
      fv(i) = double (rand()) / RAND_MAX;
  }
  
  
  void BaseVector :: GetIndirect (const ARRAY<int> & ind, 
				  FlatVector<double> & v) const 
  { 
    FlatVector<double> fv = FVDouble();
    int es = EntrySize();
    int ii = 0;
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	{
	  int base = es * ind[i];
	  for (int j = 0; j < es; j++)
	    v[ii++] = fv[base++];
	}
      else
	{
	  for (int j = 0; j < es; j++)
	    v[ii++] = 0;
	}
  }
  
  void BaseVector :: GetIndirect (const ARRAY<int> & ind, 
				  FlatVector<Complex> & v) const 
  { 
    FlatVector<Complex> fv = FVComplex();
    int es = EntrySize() / 2;
    int ii = 0;
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	{
	  int base = es * ind[i];
	  for (int j = 0; j < es; j++)
	    v[ii++] = fv[base++];
	}
      else
	{
	  for (int j = 0; j < es; j++)
	    v[ii++] = 0;
	}
  }
  
  void BaseVector :: SetIndirect (const ARRAY<int> & ind, 
				  const FlatVector<double> & v) 
  { 
    FlatVector<double> fv = FVDouble();
    int es = EntrySize();
    int ii = 0;
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	{
	  int base = es * ind[i];
	  for (int j = 0; j < es; j++)
	    fv[base++] = v[ii++];
	}
      else
	ii += es;
  }

  void BaseVector :: SetIndirect (const ARRAY<int> & ind, 
				  const FlatVector<Complex> & v) 
  { 
    FlatVector<Complex> fv = FVComplex();
    int es = EntrySize() / 2;
    int ii = 0;
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	{
	  int base = es * ind[i];
	  for (int j = 0; j < es; j++)
	    fv[base++] = v[ii++];
	}
      else
	ii += es;
  }

  void BaseVector :: AddIndirect (const ARRAY<int> & ind, 
				  const FlatVector<double> & v) 
  { 
    FlatVector<double> fv = FVDouble();
    int es = EntrySize();
    int ii = 0;
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	{
	  int base = es * ind[i];
	  for (int j = 0; j < es; j++)
	    fv[base++] += v[ii++];
	}
      else
	ii += es;
  }

  void BaseVector :: AddIndirect (const ARRAY<int> & ind, 
				  const FlatVector<Complex> & v)
  { 
    FlatVector<Complex> fv = FVComplex();
    int es = EntrySize() / 2;
    int ii = 0;
    for (int i = 0; i < ind.Size(); i++)
      if (ind[i] != -1)
	{
	  int base = es * ind[i];
	  for (int j = 0; j < es; j++)
	    fv[base++] += v[ii++];
	}
      else
	ii += es;
  }
  

















/**
   Decision between double or Complex
 */

/*
template <class SCAL>
S_BaseVector<SCAL> :: S_BaseVector () throw()
{ 
  ;
}

template <class SCAL>
S_BaseVector<SCAL> :: ~S_BaseVector () throw()
{ 
  ;
}
*/


template <class SCAL>
S_BaseVector<SCAL> & S_BaseVector<SCAL> :: operator= (double s)
{
  SetScalar (s);
  return *this;
}

/*
template <class SCAL>
SCAL S_BaseVector<SCAL> :: InnerProduct (const BaseVector & v2) const
{
  throw Exception ("Inner Product called for S_BaseVector");
}
*/

template <class SCAL>
FlatVector<double> S_BaseVector<SCAL> :: FVDouble () const throw()
{
  return FlatVector<SCAL> (size * entrysize, Memory());
  /*
  FlatVector<SCAL> fv = FVScal();
  return FlatVector<SCAL> (fv.Size() * sizeof(SCAL)/sizeof(double),
			   reinterpret_cast<double*> (&fv(0)));
  */
}

template <class SCAL>
FlatVector<Complex> S_BaseVector<SCAL> :: FVComplex () const throw()
{
  throw Exception ("FVComplex called for real vector");
}



/*
S_BaseVector<Complex> :: S_BaseVector () throw()
{
  ;
}

S_BaseVector<Complex> :: ~S_BaseVector () throw()
{
  ;
}
*/

/*
Complex S_BaseVector<Complex> :: InnerProduct (const BaseVector & v2) const
{
  throw Exception ("Inner Product called for S_BaseVector");
}
*/

FlatVector<double> S_BaseVector<Complex> :: FVDouble () const throw()
{
  FlatVector<Complex> fv = FVScal();
  return FlatVector<double> (fv.Size() * sizeof(Complex)/sizeof(double),
			     reinterpret_cast<double*> (&fv(0)));
}


FlatVector<Complex> S_BaseVector<Complex> :: FVComplex () const throw()
{
  FlatVector<Complex> fv = FVScal();
  return FlatVector<Complex> (fv.Size() * sizeof(Complex)/sizeof(Complex),
			      reinterpret_cast<Complex*> (&fv(0)));
}


template class S_BaseVector<double>;
//SZ	template class S_BaseVector<Complex>;


}