//
// $Source: /cvsroot/gambit/gambit/sources/tools/lp/lptab.imp,v $
// $Date: 2006/01/07 06:37:34 $
// $Revision: 1.7 $
//
// DESCRIPTION:
// Implementation of LP tableaus
//
// 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 "lptab.h"

// ---------------------------------------------------------------------------
//                   LPTableau member definitions 
// ---------------------------------------------------------------------------


template <class T>
LPTableau<T>::LPTableau(const Gambit::Matrix<T> &A, const Gambit::Vector<T> &b)
  : Tableau<T>(A,b), dual(A.MinRow(),A.MaxRow()),
    unitcost(A.MinRow(),A.MaxRow()), cost(A.MinCol(),A.MaxCol())
{ };

template <class T>
LPTableau<T>::LPTableau(const Gambit::Matrix<T> &A, const Gambit::Array<int> &art, 
			const Gambit::Vector<T> &b)
  : Tableau<T>(A,art,b), dual(A.MinRow(),A.MaxRow()),
    unitcost(A.MinRow(),A.MaxRow()), cost(A.MinCol(),A.MaxCol()+art.Length())
{ };

template <class T>
LPTableau<T>::LPTableau(const LPTableau<T> &orig)
  : Tableau<T>(orig), dual(orig.dual),  unitcost(orig.unitcost),
    cost(orig.cost)
{ }

template <class T>
LPTableau<T>::~LPTableau()
{ }


template <class T>
LPTableau<T>& LPTableau<T>::operator=(const LPTableau<T> &orig)
{
  Tableau<T>::operator=(orig);
  if(this!= &orig) {
    dual=orig.dual;
    unitcost= orig.unitcost;
    cost= orig.cost;
  }
  return *this;
}

// cost-based functions

template<>
void LPTableau<Gambit::Rational>::SetCost(const Gambit::Vector<Gambit::Rational>& c)
{
  int i;
  if(cost.First()==c.First() && cost.Last()==c.Last()) {
    for(i=cost.First();i<=cost.Last();i++) cost[i] = c[i]*(Gambit::Rational)Tableau<Gambit::Rational>::TotDenom();
    for(i=unitcost.First();i<=unitcost.Last();i++) unitcost[i] = (Gambit::Rational)0;
    Refactor();
    SolveDual();
    return;
  }
  if(c.First()!=cost.First()) throw Gambit::DimensionException();
  if(c.Last()!=(cost.Last()+unitcost.Length())) throw Gambit::DimensionException();
  for(i=c.First();i<=cost.Last();i++)
    cost[i]=c[i]*(Gambit::Rational)Tableau<Gambit::Rational>::TotDenom();
  for(i=unitcost.First();i<=unitcost.Last();i++)
    unitcost[i]=c[cost.Length()+i-unitcost.First()+1];
  // gout << "\nc: " << c.First() << " " << c.Last() << " " << c;
  // gout << "\ncost: " << cost.First() << " " << cost.Last() << " " << cost;
  // gout << "\nunit: " << unitcost.First() << " " << unitcost.Last() << " " << unitcost;
  //** added for Gambit::Rational
  Refactor();
  SolveDual();
}

template<>
void LPTableau<double>::SetCost(const Gambit::Vector<double>& c)
{
  int i;
  if(cost.First()==c.First() && cost.Last()==c.Last()) {
    for(i=cost.First();i<=cost.Last();i++) cost[i] = c[i];
    for(i=unitcost.First();i<=unitcost.Last();i++) unitcost[i] = (double)0;
    Refactor();
    SolveDual();
    return;
  }
  if(c.First()!=cost.First()) throw Gambit::DimensionException();
  if(c.Last()!=(cost.Last()+unitcost.Length())) throw Gambit::DimensionException();
  for(i=c.First();i<=cost.Last();i++)
    cost[i]=c[i];
  for(i=unitcost.First();i<=unitcost.Last();i++)
    unitcost[i]=c[cost.Length()+i-unitcost.First()+1];
  // gout << "\nc: " << c.First() << " " << c.Last() << " " << c;
  // gout << "\ncost: " << cost.First() << " " << cost.Last() << " " << cost;
  // gout << "\nunit: " << unitcost.First() << " " << unitcost.Last() << " " << unitcost;
  //** added for Gambit::Rational
  Refactor();
  SolveDual();
}

