Matrix.inl

Go to the documentation of this file.

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



namespace CIMPL
{


// Template Implementation


template< class T >
Matrix<T>::Matrix()
        : ndims(2), length(0)
{
        data = 0;
        xDim = 0;
        yDim = 0;
        columns = 0;
        clean = 0;

}


template< class T >
Matrix<T>::Matrix(const int rows, const int cols)
{
        if(rows < 1 || cols < 1)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
        }

        ndims = 2;
        length = rows*cols;
        xDim = cols;
        yDim = rows;
        data = new (std::nothrow) T[length];
        Utility::CheckPointer(data);

        columns = new (std::nothrow) Vector<T>[cols];
        Utility::CheckPointer(columns);
        for(int i=0; i<cols; i++)
        {
                columns[i].Set(&(data[i*rows]), rows);
        }
        
        clean = new (std::nothrow) Cleaner<T>(data, columns);
        Utility::CheckPointer(clean);
        
}

template< class T >
Matrix<T>::Matrix(const int rows, const int cols, T init)
{
        if(rows < 1 || cols < 1)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
        }

        ndims = 2;
        length = rows*cols;
        xDim = cols;
        yDim = rows;
        data = new (std::nothrow) T[length];
        Utility::CheckPointer(data);

        columns = new (std::nothrow) Vector<T>[cols];
        Utility::CheckPointer(columns);
        for(int i=0; i<cols; i++)
        {
                columns[i].Set(&(data[i*rows]), rows);
        }
        
        clean = new (std::nothrow) Cleaner<T>(data, columns);
        Utility::CheckPointer(clean);
        
        for(int i=0;i<length;i++)
        {
                data[i] = init;
        }

}




template< class T >
Matrix<T>::Matrix(string str)
{
        if(str.substr(0,1) == "[" && str.substr(str.size()-1,1) == "]")
        {
                str = str.substr(1,str.size()-2);
                vector<string> row_strs = Utility::Split(str, ";");
                if(row_strs.size() < 1 )
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
                }
                int rsize = (int)row_strs.size();
                
                vector<string> dummy = Utility::Split(row_strs[0]);
                if(dummy.size() < 1 )
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
                }
                int csize = (int)dummy.size();
                ndims = 2;
                length = rsize*csize;
                xDim = csize;
                yDim = rsize;
                data = new (std::nothrow) T[length];
                Utility::CheckPointer(data);
                columns = new (std::nothrow) Vector<T>[csize];
                Utility::CheckPointer(columns);
                for(int i=0; i<csize; i++)
                {
                        columns[i].Set(&(data[i*rsize]), rsize);
                }
                clean = new (std::nothrow) Cleaner<T>(data, columns);
                Utility::CheckPointer(clean);
                
                for(int i=0; i<rsize; i++)
                {
                        vector<string> elem_strs = Utility::Split(row_strs[i]);
                        if(elem_strs.size() != xDim )
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("String parse error. Matrix dimensions are not consistent!");
                        }
                        for(int j=0; j<csize; j++)
                        {
                                columns[j].data[i] = (T)Utility::ToDouble(elem_strs[j]);
                        }
                }
        }
        else
        {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Incorrect initialization. String cannot be parsed!");
        }
}





template< class T >
Matrix<T>::Matrix(Matrix<T> &m) 
{
        ndims = m.ndims;
        length = m.length;
        xDim = m.xDim;
        yDim = m.yDim;

        data = m.data;
        columns = m.columns;

        clean = new (std::nothrow) Cleaner<T>(data, columns);
        Utility::CheckPointer(clean);

}



template< class T >
Matrix<T>::~Matrix()
{
        
        if(clean != 0)
        {
                delete clean;
        }
}


template< class T >
void Matrix<T>::Clean()
{
        if(clean != 0)
        {
                delete clean;
                data = 0;
                columns = 0;
                clean = 0;
        }
}


