#ifndef _MAdoubleRIX_H
#define _MAdoubleRIX_H

#include <iostream>
#include <string>
#include <cmath>

namespace math_tools {

using namespace std;

struct Exception { 
   const char * message;

   Exception() { message = ""; }
   Exception(const char * msg) { message = msg; }
};

struct IOException : Exception { };

typedef string String;

struct PrintWriter {
   void println(String arg) { cout << arg << endl; } 

   void flush() { cout << ::flush; }
};

typedef bool boolean;

/**
 * A matrix of floating point numbers with a set of 
 * linear algebra methods.
 * 
 * @author Matthew W. Coan
 * Started On: 10/01/2001
 * Last Modified: 11/28/2001
 */
class Matrix {
    double ** _mat;
    int _rows;
    int _cols;

public:

   double ** alloc(int rows, int cols);

   Matrix & round();

   Matrix();

   void freeMatrix();

   void freeMatrix2();

   Matrix(const Matrix & cp);

   /**
    * Construct a matrix from a two dimensional array
    * of double.
    * @param array the array of numbers.
    * @param rows the number of rows of the array.
    * @param cols the number of columns of the array.
    */
   Matrix(double ** array, int rows, int cols);

   Matrix & operator=(const Matrix & mat);

   /**
    * Construct a new matrix.
    * @param rows the number of rows in the new matrix.
    * @param cols the number of columns in the new matrix.
    * @param startValue the starting value for each entry in
    *                   the new matrix.
    */
   Matrix(int rows, int cols, double startValue);

   /**
    * Construct a new matrix.
    * @param rows the number of rows in the new matrix.
    * @param cols the number of colums in the new matrix.
    */
   Matrix(int rows, int cols);

   ~Matrix();
      
   /**
    * Return the number of rows.
    * @return the number of rows.
    */
   int rows();
   
   /**
    * Return the number of colums.
    * @return the number of colums.
    */
   int cols();
   
   /**
    * Set the value of an entry in the matrix.
    * @param i which row.
    * @param j which colum.
    * @param value the value to set.
    */
   void set(int i, int j, double value);
   
   /**
    * Return the value at an i,j entry in the matrix.
    * @param i the row.
    * @param j the column.
    * @return the value at that i,j entry.
    */
   double get(int i, int j) const;
   
   /**
    * doubleest to see if another matrix is equal 
    * to this matrix.
    * @param A the other matrix.
    * @return true if they are equal. False, otherwise.
    */
   boolean equal(Matrix & A);
   
   /**
    * Multiply this matrix by another and return the product 
    * as a new matrix.
    * @param A a matrix.
    * @return 0 if the dimensions are wrong.  Otherwise, a 
    *              new matrix that is the product of this 
    *              matrix and A is returned.
    */
   Matrix mul(const Matrix & A);
   
   /**
    * Add a matrix to this matrix and return the result.
    * @param A the matrix to add with.
    * @return 0 if the two matrixes can not be added 
    *         together.  Otherwise, the result of the addition
    *         operation is returned.
    */
   Matrix add(const Matrix & A);

   /**
    * Subtract a matrix to this matrix and return the result.
    * @param A the matrix to subtract.
    * @return 0 if the two matrixes can not be subtracted.
    *         Otherwise, the result of the subtraction
    *         operation is returned.
    */
   Matrix sub(const Matrix & A);
   
   /**
    * Scaler multipy.
    * @param k a scaler value to multiply the matrix by.
    * @return return the result matrix of the scaler 
    *         multiplication.
    */
   Matrix scalerMul(double k);
   
   /**
    * Compute the trace of a matrix.
    * @return the value of this matrixes trace.
    * @exception an Exception is raised if the 
    *            matrix is not square.
    */
   double trace();
   
   /**
    * Compute the transpose of this matrix.
    * @return the transpose of this matrix.
    */
   Matrix transpose();

private:
   double _det2x2();
   
public:
   /**
    * Remove a row from the matrix.
    * @param k the index of a row to remove.
    */
   void removeRow(int k);
   
   /**
    * Remove a column from the matrix.
    * @param k the index of a column to remove.
    */
   void removeCol(int k);
   
   /**
    * Compute the determanent of this matrix.
    * @return the determanent of this matirx.
    */
   double det();
   
   /**
    * Compute the adjoint of this matrix.
    * @return the adjoint of this matrix.
    */
   Matrix adj();
   
   /**
    * Compute the inverse of a matrix.
    * @return the inverse of this matrix or 0 if there is no inverse
    *         for this matrix.
    */
   Matrix inverse();

   /**
    * Swap two rows of this matrix.
    * @param i1 the first row.
    * @param i2 the second row.
    */
   void swapRows(int i1, int i2);
   
   /**
    * Multiply one row by a scaler and add it to another row.
    * @param srcRow the row to multiply by the scaler.
    * @param destRow the row to add the multipl too.
    */
   void addMulOfRowdoubleoRow(int srcRow, int destRow, double scaler);
   
   /**
    * Multiply a row by a scaler.
    * @param i the row index.
    * @param s the scaler.
    */
   void mulRowByScaler(int i, double s);
   
   /**
    * Row echelon reduce a matrix.
    * @return the row reduced matrix for this matrix.
    */
   Matrix rref();
      
   /**
    * Return the identity matrix for a given dimension.
    * @param dim the dimension for identity matrix.
    */
   static Matrix identity(int dim);
};

}

#endif