template <class T>
void LPTableau<T>::SetCost(const Gambit::Vector<T> &uc, const Gambit::Vector<T> &c)
{
  if(cost.First()!=c.First() || cost.Last()!=c.Last()) 
    throw Gambit::DimensionException();
  if(unitcost.First()!=uc.First() || unitcost.Last()!=uc.Last()) 
    throw Gambit::DimensionException();
  int i;
  for(i=cost.First();i<=cost.Last();i++) cost[i]=c[i];
  for(i=unitcost.First();i<=unitcost.Last();i++) unitcost[i]=uc[i];
  SolveDual();
}


template <class T>
Gambit::Vector<T> LPTableau<T>::GetCost(void) const
{
  Gambit::Vector<T> x(cost.First(),cost.Last());
  for(int i=x.First();i<=x.Last();i++)x[i]=cost[i];
  return x;
  //  return cost; 
}


template <class T>
Gambit::Vector<T> LPTableau<T>::GetUnitCost(void) const
{
  Gambit::Vector<T> x(unitcost.First(),unitcost.Last());
  for(int i=x.First();i<=x.Last();i++)x[i]=unitcost[i];
  return x;
}

template<> double LPTableau<double>::TotalCost(void)
{
  Gambit::Vector<double> tmpcol(MinRow(),MaxRow());
  BasisSelect(unitcost,cost,tmpcol);
  return tmpcol*solution;
}

template<> Gambit::Rational LPTableau<Gambit::Rational>::TotalCost(void)
{
  Gambit::Vector<Gambit::Rational> tmpcol(MinRow(),MaxRow());
  Gambit::Vector<Gambit::Rational> sol(MinRow(),MaxRow());
  BasisSelect(unitcost,cost,tmpcol);
  BasisVector(sol);
  for(int i=tmpcol.First();i<=tmpcol.Last();i++)
    if(Label(i)>0)tmpcol[i]/=(Gambit::Rational)Tableau<Gambit::Rational>::TotDenom();

  return tmpcol*sol;
}

template<>
void LPTableau<double>::DualVector(Gambit::Vector<double> &L) const
{
  L= dual;
}

template<>
void LPTableau<Gambit::Rational>::DualVector(Gambit::Vector<Gambit::Rational> &out) const
{
  out= dual;
  //  for(int i=out.First();i<=out.Last();i++) 
  //  if(Label(i)>=0) out[i]*=TotDenom();
}

template <class T>
T LPTableau<T>::RelativeCost(int col) const
{
  Gambit::Vector<T> tmpcol(this->MinRow(),this->MaxRow());
  if( col<0 ) {
    return unitcost[-col] - dual[-col];
  }
  else {
    GetColumn(col, (Gambit::Vector<T> &)tmpcol);
    return cost[col] - dual*tmpcol;
  }
}

/*

template <class T>
void LPTableau<T>::RelativeCostVector(Gambit::Vector<T> &relunitcost,
				      Gambit::Vector<T> &relcost) const
{

  if(!A->CheckColumn(relunitcost)) throw Gambit::DimensionException();
  if(!A->CheckRow(relcost)) throw Gambit::DimensionException();
  
  relunitcost= unitcost - dual;
  relcost= cost - dual*A;  // pre multiplication not defined?  
}
*/


template <class T>
void LPTableau<T>::SolveDual()
{
  Gambit::Vector<T> tmpcol1(this->MinRow(),this->MaxRow());
  BasisSelect(unitcost,cost,tmpcol1);
  SolveT(tmpcol1,dual);
}

// Redefined functions

template <class T>
void LPTableau<T>::Refactor()
{
  // gout << "\nIn LPTableau<T>::Refactor()";  
  Tableau<T>::Refactor();
  SolveDual();
}

template <class T>
void LPTableau<T>::Pivot(int outrow,int col)
{
  // gout << "\nIn LPTableau<T>::Pivot() ";
  // gout << "outrow: " << outrow << " col: " << col;
  // BigDump(gout);
  Tableau<T>::Pivot(outrow,col);
  SolveDual();
}

