//
// $Source: /cvsroot/gambit/gambit/sources/libgambit/dvector.imp,v $
// $Date: 2006/01/07 06:37:06 $
// $Revision: 1.2 $
//
// DESCRIPTION:
// Implementation of doubly-partitioned vector class
//
// This file is part of Gambit
// Copyright (c) 2002, The Gambit Project
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#include "dvector.h"

namespace Gambit {

//--------------------------------------------------------------------------
//          DVector<T>: Private and protected member functions
//--------------------------------------------------------------------------

template <class T>
int DVector<T>::sum(int part, const PVector<int> &v) const
{
  int s = 0;

  Array<int> len(v.Lengths());
  for (int j = 1; j <= len[part]; j++)
    s += v(part, j);

  return s;
}

template <class T> void DVector<T>::setindex(void)
{
  int index = 1;
  
  for (int i = 1; i <= dvlen.Length(); i++)  {
    dvptr[i] = this->svptr + index - 1;
    dvidx[i] = index;
    index += dvlen[i];
  }
}

template <class T> bool DVector<T>::Check(const DVector<T> &v) const
{
  for (int i = 1; i <= dvlen.Length(); i++)
    if (dvlen[i] != v.dvlen[i])   return false;
  return true;
}
    
//--------------------------------------------------------------------------
//    DVector<T>: Constructors, destructor, and constructive operators
//--------------------------------------------------------------------------

template <class T> DVector<T>::DVector(void) : dvptr(0)
{ }

template <class T> DVector<T>::DVector(const PVector<int> &sig)
  : PVector<T>((Array<int>)sig), dvlen(sig.Lengths().Length()),
    dvidx(sig.Lengths().Length())
{
  dvptr = new T **[dvlen.Length()];
  dvptr -= 1;

  for (int i = 1; i <= dvlen.Length(); i++)
    dvlen[i] = sig.Lengths()[i];

  setindex();
}

template <class T> DVector<T>::DVector(const Vector<T> &val,
					   const PVector<int> &sig)
  : PVector<T>(val, sig), dvlen(sig.Lengths().Length()),
    dvidx(sig.Lengths().Length())
{
  dvptr = new T **[dvlen.Length()];
  dvptr -= 1;

  for (int i = 1; i <= dvlen.Length(); i++)
    dvlen[i] = sig.Lengths()[i];

  setindex();
}

template <class T> DVector<T>::DVector(const DVector<T> &v)
  : PVector<T>(v), dvlen(v.dvlen), dvidx(v.dvidx)
{
  dvptr = new T **[dvlen.Length()];
  dvptr -= 1;

  setindex();
}

template <class T> DVector<T>::~DVector()
{
  if (dvptr)  delete [] (dvptr + 1);
}

template <class T> DVector<T> &DVector<T>::operator=(const DVector<T> &v)
{
  if (!Check(v)) {
    throw DimensionException();
  }

  PVector<T>::operator=(v);
  return *this;
}

template <class T> DVector<T> &DVector<T>::operator=(const PVector<T> &v)
{
  PVector<T>::operator=(v);
  return *this;
}

template <class T> DVector<T> &DVector<T>::operator=(const Vector<T> &v)
{
  PVector<T>::operator=(v);
  return *this;
}

template <class T> DVector<T> &DVector<T>::operator=(T c)
{
  PVector<T>::operator=(c);
  return *this;
}



//--------------------------------------------------------------------------
//                    DVector<T>: Operator definitions
//--------------------------------------------------------------------------

template <class T> T &DVector<T>::operator()(int a, int b, int c)
{
  if (dvlen.First() > a || a > dvlen.Last()) {
    throw IndexException();
  }
  if (1 > b || b > dvlen[a]) {
    throw IndexException();
  }
  if (1 > c || c > this->svlen[dvidx[a] + b - 1]) {
    throw IndexException();
  }
  return dvptr[a][b][c];
}

template <class T> const T &DVector<T>::operator()(int a, int b, int c) const
{
  if (dvlen.First() > a || a > dvlen.Last()) {
    throw IndexException();
  }
  if (1 > b || b > dvlen[a]) {
    throw IndexException();
  }
  if (1 > c || c > this->svlen[dvidx[a] + b - 1]) {
    throw IndexException();
  }
  return dvptr[a][b][c];
}

template <class T>
DVector<T> DVector<T>::operator+(const DVector<T> &v) const
{
  if (!Check(v)) {
    throw DimensionException();
  }

  DVector<T> tmp(*this);
  tmp.PVector<T>::operator+=(v);
  return tmp;
}

template <class T>
DVector<T> &DVector<T>::operator+=(const DVector<T> &v)
{
  if (!Check(v)) {
    throw DimensionException();
  }
  PVector<T>::operator+=(v);
  return *this;
}

template <class T> DVector<T> DVector<T>::operator-(void) const
{
  DVector<T> tmp(*this);
  for (int i = this->First(); i <= this->Last(); i++)
    tmp[i] = -tmp[i];
  return tmp;
}

template <class T>
DVector<T> DVector<T>::operator-(const DVector<T> &v) const
{
  if (!Check(v)) {
    throw DimensionException();
  }
  DVector<T> tmp(*this);
  tmp.PVector<T>::operator-=(v);
  return tmp;
}

template <class T>
DVector<T> &DVector<T>::operator-=(const DVector<T> &v)
{
  if (!Check(v)) {
    throw DimensionException();
  }
  PVector<T>::operator-=(v);
  return *this;
}

template <class T> T DVector<T>::operator*(const DVector<T> &v) const
{
  if (!Check(v)) {
    throw DimensionException();
  }
  return (*this).PVector<T>::operator*(v);
}

template <class T> DVector<T> &DVector<T>::operator*=(const T &c)
{
  PVector<T>::operator*=(c);
  return *this;
}

template <class T> DVector<T> DVector<T>::operator/(const T &c) const
{
  DVector<T> tmp(*this);
  tmp = tmp.PVector<T>::operator/(c);
  return tmp;
}

template <class T> bool DVector<T>::operator==(const DVector<T> &v) const
{
  if (!Check(v)) {
    throw DimensionException();
  }
  return PVector<T>::operator==(v);
}

template <class T> bool DVector<T>::operator!=(const DVector<T> &v) const
{
  return !(*this == v);
}

//-------------------------------------------------------------------------
//                 DVector<T>: General data access
//-------------------------------------------------------------------------

template <class T>
void DVector<T>::CopySubRow(int row, int col, const DVector<T> &v)
{
  if (!Check(v)) {
    throw DimensionException();
  }
  if (dvlen.First() > row || row > dvlen.Last()) {
    throw IndexException();
  }
  if (1 > col || col > dvlen[row]) {
    throw IndexException();
  }

  for (int i = 1; i <= this->svlen[dvidx[row]+col-1]; i++)
    dvptr[row][col][i] = v.dvptr[row][col][i];
}



template <class T> const Array<int> &DVector<T>::DPLengths(void) const
{
  return dvlen;
}

} // end namespace Gambit