//********************************************************************
//Matrix.h C++ template class
//Author: Jun Hong
//web:http://www2.uic.edu/~jhong4/sample.html
//Email:jhong4@uic.edu
//********************************************************************

//Note: This matrix.h file defines three classes including Matrix<F>,
//SubMatRef<F>, and const_SubMatRef<F>. The Matrix class defines
//most of the matrix operations as overloaded operators or methods.
//For the Matrix<F> class, two smart pointers, SubMatRef<F> and
//const_SubMatRef<F>, are defined. Users can used them to develop
//other matrix opertions. Iterators could be also added later by Users as
//nested class in Matrix<F>

#include <iostream.h>
#include <string.h>
#include <utility>
#include <iomanip.h>
#include <math.h>
#include <strstream.h>
#include <algorithm>

// declarations of classes and functions
template <class F> class Matrix;
template <class F> class const_SubMatRef;
template <class F> class SubMatRef;

template <class F> bool operator == (const Matrix<F>& m1, const
Matrix<F>& m2);
template <class F> bool operator != (const Matrix<F>& m1, const
Matrix<F>& m2);
template <class F> Matrix<F> operator * (const Matrix<F>& m, const
F& no);
template <class F> Matrix<F> operator * (const F& no, const
Matrix<F>& m);
template <class F> Matrix<F> operator * (const Matrix<F>& m1, const
Matrix<F>& m2);
template <class F> Matrix<F> operator * (int no, const Matrix<F>&
m);
template <class F> Matrix<F> operator * (const Matrix<F>& m, int
no);

//=================Define Template Class Matrix=========================

template <class F>
class Matrix {

   friend class SubMatRef<F>;
   friend class const_SubMatRef<F>;

public:

   // Constructors and destructor
   Matrix (const Matrix<F>& m);
   Matrix (size_t row , size_t col );
   Matrix (size_t n);
   Matrix (size_t row, size_t col, char* s);

   ~Matrix ();

   // Size, this function returns a pair s, s.first and s.second
   pair<size_t, size_t> size() const { return pair<size_t,size_t>
        (md->RowSiz, md->ColSiz);}

   //RowNo and ColNo
   size_t RowNo() const {return md->Row;}
   size_t ColNo() const {return md->Col;}

   // assignment operators overloaded =
   Matrix<F>& operator = (const Matrix<F>& m);

   // unary operators + and -
   Matrix<F> operator + () { return *this;}
   Matrix<F> operator - ();

   // subscript operators
   F& operator ()(size_t row, size_t col) const;
   F& operator ()(size_t row, size_t col);

   // transpose
   Matrix<F> transpose ();

   // Null
   void Null();

   // Combined assignment +=,-=, *=
   Matrix<F>& operator += (const Matrix<F>& m);
   Matrix<F>& operator -= (const Matrix<F>& m);
   Matrix<F>& operator *= (const Matrix<F>& m);

   // logical operators
   friend bool operator == <>(const Matrix<F>& m1, const Matrix<F>& m2);
   friend bool operator != <>(const Matrix<F>& m1, const Matrix<F>& m2);

   // binary calculation operators + - *
   Matrix<F> operator+ (const Matrix<F>&m2);
   Matrix<F> operator- (const Matrix<F>&m2);

   friend Matrix<F> operator * <>(const Matrix<F>& m1, const Matrix<F>&m2);
   friend Matrix<F> operator * <>(const Matrix<F>& m, const F& no);
   friend Matrix<F> operator * <>(const F& no, const Matrix<F>& m);
   friend Matrix<F> operator * <>(int no, const Matrix<F>& m);
   friend Matrix<F> operator * <>(const Matrix<F>& m, int no);

   // Miscellaneous methods
   F Det() const;
   void Unit ();

   // identity
   static Matrix<F> identity (const size_t& n);

   // inverse
   pair <F, Matrix<F> > inverse () const {
        return inverse( Field_traits<F>::arith_type() );
   }

   pair <F, Matrix<F> > inverse(arith_exact) const;
   pair <F, Matrix<F> > inverse(arith_approx) const;

   // rowBasis
   Matrix<F> rowBasis() const {
        return rowBasis( Field_traits<F>::arith_type() );
   }

   Matrix<F> rowBasis(arith_exact) const;
   Matrix<F> rowBasis(arith_approx) const;