template <class T>
void LPTableau<T>::ReversePivots(Gambit::List<Gambit::Array<int> > &PivotList)
{
  Gambit::Vector<T> tmpcol(this->MinRow(),this->MaxRow());
  // gout << "\nIn LPTableau<T>::ReversePivots";
  bool flag;
  int i,j,k,enter;
  T ratio,a_ij,a_ik,b_i,b_k,c_j,c_k,c_jo,x;
  Gambit::List<int> BestSet;
  Gambit::Array<int> pivot(2);
  Gambit::Vector<T> tmpdual(this->MinRow(),this->MaxRow());

  Gambit::Vector<T> solution(tmpcol);  //$$
  BasisVector(solution);        //$$

  // BigDump(gout);
  // gout << "\ncost: " << GetCost();
  // gout << "\nunitcost: " << GetUnitCost() << "\n";
  // for(i=MinCol();i<=MaxCol();i++) gout << " " << RelativeCost(i);
  // for(i=MinRow();i<=MaxRow();i++) gout << " " << RelativeCost(-i);

  for(j=-this->MaxRow();j<=this->MaxCol();j++) if(j && !this->Member(j)  && !this->IsBlocked(j)) {
    SolveColumn(j,tmpcol);
    // gout << "\nColumn " << j;
    // gout << "\nPivCol = " << tmpcol;
    // gout << "\ncurrentSolCol = " << solution;
    
    // find all i where prior tableau is primal feasible
    
    BestSet = Gambit::List<int>();
    for(i=this->MinRow();i<=this->MaxRow();i++)
      if(GtZero(tmpcol[i])) BestSet.Append(i);
    if(BestSet.Length()>0) {
      ratio = solution[BestSet[1]]/tmpcol[BestSet[1]];
      // find max ratio
      for(i=2;i<=BestSet.Length();i++) {
	x = solution[BestSet[i]]/tmpcol[BestSet[i]];
	if(GtZero(x-ratio)) ratio = x;
      }
      // eliminate nonmaximizers
      for(i=BestSet.Length();i>=1;i--) {
	x = solution[BestSet[i]]/tmpcol[BestSet[i]];
	if(LtZero(x-ratio)) BestSet.Remove(i);
      }	

      // check that j would be the row to exit in prior tableau

      // first check that prior pivot entry > 0 
      for(i=BestSet.Length();i>=1;i--) {
	a_ij = (T)1/tmpcol[BestSet[i]];
	if(LeZero(a_ij)) {
	  // gout << "\nj not row to exit in prior tableau: a_ij <= 0";
	  BestSet.Remove(i);
	}
	else {
	  // next check that prior pivot entry attains max ratio
	  b_i = solution[BestSet[i]]/tmpcol[BestSet[i]];
	  ratio = b_i/a_ij;
  
	  flag = 0;
	  for(k=tmpcol.First();k<=tmpcol.Last() && !flag;k++) 
	    if(k!=BestSet[i]) {
	      a_ik = - a_ij * tmpcol[k];
	      b_k = solution[k] - b_i*tmpcol[k];
	      if(GtZero(a_ik) && GtZero(b_k/a_ik -ratio)) {
		// gout << "\nj not row to exit in prior tableau: ";
		// gout << "higher ratio at row= " << k;
		BestSet.Remove(i);
		flag = 1;
	      }
	      else if(GtZero(a_ik) && EqZero(b_k/a_ik-ratio) && this->Label(k)<j) {
		// gout << "\nj not row to exit in prior tableau: ";
		// gout << "same ratio,lower lex at k= " << k;
		BestSet.Remove(i);
		flag = 1;
	      }
	    }
	}
      }
    }
    // gout << "\nafter checking rows, BestSet = ";
    // BestSet.Dump(gout);

    // check that i would be the column to enter in prior tableau

    for(i=BestSet.Length();i>=1;i--) {
      enter = this->Label(BestSet[i]);
      // gout << "\nenter = " << enter;
      
      tmpcol = (T)0;
      tmpcol[BestSet[i]]=(T)1;
      // gout << "\ntmpcol, loc 1: " << tmpcol;
      SolveT(tmpcol,tmpdual);
      // gout << "\ntmpcol, loc 2: " << tmpcol;
      // gout << "\ntmpdual, loc 1: " << tmpdual;
      
/*      if( j<0 )
	{ tmpcol=(T)0; tmpcol[-j]=(T)1; }
      else
	A->GetColumn(j,tmpcol);
*/
      GetColumn(j,tmpcol);      
      // gout << "\ncol " << j << ": " << tmpcol;
      a_ij = tmpdual*tmpcol;
      c_j = RelativeCost(j);
      if(EqZero(a_ij)) {
	// gout << "\ni not col to enter in prior tableau: ";
	// gout << "a_ij=0";
	BestSet.Remove(i);
      }
      else {
	ratio = c_j/a_ij;
	// gout << " ratio: " << ratio;
	if(enter<0) 
	  a_ik = tmpdual[-enter];
	else {
	  GetColumn(enter,tmpcol);
//	  A->GetColumn(enter,tmpcol);
	  a_ik = tmpdual*tmpcol;
	}
	c_k = RelativeCost(enter);
	c_jo = c_k - a_ik * ratio; 
	// gout << "\ntmpdual = " << tmpdual << "\n";
	// gout << " c_j:" << c_j; 
	// gout << " c_k:" << c_k; 
	// gout << " c_jo:" << c_jo; 
	// gout << " a_ij:" << a_ij; 
	// gout << " a_ik:" << a_ik; 
	if(GeZero(c_jo)) {
	  // gout << "\ni not col to enter in prior tableau: ";
	  // gout << "c_jo<0";
	  BestSet.Remove(i);
	}
	else {
	  flag=0;
	  for(k=-this->b->Last();k<enter && !flag;k++) if(k!=0) {
	    if(k<0) 
	      a_ik=tmpdual[-k];
	    else {
//	      A->GetColumn(k,tmpcol);
	      GetColumn(k,tmpcol);
	      a_ik = tmpdual*tmpcol;
	    }
	    c_k = RelativeCost(k);
	    c_jo = c_k - a_ik * ratio; 
	    
	    if(LtZero(c_jo)) { 
	      // gout << "\ni not col to enter in prior tableau: ";
	      // gout << "c_jo < 0 for k = " << k;
	      BestSet.Remove(i);
	      flag=1;
	    }
	  }
	}
      }
    }
    // gout << "\nafter checking cols, BestSet = ";
    // BestSet.Dump(gout);

    if(BestSet.Length()>0) 
      for(i=1;i<=BestSet.Length();i++) {
	pivot[1] = BestSet[i];
	pivot[2] = j;
	PivotList.Append(pivot);
      }
  }
}

