//
// Copyright (C) 1991 Texas Instruments Incorporated.
//
// Permission is granted to any individual or institution to use, copy, modify,
// and distribute this software, provided that this complete copyright and
// permission notice is maintained, intact, in all copies and supporting
// documentation.
//
// Texas Instruments Incorporated provides this software "as is" without
// express or implied warranty.
//
// Created: MBN 04/21/89 -- Initial design and implementation
// Updated: MBN 06/22/89 -- Removed non-destructive methods
// Updated: LGO 08/09/89 -- Inherit from Generic
// Updated: MBN 08/20/89 -- Changed template usage to reflect new syntax
// Updated: MBN 09/11/89 -- Added conditional exception handling and base class
// Updated: LGO 10/05/89 -- Don't re-allocate data in operator= when same size
// Updated: LGO 10/19/89 -- Add extra parameter to varargs constructor
// Updated: MBN 10/19/89 -- Added optional argument to set_compare method
// Updated: LGO 12/08/89 -- Allocate column data in one chunk
// Updated: LGO 12/08/89 -- Clean-up get and put, add const everywhere.
// Updated: LGO 12/19/89 -- Remove the map and reduce methods
// Updated: MBN 02/22/90 -- Changed size arguments from int to unsigned int
// Updated: MJF 06/30/90 -- Added base class name to constructor initializer
//
// The parameterized Matrix<Type>  class is publicly   derived from the  Matrix
// class and implements two dimensional arithmetic matrices of a user specified
// type.   This is accompilshed by using  the parameterized  type capability of
// C++.  The only constraint placed on the type  is  that it must  overload the
// following operators: +, -,  *,  and /. Thus, it will  be possible to have  a
// matrix of  type Complex.  The Matrix<Type> class  is static in size, that is
// once a  Matrix<Type> of  a particular  size has been   declared, there is no
// dynamic growth or resize method available.
//
// Each matrix contains  a private  data section  that has a  Type** slot  that
// points to the  physical memory allocated  for the two  dimensional array. In
// addition, two integers  specify   the number  of  rows  and columns  for the
// matrix.  These values  are provided in the  constructors. A single protected
// slot  contains a pointer  to a compare  function  to   be used  in  equality
// operations. The default function used is the built-in == operator.
//
// Four  different constructors are provided.  The  first constructor takes two
// integer arguments  specifying the  row  and column  size.   Enough memory is
// allocated to hold row*column elements  of type Type.  The second constructor
// takes the  same two  first arguments, but  also accepts  an additional third
// argument that is  a reference to  an  object of  the appropriate  type whose
// value is used as an initial fill value.  The third constructor is similar to
// the third, except that it accpets a variable number of initialization values
// for the Matrix.  If there are  fewer values than elements,  the rest are set
// to zero. Finally, the last constructor takes a single argument consisting of
// a reference to a Matrix and duplicates its size and element values.
//
// Methods   are  provided   for destructive   scalar   and Matrix    addition,
// multiplication, check for equality  and inequality, fill, reduce, and access
// and set individual elements.  Finally, both  the  input and output operators
// are overloaded to allow for fomatted input and output of matrix elements.

#ifndef MATRIXH					// If no Matrix class,
#define MATRIXH					// define it

#ifndef BASE_MATRIXH				// If no base class definition
#include <cool/Base_Matrix.h>			// include it
#endif	

#ifndef STDARGH
#if defined(DOS)
extern "C" {
#include <stdarg.h>				// for variable arglists
}
#else
#include <stdarg.h>				// for variable arglists
#endif
#define STDARGH
#endif