template< class T >
void Matrix<T>::Set(const int rows, const int cols, T *_data)
{
        if(rows < 1 || cols < 1)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("All Matrix dimensions should be larger than 0!");
        }

        ndims = 2;
        length = rows*cols;
        xDim = cols;
        yDim = rows;
        data = _data;

        columns = new (std::nothrow) Vector<T>[cols];
        Utility::CheckPointer(columns);
        for(int i=0; i<cols; i++)
        {
                columns[i].Set(&(data[i*rows]), rows);
        }
        
        delete clean;
        clean = new (std::nothrow) Cleaner<T>(data, columns);
        Utility::CheckPointer(clean);
        
}

template< class T >
inline const T* Matrix<T>::DataPtr()
{
        return data;
}

template< class T >
inline T* Matrix<T>::Data()
{
        return data;
}


template< class T >
Matrix<T> Matrix<T>::Clone() const
{
        Matrix<T> temp(Rows(), Columns());
        
        memcpy(temp.data, data, sizeof(T)*temp.length);

        return temp;
}


template< class T >
void Matrix<T>::SwitchColumns(int i, int j)
{
        T *temp_i = columns[i].data;
        columns[i].data = columns[j].data;
        columns[j].data = temp_i;
}


template< class T >
void Matrix<T>::SyncData2Columns()
{
        T* temp_data = new (std::nothrow) T[length];
        Utility::CheckPointer(temp_data);
        
        for(int i=0; i<xDim; i++)
        {
                memcpy(&(temp_data[i*yDim]), columns[i].data, sizeof(T)*yDim);
        }
        data = temp_data;
        
        columns = new (std::nothrow) Vector<T>[xDim];
        Utility::CheckPointer(columns);
        
        for(int i=0; i<xDim; i++)
        {
                columns[i].data = &(data[i*yDim]);
        }
        
        delete clean;
        clean = new (std::nothrow) Cleaner<T>(data, columns);
        Utility::CheckPointer(clean);
}


template< class T >
Matrix<T> Matrix<T>::Slice(const int row1, const int row2, const int col1, const int col2)
{
        if(row1<0 || row1>=Rows() || row2<0 || row2>=Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Index outside bounds!");
        }
        if(row1 > row2)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
        }


        if(col1<0 || col1>=Columns() || col2<0 || col2>=Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Index outside bounds!");
        }
        if(col1 > col2)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
        }

        
        Matrix<T> temp(row2-row1+1, col2-col1+1);
        
        for(int j=col1; j<=col2; j++)
        {
                memcpy(temp[ j-col1 ].data, &(columns[j].data[row1]), sizeof(T)*temp.Rows());
        }

        return temp;
}


//template< class T >
//SubMatrix<T> Matrix<T>::Slice(int row1, int row2, int col1, int col2)
//{
//      if(row1<0 || row1>=YDim() || row2<0 || row2>=YDim())
//      {
//              cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//              Utility::RunTimeError("Index outside bounds!");
//      }
//      if(row1 > row2)
//      {
//              cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//              Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
//      }
//
//
//      if(col1<0 || col1>=XDim() || col2<0 || col2>=XDim())
//      {
//              cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//              Utility::RunTimeError("Index outside bounds!");
//      }
//      if(col1 > col2)
//      {
//              cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
//              Utility::RunTimeError("Second slice parameter cannot be less than the first parameter!");
//      }
//
//      
//      SubMatrix<T> temp;
//      temp.xDim = col2-col1+1;
//      temp.yDim = row2-row1+1;
//      temp.columns = new (std::nothrow) Vector<T>[temp.xDim];
//      Utility::CheckPointer(temp.columns);
//      temp.clean = new (std::nothrow) Cleaner<T>(temp.columns);
//      Utility::CheckPointer(temp.clean);
//
//      for(int j=col1; j<=col2; j++)
//      {
//              temp.columns[j-col1].Set(&(columns[j].data[row1]),  row2-row1+1);
//      }
//
//      return temp;
//}



template< class T >
Matrix<T> Matrix<T>::Slice(string str1, string str2)
{
        int row1, row2, col1, col2;
        vector<string> bounds1 = Utility::Split(str1, ":");
        vector<string> bounds2 = Utility::Split(str2, ":");
        
        if(str1 == ":")
        {
                row1 = 0;
                row2 = yDim-1;
        }
        else if(bounds1.size() == 2)
        {
                row1 = Utility::ToInt(bounds1[0]);
                row2 = Utility::ToInt(bounds1[1]);
        }
        else
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrect slice argument. String cannot be parsed!");
        }

        if(str2 == ":")
        {
                col1 = 0;
                col2 = xDim-1;
        }
        else if(bounds2.size() == 2)
        {
                col1 = Utility::ToInt(bounds2[0]);
                col2 = Utility::ToInt(bounds2[1]);
        }
        else
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrect slice argument. String cannot be parsed!");
        }
        
        return this->Slice(row1, row2, col1, col2);

}