template <class T>
bool LPTableau<T>::IsReversePivot(int i, int j)
{
  Gambit::Vector<T> tmpcol(this->MinRow(),this->MaxRow());

  // first check that pivot preserves primal feasibility
  
  // gout << "\nin IsReversePivot, i= " << i << " j = "<< j;
  SolveColumn(j,tmpcol);
  Gambit::Vector<T> solution(tmpcol);  //$$
  BasisVector(solution);        //$$
  // gout << "\ncurrentPivCol = " << tmpcol;
  // gout << "\ncurrentSolCol = " << solution;
  if(LeZero(tmpcol[i])) { 
    // gout << "\nPrior tableau not primal feasible: currentPivCol[i] <= 0";
    return 0;
  }
  int k;
  T ratio = solution[i]/tmpcol[i];
  // gout << "\nratio = " << ratio;
  
  for(k=tmpcol.First();k<=tmpcol.Last();k++)
    if(GtZero(tmpcol[k]) && GtZero(solution[k]/tmpcol[k]-ratio)) {
      // gout << "\nPrior tableau not primal feasible: i not min ratio";
      return 0;
    }
  // check that j would be the row to exit in prior tableau
  
  T a_ij,a_ik,b_i,b_k,c_j,c_k,c_jo;

  a_ij = (T)1/tmpcol[i];
  if(LeZero(a_ij)) {
    // gout << "\nj not row to exit in prior tableau: a_ij <= 0";
    return 0;
  }
  b_i = solution[i]/tmpcol[i];
  ratio = b_i/a_ij;
  
  for(k=tmpcol.First();k<=tmpcol.Last();k++) 
    if(k!=i) {
      a_ik = - a_ij * tmpcol[k];
      b_k = solution[k] - b_i*tmpcol[k];
      if(GtZero(a_ik) && GtZero(b_k/a_ik -ratio)) {
	// gout << "\nj not row to exit in prior tableau: ";
	// gout << "higher ratio at row= " << k;
	return 0;
      }
      if(GtZero(a_ik) && EqZero(b_k/a_ik-ratio) && this->Label(k)<j) {
	// gout << "\nj not row to exit in prior tableau: ";
	// gout << "same ratio,lower lex at k= " << k;
	return 0;
      }
    }

  // check that i would be the column to enter in prior tableau
  
  int enter = this->Label(i);
  // gout << "\nenter = " << enter;
  
  Gambit::Vector<T> tmpdual(this->MinRow(),this->MaxRow());
  tmpcol = (T)0;
  tmpcol[i]=(T)1;
  SolveT(tmpcol,tmpdual);

/*
  if( j<0 )
    { tmpcol=(T)0; tmpcol[-j]=(T)1; }
  else
    A->GetColumn(j,tmpcol);
*/
    GetColumn(j,tmpcol);
  
  // gout << "\ncol j = " << tmpcol;
  a_ij = tmpdual*tmpcol;
  c_j = RelativeCost(j);
  if(EqZero(a_ij)) {
    // gout << "\ni not col to enter in prior tableau: ";
    // gout << "a_ij=0";
    return 0;
  }
  ratio = c_j/a_ij;
  
  if(enter<0) 
    a_ik = tmpdual[-enter];
  else {
//    A->GetColumn(enter,tmpcol);
    GetColumn(enter,tmpcol);
    a_ik = tmpdual*tmpcol;
  }
  c_k = RelativeCost(enter);
  c_jo = c_k - a_ik * ratio; 
  if(GeZero(c_jo)) {
    // gout << "\ni not col to enter in prior tableau: ";
    // gout << "c_jo<0";
    return 0;
  }

  for(k=-this->b->Last();k<enter;k++) if(k!=0) {
    if(k<0) 
      a_ik=tmpdual[-k];
    else {
//      A->GetColumn(k,tmpcol);
      GetColumn(k,tmpcol);
      a_ik = tmpdual*tmpcol;
    }
    c_k = RelativeCost(k);
    c_jo = c_k - a_ik * ratio; 
    
    if(LtZero(c_jo)) { 
      // gout << "\ni not col to enter in prior tableau: ";
      // gout << "c_jo < 0 for k = " << k;
      return 0;
    }
  }
  // gout << "\nValid Reverse pivot at i = " << i << " j =  " << j;
  return 1;
}