template <class Type> Matrix {
  typedef Boolean (*Type ## _Matrix_Compare) (const Type&, const Type&); 
}

template <class Type>
class Matrix<Type> : public Matrix {
private:
  Type** data;					// Pointer to the Matrix 

protected:
  static Type ## _Matrix_Compare compare_s;	// Pointer operator== function
  friend Boolean Type ## _Matrix_is_data_equal (const Type&, const Type&);

public:
  Matrix<Type> (unsigned int, unsigned int);	// Matrix m (r,c);
  Matrix<Type> (unsigned int, unsigned int, const Type&); // Matrix m(r,c,val);
  Matrix<Type> (unsigned int, unsigned int, int, Type v1, ...);// Option values
  Matrix<Type> (const Matrix<Type>&);		// m1 = m2;
  ~Matrix<Type>();				// Destructor
  
  inline void put (unsigned int, unsigned int, Type); // Assign value
  inline Type get (unsigned int, unsigned int) CONST; // Get value

  Matrix<Type>& operator= (const Type&);	 // Assignment: m = 2;
  Matrix<Type>& operator= (const Matrix<Type>&); // Assignment: m = n;

  inline friend ostream& operator<< (ostream&, const Matrix<Type>*);
  friend ostream& operator<< (ostream&, const Matrix<Type>&);

  Boolean operator== (const Matrix<Type>&) CONST; // Matrix equality test
  inline Boolean operator!= (const Matrix<Type>&) CONST;  // inequality test
  inline void set_compare (Type ## _Matrix_Compare = NULL);// Compare function

  Matrix<Type> operator+ (const Type&) CONST;	// Scalar addition to a martix
  Matrix<Type> operator+ (const Matrix<Type>&) CONST;	// Matrix addition
  Matrix<Type> operator* (const Type&) CONST;	        // Scalar multiply
  Matrix<Type> operator* (const Matrix<Type>&) CONST;	// Matrix Multiply

  Matrix<Type>& operator+= (const Type&);	  // Addition and assignment
  Matrix<Type>& operator+= (const Matrix<Type>&); // Addition and assignment
  Matrix<Type>& operator*= (const Type&);	  // Multiply and assignment
  Matrix<Type>& operator*= (const Matrix<Type>&); // Multiply and assignment

  void fill (const Type&);			// Set elements to value
};


// Matrix -- Constructor specifiying size of matrix
// Input:    Row, Column parameters
// Output:   None

template<class Type> 
Matrix<Type>::Matrix<Type> (unsigned int row, unsigned int col)
#ifdef __cplusplus
 : Matrix(row, col)
#else
 : (row, col)
#endif
{
  this->data = new Type*[row];			// Allocate the rows
  Type* columns = new Type[col*row];		// Allocate the columns
  for (int i = 0; i < row; i ++)		// For each row in the Matrix
    this->data[i] = &columns[i*col];		// Fill in address of column
  if (this->compare_s == NULL)			// If no compare function
    this->compare_s = &Type ## _Matrix_is_data_equal; // Default is_equal
}


// Matrix -- Constructor specifiying size of matrix and initial value
// Input:    Row, Column parameters and initial value
// Output:   None

template<class Type> 
Matrix<Type>::Matrix<Type> (unsigned int row, unsigned int col,
			    const Type& value)
#ifdef __cplusplus
 : Matrix(row, col)
#else
 : (row, col)
#endif
{
  this->data = new Type*[row];			// Allocate the rows
  Type* columns = new Type[col*row];		// Allocate the columns
  for (int i = 0; i < row; i ++) {		// For each row in the Matrix
    this->data[i] = &columns[i*col];		// Fill in address of column
    for (int j = 0; j < col; j++)		// For each element in column
      this->data[i][j] = value;			// Assign initial value
  }
  if (this->compare_s == NULL)			// If no compare function
    this->compare_s = &Type ## _Matrix_is_data_equal; // Default is_equal
}


// Matrix -- Constructor that specifies number of elements and also a 
//           variable number of references to initialization values for 
//           successive elements in the Matrix
// Input:    Row, column dimensions, number of initial values, initial values
// Output:   None

template<class Type> 
Matrix<Type>::Matrix<Type> (unsigned int row, unsigned int col, int n,
			    Type v1, ...)
#ifdef __cplusplus
 : Matrix(row, col)
#else
 : (row, col)
#endif
{
  IGNORE(v1); // v1 added to disambiguate from Matrix(int, int, int&)
  this->data = new Type*[row];			// Allocate the rows
  Type* columns = new Type[col*row];		// Allocate the columns
  for (int i = 0; i < row; i ++)		// For each row in the Matrix
    this->data[i] = &columns[i*col];		// Fill in address of column
  if (n > 0) {					// If user specified values
    va_list argp;				// Declare argument list
    va_start (argp, row);			// Initialize macro
    Type* temp = (Type*)va_arg (argp, int);	// Skip over "col" argument
    temp = (Type*)va_arg (argp, int);		// Skip over "n" argument
    for (i = 0; i < row && n; i++)		// For remaining values given
      for (int j = 0; j < col && n; j++, n--)	// Moving sequentially in Matrix
	this->data[i][j] = (Type&) va_arg (argp, Type); // Extract and assign
    va_end(argp);
  }
  if (this->compare_s == NULL)			// If no compare function
    this->compare_s = &Type ## _Matrix_is_data_equal; // Default is_equal
}


// Matrix -- Constructor for reference to another Matrix object
// Input:    Matrix reference
// Output:   None

template<class Type> 
Matrix<Type>::Matrix<Type> (const Matrix<Type>& m)
#ifdef __cplusplus
 : Matrix(m)
#else
 : (m) 
#endif
{
  this->num_rows = m.num_rows;			// Save row specification
  this->num_cols = m.num_cols;			// Save column specification
  this->data = new Type*[this->num_rows];	// Allocate the rows
  Type* columns = new Type[num_cols*num_rows];	// Allocate the columns
  for (int i = 0; i < num_rows; i ++) {		// For each row in the Matrix
    this->data[i] = &columns[i*num_cols];	// Fill in address of column
    for (int j = 0; j < this->num_cols; j++)	// For each element in column
      this->data[i][j] = m.data[i][j];		// Copy value
  }
  if (this->compare_s == NULL)			// If no compare function
    this->compare_s = &Type ## _Matrix_is_data_equal; // Default is_equal
}


// ~Matrix -- Destructor for Matrix class that frees up storage
// Input:     *this
// Output:    None

template<class Type> 
Matrix<Type>::~Matrix<Type>() {
  delete this->data[0];				// Free up the column memory
  delete this->data;				// Free up the row memory
}


// get -- Get the element at specified index and return value
// Input: this*, row, column
// Output: Element value

template<class Type> 
inline Type Matrix<Type>::get (unsigned int row, unsigned int column) CONST {
#if ERROR_CHECKING
  if (row >= this->num_rows)			// If invalid size specified
    this->get_row_error (#Type, row);		// Raise exception
  if (column >= this->num_cols)			// If invalid size specified
    this->get_col_error (#Type, column);	// Raise exception
#endif
  return this->data[row][column];
}


// put -- Put the element value at specified index
// Input: *this, row, column, value
// Output: Element value

template<class Type> 
inline void Matrix<Type>::put (unsigned int row, unsigned int column, Type value) {
#if ERROR_CHECKING
  if (row >= this->num_rows)			// If invalid size specified
    this->put_row_error (#Type, row);		// Raise exception
  if (column >= this->num_cols)			// If invalid size specified
    this->put_col_error (#Type, column);	// Raise exception
#endif
  this->data[row][column] = value;		// Assign data value
}


// operator= -- Overload the assignment operator to assign a single 
//              value to the elements of a Matrix. 
// Input:       *this, reference to a value
// Output:      Reference to updated Matrix object

template<class Type> 
Matrix<Type>& Matrix<Type>::operator= (const Type& value) {
  for (int i = 0; i < this->num_rows; i++)	// For each row in Matrix
    for (int j = 0; j < this->num_cols; j++)	// For each column in Matrix
      this->data[i][j] = value;			// Assign value
  return *this;					// Return Matrix reference
}


// operator= -- Overload the assignment operator to copy the elements
//              in one Matrix to another. The existing storage for the 
//              destination matrix is freed up and new storage of the same 
//              size as the source is allocated.
// Input:       *this, reference to Matrix
// Output:      Reference to copied Matrix object

template<class Type> 
Matrix<Type>& Matrix<Type>::operator= (const Matrix<Type>& m) {
  if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) {
    delete this->data[0];			// Free up the column memory
    delete this->data;				// Free up the row memory
    this->num_rows = m.num_rows;		// Copy row specification
    this->num_cols = m.num_cols;		// Copy columns specification
    this->data = new Type*[this->num_rows];	// Allocate the rows
    Type* columns = new Type[num_cols*num_rows]; // Allocate the columns
    for (int i = 0; i < this->num_rows; i++)	// For each row
      this->data[i] = &columns[i*num_cols];	// Fill in address of column
  }
  for (int i = 0; i < this->num_rows; i++)	// For each row in the Matrix
    for (int j = 0; j < this->num_cols; j++)	// For each element in column
      this->data[i][j] = m.data[i][j];		// Copy value
  return *this;					// Return matrix reference
}


// operator<< -- Overload the output operator to print a matrix
// Input:        ostream reference, Matrix pointer
// Output:       ostream reference

template<class Type> Matrix {
inline ostream& operator<< (ostream& os, const Matrix<Type>* m) {
  return operator<< (os, *m);
}
}


// operator<< -- Overload the output operator to print a matrix
// Input:        ostream reference, Matrix reference
// Output:       ostream reference

template<class Type> Matrix {
ostream& operator<< (ostream& s, const Matrix<Type>& m) {
  for (int i = 0; i < m.rows(); i++) {	// For each row in matrix
    for (int j = 0; j < m.columns(); j++)	// For each column in matrix
      s << m.data[i][j] << " ";		// Output data element
    s << "\n";				// Output newline
  }
  return (s);					// Return ostream reference
}
}


// operator!= -- Perform not equal comparison test
// Input:        this*, matrix reference
// Output:       TRUE/FALSE

template<class Type> 
inline Boolean Matrix<Type>::operator!= (const Matrix<Type>& m) CONST {
  return (!operator== (m));
}


// operator== -- Compare the elements of two Matrices of Type Type using
//               the Compare pointer to funtion (default is ==). If one 
//               Matrix has more rows and/or columns than the other, the
//               result is FALSE
// Input:        Reference to Matrix of Type Type
// Output:       TRUE/FALSE

template<class Type> 
Boolean Matrix<Type>::operator== (const Matrix<Type>& m) CONST {
  if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
    return FALSE;				// Then not equal
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++)	// For each columne
    if ((*this->compare_s)(this->data[i][j],m.data[i][j]) == FALSE) // Same?
      return FALSE;				// Then no match
  return TRUE;					// Else same, so return TRUE
}


// is_data_equal -- Default data comparison function if user has not provided
//                  another one. Note that this is not inline because we need
//                  to take the address of it for the compare static variable
// Input:           Two Type references
// Output:          TRUE/FALSE

template<class Type> Matrix {
  Boolean Type ## _Matrix_is_data_equal (const Type& t1, const Type& t2) {
    return ((t1 == t2) ? TRUE : FALSE);
  }
}


// set_compare -- Specify the comparison function to be used
//                in logical tests of vector elements
// Input:         Pointer to a compare function
// Output:        None

template<class Type> 
inline void Matrix<Type>::set_compare (Type ## _Matrix_Compare c) {
  if (c == NULL)				// If no argument supplied
    this->compare_s = &Type ## _Matrix_is_data_equal; // Default is_equal
  else
    this->compare_s = c;			// Else set to user function
}


// operator+ -- Non-destructive matrix addition.  Note that the dimensions of
//              each matrix must be identical
// Input:       this*, matrix reference
// Output:      New matrix 

template<class Type> 
Matrix<Type> Matrix<Type>::operator+ (const Matrix<Type>& m) CONST {
  Matrix<Type> temp(*this);
  return temp += m;
}


// operator+ -- Non-destructive matrix addition of a scalar.
// Input:       this*, scalar value
// Output:      New matrix 

template<class Type> 
Matrix<Type> Matrix<Type>::operator+ (const Type& value) CONST {
  Matrix<Type> temp(this->num_rows, this->num_cols);
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++) 	// For each element in column
      temp.data[i][j] = (this->data[i][j] + value); // Add scalar
  return temp;
}


// operator* -- Non-destructive matrix multiply. Note that the dimensions of
//              each matrix must be identical
// Input:       this*, matrix reference
// Output:      New matrix 

template<class Type> 
Matrix<Type> Matrix<Type>::operator* (const Matrix<Type>& m) CONST {
  Matrix<Type> temp(this->num_rows, this->num_cols);
  if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
    this->op_mult_error (#Type, this->num_rows, this->num_cols, m.num_rows,
			 m.num_cols);
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++) {	// For each element in column
      temp.data[i][j] = (Type) 0.0;		// Initialize element value
      for (int k = 0; k < this->num_rows; k++)	// Loop over rows values
	temp.data[i][j] += (this->data[i][k] * m.data[k][j]); // Multiply
    }
  return temp;
}


// operator* -- Non-destructive matrix multiply by a scalar.
// Input:       this*, scalar value
// Output:      New matrix 

template<class Type> 
Matrix<Type> Matrix<Type>::operator* (const Type& value) CONST {
  Matrix<Type> temp(this->num_rows, this->num_cols);
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++) 	// For each element in column
      temp.data[i][j] = (this->data[i][j] * value); // Multiply
  return temp;
}