   // conversion
   template <class G>
   Matrix<F>(const Matrix<G> & p);

   //a(r1, r2, c1, c2);
   SubMatRef<F> operator ()(size_t a, size_t b, size_t c, size_t d){
        return SubMatRef<F>(*this, a, b, c, d);
   }

   const_SubMatRef<F> operator ()(size_t a, size_t b, size_t c, size_t d)
   const{
        return const_SubMatRef<F>(*this, a, b, c, d);
   }

   //a[i]
   SubMatRef<F> operator[] (size_t i){
        return SubMatRef<F>(*this, i, i, 0, size().second-1);
   }

   const_SubMatRef<F> operator[] (size_t i) const{
        return const_SubMatRef<F>(*this, i, i, 0, size().second-1);
   }

   //Function to swap this matrix with another
   void swap (Matrix &a){
        if (size()!=a.size())
                throw "size error";
        for (size_t i=0; i<md->Row; i++)
                for (size_t j=0; j<md->Col; j++)
                std::swap(md->element[i][j], a(i,j));
   }

  //assign SubMatRef<F> to Matrix<F>
   Matrix<F>& operator= (const SubMatRef<F>& a){
        md =new membervariables (a.r2-a.r1+1, a.c2-a.c1+1, 0);
        for(size_t i=0, ii=a.r1; i<a.r2-a.r1+1; i++, ii++)
                for (size_t j=0, jj=a.c1; j<a.c2-a.c1+1; j++, jj++)
                        md->element[i][j]=a.rA(ii,jj);
        return *this;
  }

private:
    struct membervariables
    {
       F**element;
        size_t Row, Col, RowSiz, ColSiz;
        int count;

     membervariables  (size_t row, size_t col, F** v)
     {
         Row = row;
         RowSiz = row;
         Col = col;
         ColSiz = col;

         count = 1;
         element = new F* [row];
         size_t rowlen = col * sizeof(F);

         for (size_t i=0; i < row; i++){
              element[i] = new F[col];
              if (v) memcpy( element[i], v[i], rowlen);
         }
      }

      ~membervariables  (){
          for (size_t i=0; i < RowSiz; i++)
              delete [] element[i];
          delete [] element;
      }
  };

    membervariables  *md;

    void clone ();
    void realloc (size_t row, size_t col);
    int pivot (size_t row);
};

//=======================definition of template class const_SubMatRef====

   template <class F>
   class const_SubMatRef{
        friend class Matrix<F>;
        friend class SubMatRef<F>;
   protected:
        Matrix<F> &rA;
        size_t r1, r2, c1, c2;
        const_SubMatRef(const Matrix<F> &aa, size_t a, size_t b, size_t c,
                size_t d): rA(const_cast<Matrix<F> &>(aa)), r1(a), r2(b),
                c1(c), c2(d) {}

   public:

        //Conversion to Matrix<F>
        operator Matrix<F> ()  const{
                Matrix <F> a(r2-r1+1,c2-c1+1);
                for (size_t i=0, ii=r1; i<=r2-r1; i++, ii++)
                        for (size_t j=0, jj=c1; j<=c2-c1; j++, jj++)
                           a.md->element[i][j]=rA.md->element[ii][jj];
                return a;
        }

        //This function is used by ref-=mu*cref
        friend pair<F, const_SubMatRef<F> > operator* (F mu,
           const const_SubMatRef<F>& m){

           return pair<F, const_SubMatRef<F> >(mu, m);
        }

  };

