//
// $Source: /cvsroot/gambit/gambit/sources/libgambit/matrix.imp,v $
// $Date: 2006/07/20 16:23:14 $
// $Revision: 1.6 $
//
// DESCRIPTION:
// Implementation of gbtMatrix method functions
//
// 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 "matrix.h"

namespace Gambit {

//-------------------------------------------------------------------------
//       Matrix<T>: Constructors, destructors, constructive operators
//-------------------------------------------------------------------------

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

template <class T> Matrix<T>::Matrix(unsigned int rows, unsigned int cols)
  : RectArray<T>(rows, cols)
{ }

template <class T> Matrix<T>::Matrix(unsigned int rows, unsigned int cols,
				     int minrows)
  : RectArray<T>(minrows, minrows + rows - 1, 1, cols)
{ }

template <class T> Matrix<T>::Matrix(int rl, int rh, int cl, int ch)
  : RectArray<T>(rl, rh, cl, ch)
{ }

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

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

template <class T> Matrix<T> &Matrix<T>::operator=(const T &c)
{
  for (int i = this->minrow; i <= this->maxrow; i++) 
    for (int j = this->mincol; j <= this->maxcol; j++) 
      (*this)(i, j) = c;
  return *this;
}

template <class T> Matrix<T> Matrix<T>::operator-(void)
{
  Matrix<T> tmp(this->minrow, this->maxrow, this->mincol, this->maxcol);
  for (int i = this->minrow; i <= this->maxrow; i++)
    for (int j = this->mincol; j <= this->maxcol; j++)
      tmp(i,j)= -(*this)(i,j);
  return tmp;
}


//-------------------------------------------------------------------------
//                     Matrix<T>: Additive operators
//-------------------------------------------------------------------------

template <class T> Matrix<T> Matrix<T>::operator+(const Matrix<T> &M) const
{
  if (!CheckBounds(M)) {
    throw DimensionException();
  }

  Matrix<T> tmp(this->minrow, this->maxrow, this->mincol, this->maxcol);
  for (int i = this->minrow; i <= this->maxrow; i++)  {
    T *src1 = this->data[i] + this->mincol;
    T *src2 = M.data[i] + this->mincol;
    T *dst = tmp.data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++)= *(src1++) + *(src2++);
    //assert((dst - 1) == tmp.data[i] + this->maxcol );
  }
  return tmp;
}

template <class T> Matrix<T> Matrix<T>::operator-(const Matrix<T> &M) const
{
  if (!CheckBounds(M)) {
    throw DimensionException();
  }

  Matrix<T> tmp(this->minrow, this->maxrow, this->mincol, this->maxcol);
  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T *src1 = this->data[i] + this->mincol;
    T *src2 = M.data[i] + this->mincol;
    T *dst = tmp.data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++)= *(src1++) - *(src2++);
    //assert((dst - 1) == tmp.data[i] + this->maxcol);
  }
  return tmp;
}

template <class T> Matrix<T> &Matrix<T>::operator+=(const Matrix<T> &M)
{
  if (!CheckBounds(M)) {
    throw DimensionException();
  }

  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T *src = M.data[i] + this->mincol;
    T *dst = this->data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) += *(src++);
  }
  return (*this);
}

template <class T> Matrix<T> &Matrix<T>::operator-=(const Matrix<T> &M)
{
  if (!CheckBounds(M)) {
    throw DimensionException();
  }

  for (int i = this->minrow; i <= this->maxrow; i++)  {
    T *src = M.data[i] + this->mincol;
    T *dst = this->data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) -= *(src++);
    //assert((dst - 1) == this->data[i] + this->maxcol); 
  }
  return (*this);
}

//-------------------------------------------------------------------------
//                  Matrix<T>: Multiplicative operators
//-------------------------------------------------------------------------

template <class T>
void Matrix<T>::CMultiply(const Vector<T> &in, Vector<T> &out) const
{
  if (!CheckRow(in) || !CheckColumn(out))  {
    throw DimensionException();
  }

  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T sum = (T)0;

    T *src1 = this->data[i] + this->mincol;
    T *src2 = in.data + this->mincol;
    int j = this->maxcol - this->mincol +1;
    while (j--)
      sum += *(src1++) * *(src2++);
    out[i] = sum;
  }
}