template< class T >
Vector<T> Matrix<T>::Slice(string str)
{
        Vector<T> temp = (Vector<T>)*this;
        return temp.Slice(str);
}





template< class T >
inline const int Matrix<T>::Columns() const
{
        return xDim;
}

template< class T >
inline const int Matrix<T>::Rows() const
{
        return yDim;
}


template< class T >
inline const int Matrix<T>::XDim() const
{
        return xDim;
}

template< class T >
inline const int Matrix<T>::YDim() const
{
        return yDim;
}




template< class T >
inline const int Matrix<T>::Length() const
{
        return length;
}

template< class T >
inline const int Matrix<T>::Numel() const
{
        return length;
}

template< class T >
inline const int Matrix<T>::NDims() const
{
        return ndims;
}

template< class T >
inline void Matrix<T>::Init(const T init)
{
        for(int i=0;i<length;i++)
        {
                data[i] = init;
        }
}

//template< class T >
//inline Matrix<T>&  Matrix<T>::Rand()
//{
//      if(!RandomGen::Initialized())
//      {
//              RandomGen::Initialize();
//      }
//      for(int i=0;i<length;i++)
//      {
//              data[i] = (T)rand();
//      }
//      return *this;
//}

template< class T >
inline Matrix<T>&  Matrix<T>::Rand(const double max)
{
        if(!RandomGen::Initialized())
        {
                RandomGen::Initialize();
        }
        for(int i=0;i<length;i++)
        {
                data[i] = (T)(rand()/(double)RAND_MAX*max);
        }
        return *this;
}


//template< class T >
//Matrix<T> Matrix<T>::Rand(const int rows, const int cols)
//{
//      Matrix<T> m(rows, cols);
//      if(!RandomGen::Initialized())
//      {
//              RandomGen::Initialize();
//      }
//      for(int i=0;i<rows*cols;i++)
//      {
//              m.data[i] = (T)rand();
//      }
//      return m;
//}

template< class T >
Matrix<T> Matrix<T>::Rand(const int rows, const int cols, const double max)
{
        Matrix<T> m(rows, cols);
        if(!RandomGen::Initialized())
        {
                RandomGen::Initialize();
        }
        for(int i=0;i<rows*cols;i++)
        {
                m.data[i] = (T)(rand()/(double)RAND_MAX*max);
        }
        return m;
}



template< class T >
void Matrix<T>::ReadFromMatrix(const Matrix<T>& m)
{
        if(xDim < m.xDim || yDim < m.yDim)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("You can't read from a larger Matrix to a smaller Matrix!");
        }

        for(int i=0;i<m.xDim;i++)
        {
                for(int j=0;j<m.yDim;j++)
                {
                        columns[i].data[j] = m.columns[i].data[j];
                }
        }

}

template< class T >
void Matrix<T>::ReadFromMatrix(const Matrix<T>& m, const int rowStart, const int colStart)
{
        if(xDim < m.xDim+colStart || yDim < m.yDim+rowStart)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("You can't read from a larger Matrix (from the copy start point) to a smaller Matrix!");
        }

        for(int i=0;i<m.xDim;i++)
        {
                for(int j=0;j<m.yDim;j++)
                {
                        columns[i+colStart].data[j+rowStart] = m.columns[i].data[j];
                }
        }

}