//====================definition of template class SubMatRef===========

   template <class F>
   class SubMatRef: public const_SubMatRef<F> {
        friend class Matrix<F>;
   private:
        SubMatRef(const Matrix<F>& aa, size_t u, size_t v, size_t w,
                size_t x):
                const_SubMatRef<F> (aa, u, v, w, x){}
   public:

        //Conversion to Matrix<F>
        operator Matrix<F> ()  const{
                Matrix <F> a(r2-r1+1, c2-c1+1);
                for (size_t i=0, ii=r1; i<=r2-r1; i++, ii++)
                     for (size_t j=0, jj=c1; j<=c2-c1; j++, jj++)
                        a.md->element[i][j]=rA.md->element[ii][jj];
                return a;
        }

        SubMatRef& operator=(const Matrix<F> &a){
                if(a.size().first !=r2-r1+1 || a.size().second !=c2-c1+1 )
                        throw "Size error";

                for(size_t r=r1, i=0; r<=r2; r++, i++)
                        for (size_t c=c1, j=0; c<=c2; c++, j++)
                                rA(r,c)=a(i,j);
                return *this;
        }

        SubMatRef& operator=(const const_SubMatRef<F> &s){
                if ( r2-r1!=s.r2-s.r1 || c2-c1!=s.c2-s.c1)
                        throw "Size error";

                for(size_t r=r1, i=s.r1; r<=r2; ++r, ++i)
                    for (size_t c=c1, j=s.c1; c<=c2; ++c, ++j)
                        rA(r,c)=s.rA(i,j);
                return *this;
        }

        //Function to swap this SubMatRef with another
        void swap (SubMatRef<F> a){
                if (r2-r1!=a.r2-a.r1||c2-c1!=a.c2-a.c1)
                        throw "size error";
                for (size_t i=r1, ii=a.r1; i<=r2; i++, ii++)
                        for (size_t j=c1, jj=a.c1; j<=c2; j++, jj++)
                                std::swap(rA(i,j), a.rA(ii,jj));
        }

        //ref*=mu
        SubMatRef& operator*= (F mu){
                for (size_t i=r1; i<=r2; i++)
                        for (size_t j=c1; j<=c2; j++)
                                rA(i,j)*=mu;
                return *this;
        }

        //ref+ref
        Matrix<F> operator + (const SubMatRef &a){
                Matrix<F> A(a.r2-a.r1+1, a.c2-c1+1);
                A=a;
                Matrix<F> B(r2-r1+1, c2-c1+1);
                B=*this;
                Matrix<F> C(r2-r1+1, c2-c1+1);
                C=A+B;
                return C;
        }

        //ref-ref
        Matrix<F> operator - (const SubMatRef &a){
            Matrix<F> A(a.r2-a.r1+1, a.c2-c1+1);
            A=a;
            Matrix<F> B(r2-r1+1, c2-c1+1);
            B=*this;
            Matrix<F> C(r2-r1+1, c2-c1+1);
            C=B-A;
            return C;
       }

        //ref+Matrix
        Matrix<F> operator + (const Matrix<F>& a){
            Matrix<F> A(r2-r1+1, c2-c1+1);
            A=*this;
            Matrix<F> C(r2-r1+1, c2-c1+1);
            C=A+a;
            return C;
        }

        //ref-=mu*cref
        SubMatRef& operator -= ( pair<F, const_SubMatRef<F> > a){
                for (size_t i=r1, ii=a.second.r1; i<=r2; i++, ii++)
                  for (size_t j=c1, jj=a.second.c1; j<=c2; j++, jj++)
                    rA(i,j) -= (a.second.rA(ii,jj)*a.first);
                return *this;
        }

        SubMatRef& operator=( const SubMatRef <F> &s){
                for (size_t i=r1, ii=s.r1; i<=r2; i++, ii++)
                   for (size_t j=c1, jj=s.c1; j<=c2; j++, jj++)
                        rA(i,j)=s.rA(ii, jj);

                return *this;
        }
   };

//=======================Definiation of Class Member Functions==========

// Constructor
template <class F> inline
Matrix<F>::Matrix (size_t row, size_t col)
{
        md = new membervariables ( row, col, 0);
        Null();
}

// Constructor
template <class F> inline
Matrix<F>::Matrix (size_t n){
        md=new membervariables (n,n,0);
        Null();
}

// Copy constructor
template <class F> inline
Matrix<F>::Matrix (const Matrix<F>& m)
{
    md = m.md;
    md->count++;
}

// Constructor
template <class F> inline
Matrix<F>::Matrix (size_t row, size_t col, char* s)
{
        F* elt=new F[row*col];
        istrstream iss(s);
        for (size_t i=0; i<row*col; ++i)
        iss >> elt[i];

        md=new membervariables (row, col, 0);

        if (md->count>1) clone();
        size_t a=0;
        for (size_t k=0; k<md->Row; ++k)
                for (size_t j=0; j<md->Col; ++j)
                        (md->element)[k][j] = elt[a++];
}

