//
// $Source: /cvsroot/gambit/gambit/sources/tools/lp/ludecomp.imp,v $
// $Date: 2006/01/07 06:37:34 $
// $Revision: 1.6 $
//
// DESCRIPTION:
// Implementation of LU decomposition
//
// 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 "ludecomp.h"
#include "tableau.h"

// ---------------------------------------------------------------------------
// Class EtaMatrix
// ---------------------------------------------------------------------------

template <class T>
bool EtaMatrix<T>::operator==(const EtaMatrix<T> &a) const
{
  return ( col == a.col && etadata == a.etadata );
}

template <class T>
bool EtaMatrix<T>::operator!=(const EtaMatrix<T> &a) const
{
  return ( col != a.col || etadata != a.etadata );
}

// ---------------------------------------------------------------------------
// Class LUdecomp
// ---------------------------------------------------------------------------
// -------------------------
//  C-tors, D-tor, Operators
// -------------------------

// copy constructor

template <class T> 
LUdecomp<T>::LUdecomp( const LUdecomp<T> &a, Tableau<T> &t)
: tab(t), basis(t.GetBasis()), 
  scratch1(basis.First(), basis.Last()), 
  scratch2(basis.First(), basis.Last()),
  refactor_number( a.refactor_number ), iterations(a.iterations),
  total_operations( a.total_operations),
  parent(&a), copycount(0)

{ 
  ((LUdecomp<T> &)*parent).copycount++;
}

// Decomposes given matrix

template <class T> 
LUdecomp<T>::LUdecomp( Tableau<T> &t,
		      int rfac/* = 0 */)	
: tab(t), basis(t.GetBasis()),  
  scratch1(basis.First(), basis.Last()), 
  scratch2(basis.First(), basis.Last()),
  refactor_number(rfac), iterations(0), parent(NULL), copycount(0)
{
  int m = basis.Last() - basis.First() +1;
  total_operations = (m - 1) * m * (2 * m - 1) / 6;
}

// Destructor
template <class T> LUdecomp<T>::~LUdecomp() 
{ 
  if ( parent != NULL )
    ((LUdecomp<T> &) *parent).copycount--;
  if(copycount != 0) throw BadCount();
}



// -------------------------
//  Public Members
// -------------------------

// use this to copy ludecomps.... 
template <class T>
void LUdecomp<T>::Copy(const LUdecomp<T> &orig, Tableau<T> &t)
{
  if(this != &orig) {
    if (parent != NULL)
      ((LUdecomp<T> &) *parent).copycount--;
 
    tab = t;
    basis = t.GetBasis();
    
    L = Gambit::List<EtaMatrix<T> >();
    P = Gambit::List<int>();
    E = Gambit::List<EtaMatrix<T> >();
    U = Gambit::List<EtaMatrix<T> >();

    refactor_number = orig.refactor_number;
    iterations = orig.iterations;
    total_operations = orig.total_operations;
    parent = &orig;
    copycount = 0;
    ((LUdecomp<T> &)*parent).copycount++;
  }
}

template <class T>
void LUdecomp<T>::update( int col, int matcol )
{

  if( copycount != 0 ) throw BadCount();
  int m = basis.Last() - basis.First() + 1;

  iterations++;
  if (( refactor_number > 0 && iterations >= refactor_number ) ||
      ( refactor_number == 0 && RefactorCheck()) )  
    refactor();
  else {
    tab.GetColumn( matcol, scratch1); 
    solve( scratch1, scratch1 );
    if ( scratch1[col] == (T) 0 ) throw BadPivot();
    E.Append( EtaMatrix<T>( col, scratch1 ) );
    
    total_operations += iterations * m + 2 * m * m;    
  }
  
}


template <class T> 
void LUdecomp<T>::refactor( ) 
{

  L = Gambit::List<EtaMatrix<T> >();
  U = Gambit::List<EtaMatrix<T> >();
  E = Gambit::List<EtaMatrix<T> >();
  P = Gambit::List<int>();

  if ( !basis.IsIdent() ) FactorBasis();

  iterations = 0;
  int m = basis.Last() - basis.First() + 1;
  total_operations = (m - 1) * m * (2 * m - 1) / 6;
  if (parent != NULL) ((LUdecomp<T> &)*parent).copycount--;
  parent = NULL;
  
}

