tableau.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/numerical/tableau.imp,v $
// $Date: 2002/09/26 17:50:57 $
// $Revision: 1.1.2.1 $
//
// DESCRIPTION:
// Implementation of tableau 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 "tableau.h"
#include "base/base.h"

// ---------------------------------------------------------------------------
//                   Tableau<double> method definitions
// ---------------------------------------------------------------------------

// Constructors and Destructor
 
Tableau<double>::Tableau(const gMatrix<double> &A, const gVector<double> &b)
  : TableauInterface<double>(A,b), B(*this), tmpcol(b.First(),b.Last())
{
  Solve(b, solution);
};

Tableau<double>::Tableau(const gMatrix<double> &A, const gBlock<int> &art, 
			 const gVector<double> &b)
  : TableauInterface<double>(A,art,b), B(*this), tmpcol(b.First(),b.Last())
{
  Solve(b, solution);
};


Tableau<double>::Tableau(const Tableau<double> &orig)
  : TableauInterface<double>(orig), B(orig.B,*this), tmpcol(orig.tmpcol)
{ }

Tableau<double>::~Tableau()
{ }

Tableau<double>& Tableau<double>::operator=(const Tableau<double> &orig)
{
  TableauInterface<double>::operator=(orig);
  if(this!= &orig) {
    B.Copy(orig.B,*this);
    tmpcol = orig.tmpcol;
  }
  return *this;
}

//
// pivoting operations
//

int Tableau<double>::CanPivot(int outlabel, int col)
{
  SolveColumn(col,tmpcol);
  double val = tmpcol[basis.Find(outlabel)];
  if(val <=eps2 && val >= -eps2) return 0;
  return 1;  
}

void Tableau<double>::Pivot(int outrow,int col)
{
  if(!RowIndex(outrow) || !ValidIndex(col)) throw BadPivot();

  // int outlabel = Label(outrow);
  // gout << "\noutrow:" << outrow;
  // gout << " outlabel: " << outlabel;
  // gout << " inlabel: " << col;
  // BigDump(gout);
  basis.Pivot(outrow,col);
  
  B.update(outrow, col);
  Solve(*b, solution);
  npivots++;
  // BigDump(gout);
}

void Tableau<double>::SolveColumn(int col, gVector<double> &out)
{
  //** can we use tmpcol here, instead of allocating new vector?
  gVector<double> tmpcol2(MinRow(),MaxRow());
  GetColumn(col,tmpcol2);
  Solve(tmpcol2,out);
}

void Tableau<double>::BasisVector(gVector<double> &out) const
{
  out= solution;
}

//
// raw Tableau functions
//

void Tableau<double>::Refactor()
{
  B.refactor();
  //** is re-solve necessary here?
  Solve(*b, solution);
}

void Tableau<double>::SetRefactor(int n)
{
  B.SetRefactor(n);
}

void Tableau<double>::SetConst(const gVector<double> &bnew)
{
  if(bnew.First()!=b->First() || bnew.Last()!=b->Last())
    throw BadDim();
  b=&bnew;
  Solve(*b, solution);
}

//** this function is not currently used.  Drop it?
void Tableau<double>::SetBasis(const Basis &in)
{
  basis= in;
  B.refactor();
  Solve(*b, solution);
}

void Tableau<double>::Solve(const gVector<double> &b, gVector<double> &x)
{
  B.solve(b,x);
}

void Tableau<double>::SolveT(const gVector<double> &c, gVector<double> &y)
{
  B.solveT(c,y);
  //** gout << "\nTableau<double>::SolveT(), y: " << y;
  //   gout << "\nc: " << c;
}

bool Tableau<double>::IsFeasible()
{
  //** is it really necessary to solve first here?
  Solve(*b, solution);
  for(int i=solution.First();i<=solution.Last();i++)
    if(solution[i]>=eps2) return false;
  return true;
}

bool Tableau<double>::IsLexMin()
{
  int i,j;
  for(i=MinRow();i<=MaxRow();i++)
    if(EqZero(solution[i]))
      for(j=-MaxRow();j<Label(i);j++) if(j!=0){
	SolveColumn(j,tmpcol);
	if(LtZero(tmpcol[i]))
	  return 0;
      }
  return 1;
}

gOutput &operator<<(gOutput &to, const Tableau<double> &v)
{
  v.Dump(to); return to;
}