// clone constructor
template <class F> inline void
Matrix<F>::clone ()
{
    md->count--;
    md = new membervariables ( md->Row, md->Col, md->element);
}

// Destructor
template <class F> inline
Matrix<F>::~Matrix ()
{
   if (--md->count == 0) delete md;
}

// assignment operator
template <class F> inline Matrix<F>&
Matrix<F>::operator = (const Matrix<F>& m)
{
    m.md->count++;
    if (--md->count == 0) delete md;
    md = m.md;

    return *this;
}

//  reallocation method
template <class F> inline void
Matrix<F>::realloc (size_t row, size_t col)
{
   if (row == md->RowSiz && col == md->ColSiz)
   {
      md->Row = md->RowSiz;
      md->Col = md->ColSiz;
      return;
   }

   membervariables  *m1 = new membervariables ( row, col, NULL);
   size_t colSize = (md->Col<col? md->Col:col) * sizeof(T);
   size_t minRow = (md->Row<row? md->Row:row);

   for (size_t i=0; i < minRow; i++)
      memcpy( m1->element[i], md->element[i], colSize);

   if (--md->count == 0)
       delete md;
   md = m1;

   return;
}

// a==b
template <class F> inline bool
operator == (const Matrix<F>&  m1, const Matrix<F>& m2)
{
   if (m1.RowNo() != m2.RowNo() || m1.ColNo() != m2.ColNo())
      return false;

   for (size_t i=0; i < m1.RowNo(); i++)
      for (size_t j=0; j < m1.ColNo(); j++)
         if (m1.md->element[i][j] != m2.md->element[i][j])
            return false;

   return true;
}

// a!=b
template <class F> inline bool
operator != (const Matrix<F>& m1, const Matrix<F>& m2)
{
    return (m1 == m2) ? false : true;
}

// a+=b
template <class F> inline Matrix<F>&
Matrix<F>::operator += (const Matrix<F>& m)
{
   if (md->Row != m.md->Row || md->Col != m.md->Col)
      cout<< "Inconsistent Matrix sizes in += !";
   if (md->count > 1) clone();
   for (size_t i=0; i < m.md->Row; i++)
      for (size_t j=0; j < m.md->Col; j++)
         md->element[i][j] += m.md->element[i][j];
   return *this;
}

// a-=b
template <class F> inline Matrix<F>&
Matrix<F>::operator -= (const Matrix<F>& m)
{
   if (md->Row != m.md->Row || md->Col != m.md->Col)
      cout<< "Inconsistent Matrix sizes in -= !";
   if (md->count > 1) clone();
   for (size_t i=0; i < m.md->Row; i++)
      for (size_t j=0; j < m.md->Col; j++)
         md->element[i][j] -= m.md->element[i][j];
   return *this;
}

// a *=b
template <class F> inline Matrix<F>&
Matrix<F>::operator *= (const Matrix<F>& m)
{
   if (md->Col != m.md->Row)
      cout<< "Inconsistent Matrix sizes in *= !";

   *this = *this * m;

   return *this;
}

// -a
template <class F> inline Matrix<F>
Matrix<F>::operator - ()
{
   Matrix<F> temp(md->Row,md->Col);

   for (size_t i=0; i < md->Row; i++)
      for (size_t j=0; j < md->Col; j++)
         temp.md->element[i][j] = - md->element[i][j];

   return temp;
}

// a+b
template <class F> inline Matrix<F> Matrix<F>::operator + (const Matrix<F>& m2)
{
   if (md->Row != m2.md->Row || md->Col != m2.md->Col)
   cout<< "Inconsistent Matrix size in addition!";

   Matrix<F> temp(md->Row, md->Col);

   for (size_t i=0; i < m2.md->Row; i++)
      for (size_t j=0; j < m2.md->Col; j++)
         temp.md->element[i][j] = md->element[i][j]+m2.md->element[i][j];

   return temp;
}

// a-b
template <class F>
inline Matrix<F> Matrix<F>::operator - (const Matrix<F>& m2){

   if (md->Row != m2.md->Row || md->Col != m2.md->Col)
     cout<< "Inconsistent Matrix size in subtraction!";

   Matrix<F> temp(md->Row, md->Col);

   for (size_t i=0; i < md->Row; i++)
      for (size_t j=0; j < md->Col; j++)
         temp.md->element[i][j] =md->element[i][j]-m2.md->element[i][j];

   return temp;
}