// operator+= -- Destructive matrix addition of a scalar.
// Input:        this*, scalar value
// Output:       New matrix reference

template<class Type> 
Matrix<Type>& Matrix<Type>::operator+= (const Type& value) {
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++) 	// For each element in column
      this->data[i][j] += value;		// Add scalar
  return *this;
}


// operator+= -- Destructive matrix addition with assignment. Note that the
//               dimensions of each matrix must be identical
// Input:        this*, matrix reference
// Output:       Updated this* matrix reference

template<class Type> 
Matrix<Type>& Matrix<Type>::operator+= (const Matrix<Type>& m) {
  if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
    this->op_plus_assign_error (#Type, this->num_rows, this->num_cols,
				m.num_rows, m.num_cols);
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++)	// For each element in column
      this->data[i][j] = this->data[i][j] + m.data[i][j]; // Add elements
  return *this;
}


// operator*= -- Destructive matrix multiplication by a scalar.
// Input:        this*, scalar value
// Output:       New matrix reference

template<class Type> 
Matrix<Type>& Matrix<Type>::operator*= (const Type& value) {
  for (int i = 0; i < this->num_rows; i++)	// For each row
    for (int j = 0; j < this->num_cols; j++) 	// For each element in column
      this->data[i][j] *= value;		// Multiply by scalar
  return *this;
}


