//
// $Source: /cvsroot/gambit/gambit/sources/libgambit/sqmatrix.imp,v $
// $Date: 2006/07/20 16:23:14 $
// $Revision: 1.4 $
//
// DESCRIPTION:
// Implementation of square matrices
//
// 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 "sqmatrix.h"

namespace Gambit {

//-----------------------------------------------------------------------------
//     SquareMatrix<T>: Constructors, destructor, constructive operators
//-----------------------------------------------------------------------------

template <class T> SquareMatrix<T>::SquareMatrix(void)   { }

template <class T> SquareMatrix<T>::SquareMatrix(int size)
  : Matrix<T>(size, size)
{ }

template <class T> SquareMatrix<T>::SquareMatrix(const Matrix<T> &M)
  : Matrix<T>(M)
{ }

template <class T> SquareMatrix<T>::SquareMatrix(const SquareMatrix<T> &M)
  : Matrix<T>(M)
{ }

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

template <class T>
SquareMatrix<T> &SquareMatrix<T>::operator=(const SquareMatrix<T> &M)
{
  Matrix<T>::operator=(M);
  return *this;
}

//-----------------------------------------------------------------------------
//                SquareMatrix<T>: Public member functions
//-----------------------------------------------------------------------------

template <class T> SquareMatrix<T> SquareMatrix<T>::Inverse(void) const
{
  if (this->mincol != this->minrow || this->maxcol != this->maxrow) {
    throw DimensionException();
  }

  SquareMatrix<T> copy(*this);
  SquareMatrix<T> inv(this->maxrow);

  // initialize inverse matrix and prescale row vectors
  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T max= (T) 0;
    for (int j = this->mincol; j <= this->maxcol; j++)  {
      T abs = copy.data[i][j];
      if (abs < (T) 0)
	abs = -abs;
      if (abs > max)
	max = abs;
    }

    if (max == (T) 0) {
      throw SingularMatrixException();
    }

    T scale = (T) 1 / max;
    for (int j = this->mincol; j <= this->maxcol; j++)  {
      copy.data[i][j] *= scale;
      if (i == j)
	inv.data[i][j] = scale;
      else
	inv.data[i][j] = (T) 0;
    }
  }

  for (int i = this->mincol; i <= this->maxcol; i++)  {
    // find pivot row
    T max = copy.data[i][i];
    if (max < (T) 0)
      max = -max;
    int row = i;
    for (int j = i + 1; j <= this->maxrow; j++)  {
      T abs = copy.data[j][i];
      if (abs < (T) 0)
	abs = -abs;
      if (abs > max)  {
	max = abs;
	row = j;
      }
    }
    
    if (max <= (T) 0) {
      throw SingularMatrixException();
    }

    copy.SwitchRows(i, row);
    inv.SwitchRows(i, row);
    // scale pivot row
    T factor = (T) 1 / copy.data[i][i];
    for (int k = this->mincol; k <= this->maxcol; k++)  {
      copy.data[i][k] *= factor;
      inv.data[i][k] *= factor;
    }

    // reduce other rows
    for (int j = this->minrow; j <= this->maxrow; j++)  {
      if (j != i)  {
	T mult = copy.data[j][i];
	for (int k = this->mincol; k <= this->maxcol; k++)  {
	  copy.data[j][k] -= copy.data[i][k] * mult;
	  inv.data[j][k] -= inv.data[i][k] * mult;
	}
      }
    }
  } 

  return inv;
}

template <class T> T SquareMatrix<T>::Determinant(void) const
{
  T factor = (T) 1;
  SquareMatrix<T> M(*this);
  
  for (int row = this->minrow; row <= this->maxrow; row++)  {

    // Experience (as of 3/22/99) suggests that, in the interest of
    // numerical stability, it might be best to do Gaussian
    // elimination with respect to the row (of those feasible)
    // whose entry has the largest absolute value.
    int swap_row = row;
    for (int i = row+1; i <= this->maxrow; i++) {
      if (abs(M.data[i][row]) > abs(M.data[swap_row][row]))
	swap_row = i;
    }

    if (swap_row != row)  {
      M.SwitchRows(row, swap_row);
      for (int j = this->mincol; j <= this->maxcol; j++)
	M.data[row][j] *= (T) -1;
    }

    if (M.data[row][row] == (T)0) return (T)0;

    // now do row operations to clear the row'th column
    // below the diagonal
    for (int row1 = row+1; row1 <= this->maxrow; row1++)
      {
	factor = -M.data[row1][row]/M.data[row][row];
	for (int i = this->mincol; i <= this->maxcol; i++)
	  M.data[row1][i] += M.data[row][i]*factor;
      }
  }

  // finally we multiply the diagonal elements
  T det = (T) 1;
  for (int row = this->minrow; row <= this->maxrow; row++)  {
    det *= M.data[row][row];
  }
  return det;
}


}  // end namespace Gambit