// mu*a int
template <class F> inline Matrix<F>
operator * (const Matrix<F>& m, int no)
{
   Matrix<F> temp(m.md->Row, m.md->Col);

   for (size_t i=0; i < m.md->Row; i++)
      for (size_t j=0; j < m.md->Col; j++)
         temp.md->element[i][j] = no * m.md->element[i][j];

   return temp;
}

// mu*a
template <class F> inline Matrix<F>
operator * (const Matrix<F>& m, const F& no)
{
   Matrix<F> temp(m.md->Row, m.md->Col);

   for (size_t i=0; i < m.md->Row; i++)
      for (size_t j=0; j < m.md->Col; j++)
         temp.md->element[i][j] = no * m.md->element[i][j];

   return temp;
}

// a*mu
template <class F> inline Matrix<F>
operator * (const F& no, const Matrix<F>& m)
{
   return (m * no);
}

// int a*mu
template <class F> inline Matrix<F>
operator * (int no, const Matrix<F>& m)
{
   return (m * no);
}

// a*b
template <class F> inline Matrix<F>
operator * (const Matrix<F>& m1, const Matrix<F>& m2)
{
   if (m1.md->Col != m2.md->Row)
      cout<<"Matrix<F>::operator*: Inconsistent Matrix size in multiplication!";

   Matrix<F> temp(m1.md->Row, m2.md->Col);

   for (size_t i=0; i < m1.md->Row; i++)
      for (size_t j=0; j < m2.md->Col; j++)
         {
         temp.md->element[i][j] = F(0);
         for (size_t k=0; k < m1.md->Col; k++)
            temp.md->element[i][j] += m1.md->element[i][k] * m2.md->element[k][j];
         }
   return temp;
}

// Transpose
template <class F> inline Matrix<F> Matrix<F>::transpose(){
   Matrix<F> temp (md->Col,md->Row);
   for (size_t i=0; i < md->Row; i++){
      for (size_t j=0; j <md->Col; j++)
                  temp.md->element[j][i] = md->element [i][j];}
   return temp;
}

//Null set all elements in a Matrix to F(0)
template <class F> inline void
Matrix<F>::Null() {

    if (md->count > 1) clone();
    for (size_t i=0; i < md->Row; i++)
        for (size_t j=0; j < md->Col; j++)
                md->element[i][j] = F(0);
    return;
}

// Identity
template <class F> inline Matrix<F> Matrix<F>::identity (const size_t& n){

        Matrix<F> temp(n);
    for (size_t i=0; i < temp.md->Row; i++)
        for (size_t j=0; j < temp.md->Col; j++)
            temp.md->element[i][j] = i == j ? F(1) : F(0);
        return temp;
}

//Unit
template <class F> inline void Matrix <F>::Unit ()
{
    if (md->count > 1) clone();
    size_t row = (md->Row<md->Col? md->Row:md->Col);
    md->Row = md->Col = row;

    for (size_t i=0; i < md->Row; i++)
        for (size_t j=0; j < md->Col; j++)
            md->element[i][j] = i == j ? F(1) : F(0);
    return;
}

//Subscript
template <class F> inline F& Matrix <F>::operator () (size_t row, size_t
col) {
   if (row >= md->Row || col >= md->Col)
        cout<< "Index out of range!";
   if (md->count>1) clone();
   return md->element[row][col];
}

//Subscript for const
template <class F> inline F& Matrix <F>:: operator ()(size_t row, size_t
col) const{
   if (row >= md->Row || col >= md->Col)
        cout<< "Index out of range!";
   return md->element[row][col];
}

//Converstion
template <class F>
template <class G> inline Matrix<F>::
Matrix<F>(const Matrix<G> & p) {
     md=new membervariables(p.RowNo(), p.ColNo(), 0);
     for (size_t i=0; i<p.RowNo(); i++){
        md->element[i]= new F[p.ColNo()];
             for (size_t j=0; j<p.ColNo(); j++)
                     md->element[i][j] = p(i,j);}
}