template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Matrix<T>& m1, Matrix<T>& m2)
{
        Matrix<T> temp;
        if(dimension == 1)
        {
                if(m1.Columns() != m2.Columns())
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Matrix sizes are not compatible for concatenation!");
                }
                
                temp = Matrix<T>(m1.Rows()+m2.Rows(), m1.Columns());
                temp.ReadFromMatrix(m1);
                temp.ReadFromMatrix(m2, m1.Rows(), 0);
                
        }
        else if(dimension == 2)
        {
                if(m1.Rows() != m2.Rows())
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Matrix sizes are not compatible for concatenation!");
                }

                temp = Matrix<T>(m1.Rows(), m1.Columns()+m2.Columns());
                temp.ReadFromMatrix(m1);
                temp.ReadFromMatrix(m2, 0, m1.Columns());
        }
        else
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
        }
        
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Matrix<T>& m1, Vector<T>& m2)
{
        return Matrix<T>::Cat(dimension, m1, (Matrix<T>)m2);
}

template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Vector<T>& m1, Matrix<T>& m2)
{
        return Matrix<T>::Cat(dimension, (Matrix<T>)m1, m2);
}

template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Vector<T>& m1, Vector<T>& m2)
{
        return Matrix<T>::Cat(dimension, (Matrix<T>)m1, (Matrix<T>)m2);
}


template< class T >
Matrix<T> Matrix<T>::Ones(int side)
{
        Matrix<T> temp(side,side,1);
        return temp;
}


template< class T >
Matrix<T> Matrix<T>::Ones(int rows, int cols)
{
        Matrix<T> temp(rows,cols,1);
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Zeros(int side)
{
        Matrix<T> temp(side,side,0);
        return temp;
}


template< class T >
Matrix<T> Matrix<T>::Zeros(int rows, int cols)
{
        Matrix<T> temp(rows,cols,0);
        return temp;
}







template< class T >
bool Matrix<T>::IsSquare(Matrix<T>& m)
{
        if(m.xDim == m.yDim)
        {
                return true;
        }
        else
        {
                return false;
        }
}





template< class T >
bool Matrix<T>::IsM2MCompatible(Matrix<T>& m1, Matrix<T>& m2)
{
        if(m1.xDim == m2.yDim)
        {
                return true;
        }
        else
        {
                return false;
        }
}


template< class T >
Matrix<T> Matrix<T>::MMultiply(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsM2MCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!");
        }

        Matrix<T> temp(m1.yDim, m2.xDim);
        temp.Init(0);
        for(int i=0; i<m1.yDim; i++)
        {
                for(int j=0; j<m2.xDim; j++)
                {
                        for(int k = 0; k<m1.xDim; k++)
                        {
                                temp[j][i] += m1[k][i]*m2[j][k]; 
                        }
                }
        }
        
        return temp;
}


template< class T >
Matrix<T> Matrix<T>::MMultiply(Matrix<T>& m1, Vector<T>& m2)
{
        if(m1.xDim != m2.length)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!");
        }

        Matrix<T> temp(m1.yDim, 1);
        temp.Init(0);
        for(int i=0; i<m1.yDim; i++)
        {
                for(int j=0; j<m2.length; j++)
                {
                        temp.data[i] += m1[j][i]*m2[j]; 
                }
        }
        
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::MMultiply(Vector<T>& m1, Matrix<T>& m2)
{
        if(m2.Rows() != m1.length)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!");
        }

        Matrix<T> temp(1, m2.Columns());
        temp.Init(0);
        for(int i=0; i<m2.Columns(); i++)
        {
                for(int j=0; j<m1.length; j++)
                {
                                temp.data[i] += m2[i][j]*m1[j]; 
                }
        }
        
        return temp;
}


template< class T >
Matrix<T> Matrix<T>::Transpose(Matrix<T>& m)
{
        Matrix<T> temp(m.Columns(), m.Rows());
        for(int i=0; i<m.Rows(); i++)
        {
                for(int j=0; j<m.Columns(); j++)
                {
                        temp[i][j] = m[j][i]; 
                }
        }
        return temp;
}


template< class T >
Matrix<T>& Matrix<T>::Transpose()
{
        if(!Matrix<T>::IsSquare(*this))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        for(int i=0; i<yDim; i++)
        {
                for(int j=0; j<i; j++)
                {
                        T temp = columns[j][i];
                        columns[j][i] = columns[i][j];
                        columns[i][j] = temp;
                }
        }
        return *this;

}


template< class T >
bool Matrix<T>::IsCompatible(Matrix<T>& m1, Matrix<T>& m2)
{
        if(m1.xDim != m2.xDim || m1.yDim != m2.yDim)
        {
                return false;
        }
        else
        {
                return true;
        }

}