template <class T> Matrix<T> Matrix<T>::operator*(const Matrix<T> &M) const
{
  if (this->mincol != M.minrow || this->maxcol != M.maxrow) {
    throw DimensionException();
  }

  Matrix<T> tmp(this->minrow, this->maxrow, M.mincol, M.maxcol);
  Vector<T> column(M.minrow, M.maxrow);
  Vector<T> result(this->minrow, this->maxrow);
  for (int j = M.mincol; j <= M.maxcol; j++)  {
    M.GetColumn(j, column);
    CMultiply(column, result);
    tmp.SetColumn(j, result);
  }
  return tmp;
}

template <class T> Vector<T> Matrix<T>::operator*(const Vector<T> &v) const
{
  if (!CheckRow(v))  {
    throw DimensionException();
  }

  Vector<T> tmp(this->minrow, this->maxrow);
  CMultiply(v, tmp);
  return tmp;
}

template <class T>
void Matrix<T>::RMultiply(const Vector<T> &in, Vector<T> &out) const
{
  if (!CheckColumn(in) || !CheckRow(out)) {
    throw DimensionException();
  }

  out= (T)0;
  for (int i = this->minrow; i <= this->maxrow; i++)  {
    T k = in[i];

    T *src = this->data[i] + this->mincol;
    T *dst = out.data + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) += *(src++) * k;
    //assert(src - 1 == this->data[i] + this->maxcol);
  }
}

// transposed (row) vector*matrix multiplication operator
// a friend function of Matrix
template <class T>
Vector<T> operator*(const Vector<T> &v, const Matrix<T> &M)
{
  if (!M.CheckColumn(v)) {
    throw DimensionException();
  }
  Vector<T> tmp(M.MinCol(), M.MaxCol());
  M.RMultiply(v, tmp);
  return tmp;
}

template <class T> Matrix<T> Matrix<T>::operator*(const T &s) const
{
  Matrix<T> tmp(this->minrow, this->maxrow, this->mincol, this->maxcol);
  for (int i = this->minrow; i <= this->maxrow; i++)  {
    T *src = this->data[i] + this->mincol;
    T *dst = tmp.data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) = *(src++) * s;
    //assert((src - 1) == this->data[i] + this->maxcol);
  }
  return tmp;
}

// in-place multiplication by square matrix
template <class T> Matrix<T> &Matrix<T>::operator*=(const Matrix<T> &M)
{
  if (this->mincol != M.minrow || this->maxcol != M.maxrow) {
    throw DimensionException();
  }
  if (M.minrow != M.mincol || M.maxrow != M.maxcol) {
    throw DimensionException();
  }

  Vector<T> row(this->mincol, this->maxcol);
  Vector<T> result(this->mincol, this->maxcol);
  for (int i = this->minrow; i <= this->maxrow; i++)  {
    GetRow(i, row);
    M.RMultiply(row, result);
    SetRow(i, result);
  }
  return (*this);
}

template <class T> Matrix<T> &Matrix<T>::operator*=(const T &s)
{
  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T *dst = this->data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) *= s;
    //assert((dst - 1) == this->data[i] + this->maxcol);
  }
  return (*this);
}

template <class T> Matrix<T> Matrix<T>::operator/(const T &s) const
{
  if (s == (T) 0)   throw ZeroDivideException();

  Matrix<T> tmp(this->minrow, this->maxrow, this->mincol, this->maxcol);
  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T *src = this->data[i] + this->mincol;
    T *dst = tmp.data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) = *(src++) / s;
    //assert((src - 1) == this->data[i] + this->maxcol);
  }
  return tmp;
}

template <class T> Matrix<T> &Matrix<T>::operator/=(const T &s)
{
  if (s == (T) 0)   throw ZeroDivideException();

  for (int i = this->minrow; i <= this->maxrow; i++)  {
    T *dst = this->data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      *(dst++) /= s;
    //assert((dst - 1) == this->data[i] + this->maxcol);
  }
  return (*this);
}

//-------------------------------------------------------------------------
//                  Matrix<T>: Kronecker Product
//-------------------------------------------------------------------------