template <class T>
void LUdecomp<T>::solveT( const Gambit::Vector<T> &c, Gambit::Vector<T> &y ) const  
{

  if ( c.First() != y.First() || c.Last() != y.Last() ) throw Gambit::DimensionException();
  if ( c.First() != basis.First() || c.Last() != basis.Last()) throw Gambit::DimensionException();

  y = c;
  if ( basis.IsIdent() != true ) {
    BTransE( y );
    if ( parent != NULL ) 
      (*parent).solveT( y, y );
    else {
      FTransU( y );
      yLP_Trans( y );
    }
  }
}

template <class T>
void LUdecomp<T>::solve( const Gambit::Vector<T> &a, Gambit::Vector<T> &d ) const
{

  if ( a.First() != d.First() || a.Last() != d.Last() ) throw Gambit::DimensionException();
  if ( a.First() != basis.First() || a.Last() != basis.Last()) throw Gambit::DimensionException();
  
  d = a;
  if ( basis.IsIdent() != true ) {
    if ( parent != NULL ) 
      (*parent).solve(a,d);
    else {
      LPd_Trans( d );
      BTransU( d );
    }
    FTransE( d );
  }
}

template<class T>
void LUdecomp<T>::SetRefactor( int a )
{
  refactor_number = a;
}

// -----------------
//  Private Members
// -----------------


template<class T>
void LUdecomp<T>::FactorBasis()
{
  int i, j, piv;
  T pivVal;

  Gambit::Matrix<T> B(basis.First(), basis.Last(), 
	       basis.First(), basis.Last());

  for( i = basis.First(); i <= basis.Last(); i++ ) {
    tab.GetColumn( basis.Label(i), scratch2 );
    basis.CheckBasis();
    B.SetColumn( i, scratch2 );
  }
  
  for ( i = B.MinRow(); i <= B.MaxRow(); i++) {    
    pivVal = Gambit::abs ( B(i, i));
    piv = i;
    for ( j = i+1; j <= B.MaxRow(); j++) {
      if (  B( j, i ) * B( j, i )  > pivVal * pivVal ) {
	piv = j;
	pivVal = B( j, i );
      }
    }
    P.Append(piv);
    B.SwitchRows(i,piv);
    
    scratch2 = (T) 0;
    scratch2[i] = (T) 1 / B( i, i );
    for ( j = i+1; j <= B.MaxRow(); j++ ) {
      scratch2[j] =  - B(j, i) / B(i,i);
    }
    L.Append( EtaMatrix<T>(i, scratch2) );
    GaussElem(B, i, i);

  }
  for ( j = B.MinCol(); j <= B.MaxCol(); j++ ) {
    B.GetColumn( j, scratch2 );
    U.Append( EtaMatrix<T>( j, scratch2 ));
  }
}

template<class T>
void LUdecomp<T>::GaussElem(Gambit::Matrix<T> &B, int row, int col)
{
  if( B(row, col) == (T) 0) throw BadPivot();

  int i,j;

  for ( j = col+1; j <= B.MaxCol(); j++)
    B( row, j ) = B( row, j ) / B( row, col );

  for ( i = row+1; i <= B.MaxRow(); i++ ) 
    for ( j = col+1; j <= B.MaxCol(); j++ ) {
      B( i, j ) = B( i, j ) - ( B( i, col ) * B( row, j ) );
    }

  for ( i = row+1; i <= B.MaxRow(); i++ )
    B( i , col ) = 0;

  B( row, col ) = (T) 1;

}


template<class T>
void LUdecomp<T>::BTransE( Gambit::Vector<T> &y ) const
{

  int i;
  for ( i = E.Length(); i >= 1; i-- ) {
    ((LUdecomp<T> &) *this).scratch2 = y;
    VectorEtaSolve(scratch2, E[i], y );
  }
}
  
template<class T>
void LUdecomp<T>::FTransU( Gambit::Vector<T> &y ) const
{

  int i;
  for ( i = 1; i <= U.Length(); i++ ) {
    ((LUdecomp<T> &) *this).scratch2 = y;
    VectorEtaSolve(scratch2, U[i], y );
  }
}

template<class T>
void LUdecomp<T>::VectorEtaSolve( const Gambit::Vector<T> &v,  
				 const EtaMatrix<T>  &eta, 
				 Gambit::Vector<T> &y ) const
{

  int i, j;

  if ( v.First() != y.First() || v.Last() != y.Last() ) throw Gambit::DimensionException();

  for ( i = v.First(); i <= v.Last(); i++ ) {
    y[i] = v[i];
    if ( i == eta.col ) {
      for ( j = v.First(); j <= v.Last(); j++ )
	if ( j != eta.col ) y[i] -= v[j] * eta.etadata[j];
      y[i] /= eta.etadata[i];
    }    
  }
}