// //////////////////
// Boolean Operations...
// //////////////////

template< class T >
Matrix<int> Matrix<T>::And(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] != 0 && m2.data[i] != 0) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Or(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] == 0 && m2.data[i] == 0) ? 0 : 1;
        }
        
        return temp;
}

template< class T >
Matrix<int> Matrix<T>::Lt(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] < m2.data[i]) ? 1 : 0;
        }
        
        return temp;
}

template< class T >
Matrix<int> Matrix<T>::Gt(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] > m2.data[i]) ? 1 : 0;
        }
        
        return temp;
}

template< class T >
Matrix<int> Matrix<T>::Le(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] <= m2.data[i]) ? 1 : 0;
        }
        
        return temp;
}

template< class T >
Matrix<int> Matrix<T>::Ge(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] >= m2.data[i]) ? 1 : 0;
        }
        
        return temp;
}

template< class T >
Matrix<int> Matrix<T>::Eq(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] == m2.data[i]) ? 1 : 0;
        }
        
        return temp;
}

template< class T >
Matrix<int> Matrix<T>::Ne(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] != m2.data[i]) ? 1 : 0;
        }
        
        return temp;
}




// ////////////////
// Boolean Operations with value types...
// ////////////////


template< class T >
Matrix<int> Matrix<T>::And(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] != 0 && v != 0) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Or(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] == 0 && v == 0) ? 0 : 1;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Lt(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] < v) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Gt(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] > v) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Le(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] <= v) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Ge(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] >= v) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Eq(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] == v) ? 1 : 0;
        }
        
        return temp;
}


template< class T >
Matrix<int> Matrix<T>::Ne(Matrix<T>& m1, T v)
{
        Matrix<int> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.Length();i++)
        {
                temp.Data()[i] = (m1.data[i] != v) ? 1 : 0;
        }
        
        return temp;
}











// /////////////////
// Elementwise matrix arithmetic
// ////////////////


template< class T >
Matrix<T> Matrix<T>::Add(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] + m2.data[i];
        }
        
        return temp;
}


template< class T >
Matrix<T> Matrix<T>::Subtract(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] - m2.data[i];
        }
        
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Multiply(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] * m2.data[i];
        }
        
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Divide(Matrix<T>& m1, Matrix<T>& m2)
{
        if(!Matrix<T>::IsCompatible(m1, m2))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                //FW: catch division exception instead...
                if(m2.data[i] != 0)
                {
                        temp.data[i] = m1.data[i] / m2.data[i];
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Divide by zero in matrix division!");
                }
        }
        
        return temp;
}

// ////////
// Arithmetic operations between a matrix and a value.
// ////////


template< class T >
Matrix<T> Matrix<T>::Add(Matrix<T>& m1, T v2)
{
        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] + v2;
        }
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Subtract(Matrix<T>& m1, T v2)
{
        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] - v2;
        }
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Subtract(T v2, Matrix<T>& m1)
{
        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = v2 - m1.data[i];
        }
        return temp;
}


template< class T >
Matrix<T> Matrix<T>::Multiply(Matrix<T>& m1, T v2)
{
        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] * v2;
        }
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Divide(Matrix<T>& m1, T v2)
{
        if(v2 == 0)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Divide by zero in matrix by value division!");
        }

        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                temp.data[i] = m1.data[i] / v2;
        }
        return temp;
}

template< class T >
Matrix<T> Matrix<T>::Divide(T v2, Matrix<T>& m1)
{
        Matrix<T> temp(m1.yDim, m1.xDim);
        for(int i=0;i<temp.length;i++)
        {
                if(m1.data[i] != 0)
                {
                        temp.data[i] = v2 / m1.data[i];
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Divide by zero in value by matrix division!");
                }
        }
        return temp;
}


// //////////
// Elementwise matrix arithmetic with "this" matrix
// /////////