//Inverse
template <class F> inline pair<F, Matrix<F> >
Matrix<F>::inverse( arith_exact ) const {
   Matrix<F> A(md->Row, md->Col);
        for(size_t i=0; i<md->Row; i++)
           for(size_t j=0; j<md->Col; j++)
                A.md->element[i][j]=md->element[i][j];

   Matrix<F> B=Matrix<F>::identity (ColNo());
   F det(1);

   for (size_t r=0; r<A.ColNo(); ++r){
        if (A(r,r)==F(0)){
           for(size_t q=r+1; q<A.ColNo(); ++q)
                if(A(q,r)!=F(0)){
                   A[r].swap(A[q]);
                   B[r].swap(B[q]);
                   det=-det;
                }
           if(A(r,r)==F(0)) return pair<F, Matrix<F> >(0,*this);
        }

        F temp=A(r,r);
        A[r] *=(F(1)/temp);
        B[r] *=(F(1)/temp);
        det *=temp;
     for (size_t k=0; k<A.ColNo(); ++k)
        if(k!=r){
           temp=A(k,r);
           A[k]-=temp*A[r];
           B[k]-=temp*B[r];
        }
   }
   return pair <F, Matrix<F> >(det, B);
}

template <class F> inline pair<F, Matrix<F> >
Matrix<F>::inverse(arith_approx) const {
   Matrix<F> A(md->Row,md->Col);
        for(size_t i=0; i<md->Row; i++)
           for(size_t j=0; j<md->Col; j++)
                A.md->element[i][j]=md->element[i][j];

   Matrix<F> B=Matrix<F>::identity (ColNo());
   F det(1);
   for (size_t r=0; r<ColNo(); ++r){
        if ( abs(A(r,r))<=Field_traits<F>::insignificant()){
           for(size_t q=r+1; q<ColNo(); ++q)
                if(abs(A(q,r))>Field_traits<F>::insignificant()){
                   A[r].swap(A[q]);
                   B[r].swap(B[q]);
                   det=-det;
                }
           if(abs(A(r,r))<=Field_traits<F>::insignificant())
                return pair<F, Matrix<F> >(0, *this);
        }

       F temp=A(r,r);
       A[r] *=(F(1)/temp);
       B[r] *=(F(1)/temp);
       det *=temp;
  for (size_t k=0; k<ColNo(); ++k)
       if(k!=r){
          temp=A(k,r);
          A[k]-=temp*A[r];
          B[k]-=temp*B[r];
       }
   }

   return pair <F, Matrix<F> >(det, B);
}

//RowBasis
template <class F> inline Matrix<F>
Matrix<F>::rowBasis(arith_exact) const {

   Matrix<F> A(md->Row, md->Col);
        for(size_t i=0; i<md->Row; i++)
           for(size_t j=0; j<md->Col; j++)
                A.md->element[i][j]=md->element[i][j];

   size_t r=0, c=0;

   while (r<A.RowNo() && c<A.ColNo()){
        size_t k;
        for(k=r; k<RowNo() && A(k,c)==F(0); k++);
        if(k<A.RowNo()){
           if(k>r)
                A[r].swap(A[k]);
           A[r]*=(F(1)/A(r,c));

           for (k=0; k<RowNo(); ++k)
                if (k!=r)
                   A[k]-=A(k,c)*A[r];
           ++r;
        }
        ++c;
  }
  Matrix<F> K=A(0, r-1, 0, A.ColNo()-1);
  return K;
}

template <class F> inline Matrix<F>
Matrix<F>::rowBasis(arith_approx) const {

   Matrix<F> A(md->Row, md->Col);
        for(size_t i=0; i<md->Row; i++)
           for(size_t j=0; j<md->Col; j++)
                A.md->element[i][j]=md->element[i][j];

   size_t r=0, c=0;
   while (r<A.RowNo() && c<A.ColNo()){
        size_t k;
        for(k=r; k<RowNo()&& abs(A(k,c))<=Field_traits<F>
           ::insignificant(); k++);
        if(k<A.RowNo()){
           if(k>r)
                A[r].swap(A[k]);
           A[r]*=(F(1)/A(r,c));
           for (k=0; k<RowNo(); ++k)
                if (k!=r)
                   A[k]-=A(k,c)*A[r];
           ++r;
        }
        ++c;
  }
  Matrix<F> K=A(0, r-1, 0, A.ColNo()-1);
  return K;
}