gsmatrix.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/math/gsmatrix.imp,v $
// $Date: 2002/08/26 05:50:04 $
// $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 "base/base.h"
#include "gsmatrix.h"
#include "gmath.h"   // for abs(double)

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

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

template <class T> gSquareMatrix<T>::gSquareMatrix(int size)
  : gMatrix<T>(size, size)
{ }

template <class T> gSquareMatrix<T>::gSquareMatrix(const gMatrix<T> &M)
  : gMatrix<T>(M)
{ assert ( M.NumRows() == M.NumColumns() ); }

template <class T> gSquareMatrix<T>::gSquareMatrix(const gSquareMatrix<T> &M)
  : gMatrix<T>(M)
{ }

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

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

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

template <class T> gSquareMatrix<T>::MatrixSingular::~MatrixSingular() { }

template <class T>
gText gSquareMatrix<T>::MatrixSingular::Description(void) const
{
  return "Attempted to invert a singular matrix";
}

template <class T> gSquareMatrix<T> gSquareMatrix<T>::Inverse(void) const
{
  if (mincol != minrow || maxcol != maxrow) {
    throw gRectArray<T>::BadDim();
  }

  gSquareMatrix<T> copy(*this);
  gSquareMatrix<T> inv(maxrow);

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

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

    T scale = (T) 1 / max;
    for (int j = mincol; j <= 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 = mincol; i <= 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 <= 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 MatrixSingular();
    }

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

    // reduce other rows
    for (int j = minrow; j <= maxrow; j++)  {
      if (j != i)  {
	T mult = copy.data[j][i];
	for (int k = mincol; k <= 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 gSquareMatrix<T>::Determinant(void) const
{
  T factor = (T) 1;
  gSquareMatrix<T> M(*this);
  
  for (int row = minrow; row <= 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 <= 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 = mincol; j <= 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 <= maxrow; row1++)
      {
	factor = -M.data[row1][row]/M.data[row][row];
	for (int i = mincol; i <= maxcol; i++)
	  M.data[row1][i] += M.data[row][i]*factor;
      }
  }

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

template <class T> gOutput &operator<<(gOutput &to, const gSquareMatrix<T> &M)
{
  to << (gMatrix<T>)M; return to;
}