template <class T>
void LPTableau<T>::DualReversePivots(Gambit::List<Gambit::Array<int> > &/*list*/)
{
}

template <class T>
bool LPTableau<T>::IsDualReversePivot(int i, int j)
{
  // first check that pivot preserves dual feasibility
  
  // gout << "\nin IsDualReversePivot, i= " << i << " j = "<< j;

  int k;

  Gambit::Vector<T> tmpcol(this->MinRow(),this->MaxRow());
  Gambit::Vector<T> tmpdual(this->MinRow(),this->MaxRow());
  tmpcol = (T)0;
  tmpcol[i]=(T)1;
  SolveT(tmpcol,tmpdual);

  Gambit::Vector<T> solution(tmpcol);  //$$
  BasisVector(solution);        //$$

  // gout << "\ncurrentPivCol = " << tmpcol;
  // gout << "\ncurrentSolCol = " << solution;

  T a_ij,a_ik,c_j,c_k,ratio;

/*  if( j<0 )
    { tmpcol=(T)0; tmpcol[-j]=(T)1; }
  else
    A->GetColumn(j,tmpcol);
  */

    GetColumn(j,tmpcol);
  
  a_ij = tmpdual*tmpcol;
  c_j = RelativeCost(j);
  if(GeZero(a_ij)) {
    // gout << "\nPrior tableau not dual feasible: ";
    // gout << "a_ij>=0";
    return 0;
  }
  ratio = c_j/a_ij;
  
  for(k=-this->b->Last();k<=cost.Last();k++) if(k!=0) {
    if(k<0) 
      a_ik=tmpdual[-k];
    else {
//      A->GetColumn(k,tmpcol);
      GetColumn(k,tmpcol);
      a_ik = tmpdual*tmpcol;
    }
    c_k = RelativeCost(k);
    
    if(LtZero(a_ik) && GtZero(c_k/a_ik-ratio)) { 
      // gout << "\nPrior tableau not dual feasible: ";
      // gout << "\nhigher ratio for k = " << k;
      return 0;
    }
  }

  // check that i would be the column to enter in prior tableau

  int enter = this->Label(i);
  // gout << "\nenter = " << enter;

  if(enter<0) 
    a_ik = tmpdual[-enter];
  else {
//    A->GetColumn(enter,tmpcol);
    GetColumn(enter,tmpcol);
    a_ik = tmpdual*tmpcol;
  }
  a_ik = a_ik/a_ij;
  c_k = RelativeCost(enter);
  c_k -= a_ik * c_j; 

  if(GeZero(a_ik)) {
    // gout << "\ni not col to enter in prior tableau: ";
    // gout << "a_ik>=0";
    return 0;
  }
  ratio = c_k/a_ik;

  for(k=-this->b->Last();k<=cost.Last();k++) if(k!=0) {
    if(k<0) 
      a_ik=tmpdual[-k];
    else {
//    A->GetColumn(k,tmpcol);
    GetColumn(k,tmpcol);
    a_ik = tmpdual*tmpcol;
    }
    a_ik = a_ik/a_ij;
    c_k = RelativeCost(k);
    c_k -= a_ik * c_j; 
    
    if(LtZero(a_ik) && GtZero(c_k/a_ik- ratio)) { 
      // gout << "\ni not col to enter in prior tableau: ";
      // gout << "\nhigher ratio for k = " << k;
      return 0;
    }
    if(k<enter && LtZero(a_ik) && EqZero(c_k/a_ik - ratio)) { 
      // gout << "\ni not col to enter in prior tableau: ";
      // gout << "\nsame ratio and lower lex for k = " << k;
      return 0;
    }
  }

  // check that j would be the row to exit in prior tableau

  SolveColumn(j,tmpcol);
  // gout << "\ncurrentPivCol = " << tmpcol;
  // gout << "\ncurrentSolCol = " << solution;

  T b_k,b_i;
  b_i= solution[i]/tmpcol[i];
  if(LeZero(b_i)) {
    // gout << "\nj not row to exit in prior tableau: ";
    // gout << "b_i<=0";
    return 0;
    }

  
  for(k=this->b->First();k<=this->b->Last();k++) 
    if(k!=i) {
      b_k = solution[k] -  b_i * tmpcol[k];
      if(GtZero(b_k) && this->Label(k)<j) {
	// gout << "\nj not row to exit in prior tableau: ";
	// gout << "same ratio,lower lex at k= " << k;
	return 0;
      }
    }
  // gout << "\nValid Reverse pivot at i = " << i << " j =  " << j;
  return 1;
}