template <class T> Matrix<T> Matrix<T>::operator&(const Matrix<T> &M) const
{
  Matrix<T> tmp(1, (this->maxrow - this->minrow + 1)*(M.maxrow - M.minrow + 1), 
		1, (this->maxcol - this->mincol + 1)*(M.maxcol - M.mincol + 1));
 
  for (int i = 0; i <= this->maxrow - this->minrow; i++)
    for (int j = 1; j <= M.maxrow - M.minrow + 1; j++)
      for (int k = 0; k <= this->maxcol - this->mincol; k++)
	for (int l = 1; l <= M.maxcol - M.mincol + 1; l++)

	  tmp((M.maxrow - M.minrow + 1)*i + j, (M.maxcol - M.mincol + 1)*k + l) = 
	    (*this)(i+this->minrow, k+this->mincol) * M(j+M.minrow-1,l+M.mincol-1);

  return tmp;
}

//-------------------------------------------------------------------------
//                         Matrix<T>: Transpose
//-------------------------------------------------------------------------

template <class T> Matrix<T> Matrix<T>::Transpose() const
{
  Matrix<T> tmp(this->mincol, this->maxcol, this->minrow, this->maxrow);
 
  for (int i = this->minrow; i <= this->maxrow; i++)
    for (int j = this->mincol; j <= this->maxcol; j++)
      tmp(j,i) = (*this)(i,j);

  return tmp;
}

//-------------------------------------------------------------------------
//                    Matrix<T>: Comparison operators
//-------------------------------------------------------------------------

template <class T> bool Matrix<T>::operator==(const Matrix<T> &M) const
{
  if (!CheckBounds(M)) {
    throw DimensionException();
  }

  for (int i = this->minrow; i <= this->maxrow; i++)   {
    // inner loop
    T *src1 = M.data[i] + this->mincol;
    T *src2 = this->data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--)
      if(*(src1++) != *(src2++))   return false;
    //assert(src1 - 1 == M.data[i] + this->maxcol);
  }
  return true;
}

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

template <class T> bool Matrix<T>::operator==(const T &s) const
{
  for (int i = this->minrow; i <= this->maxrow; i++)   {
    T *src = this->data[i] + this->mincol;
    int j = this->maxcol - this->mincol + 1;
    while (j--) 
      if (*(src++) != s)   return false;
    //assert(src - 1 == this->data[i] + this->maxcol);
  }
  return true;
}

template <class T> bool Matrix<T>::operator!=(const T &s) const
{ return !(*this == s);  }

// Information

template <class T> Vector<T> Matrix<T>::Row(int i) const
{
  Vector<T> answer(this->mincol, this->maxcol);
  for (int j = this->mincol; j <= this->maxcol; j++)
    answer[j] = (*this)(i,j);
  return answer;
}

template <class T> Vector<T> Matrix<T>::Column(int j) const
{
  Vector<T> answer(this->minrow, this->maxrow);
  for (int i = this->minrow; i <= this->maxrow; i++)
    answer[i] = (*this)(i,j);
  return answer;
}


// more complex functions

template <class T> void Matrix<T>::MakeIdent(void)
{
  for (int i = this->minrow; i <= this->maxrow; i++) 
    for (int j = this->mincol; j <= this->maxcol; j++) {
      if (i == j) (*this)(i,j) = (T) 1;
      else (*this)(i,j)=(T) 0;
    }
}

template <class T> void Matrix<T>::Pivot(int row, int col)
{
  if (!this->CheckRow(row) || !this->CheckColumn(col)) {
    throw IndexException();
  }
  if (this->data[row][col] == (T) 0)  throw ZeroDivideException();

  T mult= (T)1/this->data[row][col];
  for(int j=this->mincol; j<=this->maxcol; j++)
    this->data[row][j]*= mult;
  for(int i=this->minrow; i<=this->maxrow; i++)
    {
      if(i!=row)
	{
	  mult= this->data[i][col];

	  // inner loop
	  T* src= this->data[row] + this->mincol;
	  T* dst= this->data[i] + this->mincol;
	  int j= this->maxcol-this->mincol+1;
	  while( j-- )
	    *(dst++)-= *(src++) * mult;
	  //assert( dst-1 == this->data[i] + this->maxcol ); // debug
	  // end inner loop
	}
    }
}

} // end namespace Gambit