// ---------------------------------------------------------------------------
//                   Tableau<gRational> method definitions
// ---------------------------------------------------------------------------


gInteger find_lcd(const gMatrix<gRational> &mat)
{
  gInteger lcd(1);
  for(int i=mat.MinRow();i<=mat.MaxRow();i++)
    for(int j=mat.MinCol();j<=mat.MaxCol();j++) 
      lcd = lcm(mat(i,j).denominator(),lcd);
  return lcd;
}

gInteger find_lcd(const gVector<gRational> &vec)
{
  gInteger lcd(1);
  for(int i=vec.First();i<=vec.Last();i++)
    lcd = lcm(vec[i].denominator(),lcd);
  return lcd;
}

// Constructors and Destructor
 
Tableau<gRational>::Tableau(const gMatrix<gRational> &A, 
			    const gVector<gRational> &b) 
  : TableauInterface<gRational>(A,b), 
    Tabdat(A.MinRow(),A.MaxRow(),A.MinCol(),A.MaxCol()),
    Coeff(b.First(),b.Last()), denom(1), tmpcol(b.First(),b.Last()), 
    nonbasic(A.MinCol(),A.MaxCol())
{
  int j;
  for(j=MinCol();j<=MaxCol();j++) 
    nonbasic[j] = j;
  
  totdenom = lcm(find_lcd(A),find_lcd(b));
  if(totdenom<=0) throw BadDenom();
  
  for (int i = b.First();i<=b.Last();i++) {
    gRational x = b[i]*(gRational)totdenom;
    if(x.denominator() != 1) throw BadDenom();
    Coeff[i] = x.numerator();
  }
  for (int i = MinRow();i<=MaxRow();i++) 
    for (int j = MinCol();j<=MaxCol();j++) {
      gRational x = A(i,j)*(gRational)totdenom;
      if(x.denominator() != 1) throw BadDenom();
      Tabdat(i,j) = x.numerator();
    }
  for (int i = b.First();i<=b.Last();i++) 
    solution[i] = (gRational)Coeff[i];
}

Tableau<gRational>::Tableau(const gMatrix<gRational> &A, 
			    const gBlock<int> &art, 
			    const gVector<gRational> &b) 
  : TableauInterface<gRational>(A,art,b), 
    Tabdat(A.MinRow(),A.MaxRow(),A.MinCol(),A.MaxCol()+art.Length()),
    Coeff(b.First(),b.Last()), denom(1), tmpcol(b.First(),b.Last()), 
    nonbasic(A.MinCol(),A.MaxCol()+art.Length())
{
  int j;
  for(j=MinCol();j<=MaxCol();j++) 
    nonbasic[j] = j;
  
  totdenom = lcm(find_lcd(A),find_lcd(b));
  if(totdenom<=0) throw BadDenom();
  
  for (int i = b.First();i<=b.Last();i++) {
    gRational x = b[i]*(gRational)totdenom;
    if(x.denominator() != 1) throw BadDenom();
    Coeff[i] = x.numerator();
  }
  for (int i = MinRow();i<=MaxRow();i++) {
    for (int j = MinCol();j<=A.MaxCol();j++) {
      gRational x = A(i,j)*(gRational)totdenom;
      if(x.denominator() != 1) throw BadDenom();
      Tabdat(i,j) = x.numerator();
    }
    for (int j = A.MaxCol()+1;j<=MaxCol();j++)
      Tabdat(artificial[j],j) = totdenom;
  }
  for (int i = b.First();i<=b.Last();i++) 
    solution[i] = (gRational)Coeff[i];
}


Tableau<gRational>::Tableau(const Tableau<gRational> &orig) 
  : TableauInterface<gRational>(orig), Tabdat(orig.Tabdat), Coeff(orig.Coeff), 
    totdenom(orig.totdenom), denom(orig.denom), 
    tmpcol(orig.tmpcol), nonbasic(orig.nonbasic)
{ }

Tableau<gRational>::~Tableau()
{ }

Tableau<gRational>& Tableau<gRational>::operator=(const Tableau<gRational> &orig)
{
  TableauInterface<gRational>::operator=(orig);
  if(this!= &orig) {
    Tabdat = orig.Tabdat;
    Coeff = orig.Coeff;
    totdenom = orig.totdenom;
    denom = orig.denom;
    tmpcol = orig.tmpcol;
    nonbasic = orig.nonbasic;
  }