/*
template <class T>
BFS<T> LPTableau<T>::DualBFS() const
{
  BFS<T> cbfs((T) 0);
  for(int i=MinRow();i<=MaxRow();i++) {
    if(!Member(-i)) 
      cbfs.Define(-i,dual[i]);
  }
  return cbfs;
}
*/

template <class T>
BFS<T> LPTableau<T>::DualBFS() const
{
  BFS<T> cbfs((T) 0);
  Gambit::Vector<T> d(this->MinRow(),this->MaxRow());
  DualVector(d);

  for(int i=this->MinRow();i<=this->MaxRow();i++) {
    if(!this->Member(-i)) 
      cbfs.Define(-i,d[i]);
  }
  // gout << "\ndual: " << d;
  return cbfs;
}

template <class T>
int LPTableau<T>::LastLabel( void )
{
  return this->artificial.Last();
}

template <class T>
void LPTableau<T>::BasisSelect(const Gambit::Array<T> &rowv, Gambit::Vector<T> &colv) const
{
  for(int i=this->basis.First(); i<=this->basis.Last(); i++) {
    if(this->basis.Label(i)<0)
      colv[i]= 0;
    else
      colv[i]= rowv[this->basis.Label(i)];
  }
}

template <class T>
void LPTableau<T>::BasisSelect(const Gambit::Array<T> &unitv,
			   const Gambit::Array<T> &rowv,
			   Gambit::Vector<T> &colv ) const
{
  for(int i=this->basis.First(); i<=this->basis.Last(); i++) {
    if(this->basis.Label(i)<0)
      colv[i]= unitv[-this->basis.Label(i)];
    else
      colv[i]= rowv[this->basis.Label(i)];
  }
}

template <class T> LPTableau<T>::BadPivot::~BadPivot()
{ }

template <class T> std::string LPTableau<T>::BadPivot::GetDescription(void) const
{
  return "Bad Pivot in LPTableau";
}