// operator*= -- Destructive matrix multiply with assignment. Note that the
//               dimensions of each matrix must be identical
// Input:        this*, matrix reference
// Output:       Updated this* matrix reference

template<class Type> 
Matrix<Type>& Matrix<Type>::operator*= (const Matrix<Type>& m) {
  if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
    this->op_mult_assign_error (#Type, this->num_rows, this->num_cols,
				m.num_rows, m.num_cols);
  Type** temp = new Type*[this->num_rows];	// Temporary
  Type* columns = new Type[num_cols*num_rows];	// Allocate the columns
  for (int i = 0; i < this->num_rows; i++) {	// For each row
    temp[i] = &columns[i*num_cols];		// Fill in address of column
    for (int j = 0; j < this->num_cols; j++) {	// For each element in column
      temp[i][j] = (Type) 0.0;			// Initialize element value
      for (int k = 0; k < this->num_cols; k++)	// Loop over column values
	temp[i][j] += (this->data[i][k] * m.data[k][j]); // Multiply
    }
  }
  delete this->data[0];				// Free up the column memory
  delete this->data;				// Free up the row memory
  this->data = temp;				// Pointer to new storage
  return *this;					// Return matrix reference
}


// fill -- Set all elements of a matrix to a specified fill value
// Input:  this*, reference to fill value
// Output: None

template<class Type> 
void Matrix<Type>::fill (const Type& value) {
  for (int i = 0; i < this->num_rows; i++)	// For each row in the Matrix
    for (int j = 0; j < this->num_cols; j++)	// For each element in column
      this->data[i][j] = value;			// Assign fill value
}

#endif						// End of MATRIXH