template< class T >
Matrix<T>& Matrix<T>::Add(Matrix<T>& m)
{
        if(!Matrix<T>::IsCompatible(*this, m))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        for(int i=0;i<this->length;i++)
        {
                this->data[i] += m.data[i];
        }
        
        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::Subtract(Matrix<T>& m)
{
        if(!Matrix<T>::IsCompatible(*this, m))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        for(int i=0;i<this->length;i++)
        {
                this->data[i] -= m.data[i];
        }
        
        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::Multiply(Matrix<T>& m)
{
        if(!Matrix<T>::IsCompatible(*this, m))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        for(int i=0;i<this->length;i++)
        {
                this->data[i] *= m.data[i];
        }
        
        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::Divide(Matrix<T>& m)
{
        if(!Matrix<T>::IsCompatible(*this, m))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix sizes are not compatible!");
        }

        for(int i=0;i<this->length;i++)
        {
                if(m.data[i] != 0)
                {
                        this->data[i] /= m.data[i];
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Divide by zero in matrix division!");
                }
        }
        
        return *this;
}

// /////////////
// Arithmetic operation between "this" matrix and a value
// /////////////

template< class T >
Matrix<T>& Matrix<T>::Add(T v)
{
        for(int i=0;i<this->length;i++)
        {
                this->data[i] += v;
        }
        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::Subtract(T v)
{
        for(int i=0;i<this->length;i++)
        {
                this->data[i] -= v;
        }
        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::Multiply(T v)
{
        for(int i=0;i<this->length;i++)
        {
                this->data[i] *= v;
        }
        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::Divide(T v)
{
        if(v == 0)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Divide by zero in matrix by value division!");
        }
        
        for(int i=0;i<this->length;i++)
        {
                this->data[i] /= v;
        }
        return *this;
}






// /////////
// OPERATORS
// /////////


template< class T >
Matrix<T>& Matrix<T>::operator= (Matrix<T>& m)
{
        ndims = m.ndims;
        length = m.length;
        xDim = m.xDim;
        yDim = m.yDim;
        data = m.data;
        columns = m.columns;
        
        delete clean;
        clean = new (std::nothrow) Cleaner<T>(data, columns);
        Utility::CheckPointer(clean);
        
        
        return *this;
}


template< class T >
Matrix<T>& Matrix<T>::operator= (Array<T>& m)
{
        ndims = 2;
        xDim = 0;
        yDim = 0;
        data = 0;
        columns = 0;

        if(m.ndims == 2){
                length = m.length;
                data = m.data;
                xDim = m.xDim;
                yDim = m.yDim;

                columns = new (std::nothrow) Vector<T>[Columns()];
                Utility::CheckPointer(columns);
                for(int i=0; i<Columns(); i++)
                {
                        columns[i].Set(&(data[i*Rows()]), Rows());
                }

                delete clean;
                clean = new (std::nothrow) Cleaner<T>(data, columns);
                Utility::CheckPointer(clean);
        }
        else if(m.ndims > 0) // Convert to row matrux
        {
                length = m.length;
                data = m.data;
                xDim = length;
                yDim = 1;
                
                columns = new (std::nothrow) Vector<T>[1];
                Utility::CheckPointer(columns);
                columns[0].Set(&(data[0]),length);

                delete clean;
                clean = new (std::nothrow) Cleaner<T>(data, columns);
                Utility::CheckPointer(clean);
                
                Utility::Warning("Array (not of dimension 2) is converted to row Matrix. dimensionality information is lost.");
        }
        

        return *this;
}

template< class T >
Matrix<T>& Matrix<T>::operator= (Vector<T>& m)
{
        *this = (Matrix<T>)m;
        return *this;
}


template< class T >
Matrix<T>& Matrix<T>::operator= (SubMatrix<T>& m)
{
        if(xDim != m.xDim || yDim != m.yDim)
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("All Matrix and SubMatrix dimensions should agree!");
        }
        for(int i=0;i<m.xDim;i++)
        {
                for(int j=0;j<m.yDim;j++)
                {
                        columns[i].data[j] = m[i][j];
                }
        }
        
        return *this;
}



template< class T >
Matrix<T>& Matrix<T>::operator= (string str)
{
        Matrix<T> temp(str);
        *this = temp;
        return *this;
}


// ////
//Unary operators
// ////

template< class T >
inline Matrix<T> Matrix<T>::operator+ ()
{
        return *this;
}


template< class T >
inline