template<class T>
void LUdecomp<T>::FTransE( Gambit::Vector<T> &y ) const
{

  int i;
  for ( i = 1; i <= E.Length(); i++ ) {
    ((LUdecomp<T> &) *this).scratch2 = y;
    EtaVectorSolve(scratch2, E[i], y );
  }
}
  
template<class T>
void LUdecomp<T>::BTransU( Gambit::Vector<T> &y ) const
{

  int i;
  for ( i = U.Length(); i >= 1; i-- ) {
    ((LUdecomp<T> &) *this).scratch2 = y;
    EtaVectorSolve(scratch2, U[i], y );
  }
}

template<class T>
void LUdecomp<T>::EtaVectorSolve( const Gambit::Vector<T> &v,  
				 const EtaMatrix<T>  &eta, 
				 Gambit::Vector<T> &d ) const
{
  int i;
  T temp;

  if ( v.First() != d.First() || v.Last() != d.Last() ) throw Gambit::DimensionException();
  if ( eta.etadata[eta.col] == (T)0 )
    throw BadPivot(); // or we would have a singular matrix
  
  temp = v[eta.col] / eta.etadata[eta.col];

  for ( i = v.First(); i <= v.Last(); i++) {
    if ( i == eta.col ) d[i] = temp;
    else {
      d[i] = v[i] - temp * eta.etadata[i];
    }
  }
}

template<class T>
void LUdecomp<T>::yLP_Trans( Gambit::Vector<T> &y ) const
{
  int j;
  
  for (j = L.Length(); j >= 1; j--) {
    yLP_mult( y, j, ((LUdecomp<T> &) *this).scratch2 );
    y = scratch2;
  }
}


template<class T>
void LUdecomp<T>::yLP_mult( const Gambit::Vector<T> &y, int j, Gambit::Vector<T> &ans ) const
{
    
  if ( ans.First() != y.First() || ans.Last() != y.Last() ) throw Gambit::DimensionException();
  T temp;
  int i, k, l;
  
  l = j + y.First() - 1;

  for (i = y.First(); i <= y.Last(); i++) {
    if ( i != L[j].col) ans[i] = y[i];
    else {
      for ( k = ans.First(), temp = (T) 0; k <= ans.Last(); k++) {
	temp += y[k] * L[j].etadata[k];
      }
      ans[i] = temp;
    }
  }

  temp = ans[l];
  ans[l] = ans[P[j]];
  ans[P[j]] = temp;

}

template<class T>
void LUdecomp<T>::LPd_Trans( Gambit::Vector<T> &d ) const
{
  int j;
  for (j = 1; j <= L.Length(); j++) {
    LPd_mult( d, j, ((LUdecomp<T> &) *this).scratch2 );
    d = scratch2;
  }
}

template<class T>
void LUdecomp<T>::LPd_mult( Gambit::Vector<T> &d, int j, Gambit::Vector<T> &ans ) const
{
  
  if ( d.First() != ans.First() || d.Last() != ans.Last() ) throw Gambit::DimensionException();

  T temp;

  int i, k;

  k = j + d.First() - 1;
  temp = d[k];
  d[k] = d[P[j]];
  d[P[j]] = temp;

  for (i = d.First(); i <= d.Last(); i++) {
    if ( i == L[j].col ) ans[i] = d[i] * L[j].etadata[i];
    else {
      ans[i] = d[i] + d[ L[j].col ] * L[j].etadata[i];
    }
  }

  d[P[j]] = d[k];  
  d[k] = temp;

  
}

template<class T>
bool LUdecomp<T>::CheckBasis() 
{
  int i;
  bool ret = true;

  for (i = basis.First(); i <= basis.Last() && ret != false; i++)
    ret = ret && ( basis.Label(i) == -i );

  return ret;
}

template<class T>
bool LUdecomp<T>::RefactorCheck()
{
  int m = basis.Last() - basis.First() + 1;
  int i = iterations * (iterations * m + 2 * m * m );
  int k = total_operations + iterations * m + 2 * m * m;
  bool tmp;

  tmp = ( i > k );
  return tmp;
}
  
template <class T> LUdecomp<T>::BadCount::~BadCount()
{ }

template <class T> std::string LUdecomp<T>::BadCount::GetDescription(void) const
{
  return "Bad Reference count in LUdecomp";
}

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

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