www.pudn.com > Image_segment.rar > Matrix.inl

//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. 
 
 
 
// Template Implementation 
 
 
template< class T > 
Matrix::Matrix() 
	: ndims(2), length(0) 
{ 
	data = 0; 
	xDim = 0; 
	yDim = 0; 
	columns = 0; 
	clean = 0; 
 
} 
 
 
template< class T > 
Matrix::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[cols]; 
	Utility::CheckPointer(columns); 
	for(int i=0; i(data, columns); 
	Utility::CheckPointer(clean); 
	 
} 
 
template< class T > 
Matrix::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[cols]; 
	Utility::CheckPointer(columns); 
	for(int i=0; i(data, columns); 
	Utility::CheckPointer(clean); 
	 
	for(int i=0;i 
Matrix::Matrix(string str) 
{ 
	if(str.substr(0,1) == "[" && str.substr(str.size()-1,1) == "]") 
	{ 
		str = str.substr(1,str.size()-2); 
		vector 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 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[csize]; 
		Utility::CheckPointer(columns); 
		for(int i=0; i(data, columns); 
		Utility::CheckPointer(clean); 
		 
		for(int i=0; i 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 
Matrix::Matrix(Matrix &m)  
{ 
	ndims = m.ndims; 
	length = m.length; 
	xDim = m.xDim; 
	yDim = m.yDim; 
 
	data = m.data; 
	columns = m.columns; 
 
	clean = new (std::nothrow) Cleaner(data, columns); 
	Utility::CheckPointer(clean); 
 
} 
 
 
 
template< class T > 
Matrix::~Matrix() 
{ 
	 
	if(clean != 0) 
	{ 
		delete clean; 
	} 
} 
 
 
template< class T > 
void Matrix::Clean() 
{ 
	if(clean != 0) 
	{ 
		delete clean; 
		data = 0; 
		columns = 0; 
		clean = 0; 
	} 
} 
 
 
template< class T > 
void Matrix::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[cols]; 
	Utility::CheckPointer(columns); 
	for(int i=0; i(data, columns); 
	Utility::CheckPointer(clean); 
	 
} 
 
template< class T > 
inline const T* Matrix::DataPtr() 
{ 
	return data; 
} 
 
template< class T > 
inline T* Matrix::Data() 
{ 
	return data; 
} 
 
 
template< class T > 
Matrix Matrix::Clone() const 
{ 
	Matrix temp(Rows(), Columns()); 
	 
	memcpy(temp.data, data, sizeof(T)*temp.length); 
 
	return temp; 
} 
 
 
template< class T > 
void Matrix::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::SyncData2Columns() 
{ 
	T* temp_data = new (std::nothrow) T[length]; 
	Utility::CheckPointer(temp_data); 
	 
	for(int i=0; i[xDim]; 
	Utility::CheckPointer(columns); 
	 
	for(int i=0; i(data, columns); 
	Utility::CheckPointer(clean); 
} 
 
 
template< class T > 
Matrix Matrix::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 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 Matrix::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 temp; 
//	temp.xDim = col2-col1+1; 
//	temp.yDim = row2-row1+1; 
//	temp.columns = new (std::nothrow) Vector[temp.xDim]; 
//	Utility::CheckPointer(temp.columns); 
//	temp.clean = new (std::nothrow) Cleaner(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 Matrix::Slice(string str1, string str2) 
{ 
	int row1, row2, col1, col2; 
	vector bounds1 = Utility::Split(str1, ":"); 
	vector 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 Matrix::Slice(string str) 
{ 
	Vector temp = (Vector)*this; 
	return temp.Slice(str); 
} 
 
 
 
 
 
template< class T > 
inline const int Matrix::Columns() const 
{ 
	return xDim; 
} 
 
template< class T > 
inline const int Matrix::Rows() const 
{ 
	return yDim; 
} 
 
 
template< class T > 
inline const int Matrix::XDim() const 
{ 
	return xDim; 
} 
 
template< class T > 
inline const int Matrix::YDim() const 
{ 
	return yDim; 
} 
 
 
 
 
template< class T > 
inline const int Matrix::Length() const 
{ 
	return length; 
} 
 
template< class T > 
inline const int Matrix::Numel() const 
{ 
	return length; 
} 
 
template< class T > 
inline const int Matrix::NDims() const 
{ 
	return ndims; 
} 
 
template< class T > 
inline void Matrix::Init(const T init) 
{ 
	for(int i=0;i 
//inline Matrix&  Matrix::Rand() 
//{ 
//	if(!RandomGen::Initialized()) 
//	{ 
//		RandomGen::Initialize(); 
//	} 
//	for(int i=0;i 
inline Matrix&  Matrix::Rand(const double max) 
{ 
	if(!RandomGen::Initialized()) 
	{ 
		RandomGen::Initialize(); 
	} 
	for(int i=0;i 
//Matrix Matrix::Rand(const int rows, const int cols) 
//{ 
//	Matrix m(rows, cols); 
//	if(!RandomGen::Initialized()) 
//	{ 
//		RandomGen::Initialize(); 
//	} 
//	for(int i=0;i 
Matrix Matrix::Rand(const int rows, const int cols, const double max) 
{ 
	Matrix m(rows, cols); 
	if(!RandomGen::Initialized()) 
	{ 
		RandomGen::Initialize(); 
	} 
	for(int i=0;i 
void Matrix::ReadFromMatrix(const Matrix& 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 
void Matrix::ReadFromMatrix(const Matrix& 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 
Matrix Matrix::Cat(int dimension, Matrix& m1, Matrix& m2) 
{ 
	Matrix 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(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(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 Matrix::Cat(int dimension, Matrix& m1, Vector& m2) 
{ 
	return Matrix::Cat(dimension, m1, (Matrix)m2); 
} 
 
template< class T > 
Matrix Matrix::Cat(int dimension, Vector& m1, Matrix& m2) 
{ 
	return Matrix::Cat(dimension, (Matrix)m1, m2); 
} 
 
template< class T > 
Matrix Matrix::Cat(int dimension, Vector& m1, Vector& m2) 
{ 
	return Matrix::Cat(dimension, (Matrix)m1, (Matrix)m2); 
} 
 
 
template< class T > 
Matrix Matrix::Ones(int side) 
{ 
	Matrix temp(side,side,1); 
	return temp; 
} 
 
 
template< class T > 
Matrix Matrix::Ones(int rows, int cols) 
{ 
	Matrix temp(rows,cols,1); 
	return temp; 
} 
 
template< class T > 
Matrix Matrix::Zeros(int side) 
{ 
	Matrix temp(side,side,0); 
	return temp; 
} 
 
 
template< class T > 
Matrix Matrix::Zeros(int rows, int cols) 
{ 
	Matrix temp(rows,cols,0); 
	return temp; 
} 
 
 
 
 
 
 
 
template< class T > 
bool Matrix::IsSquare(Matrix& m) 
{ 
	if(m.xDim == m.yDim) 
	{ 
		return true; 
	} 
	else 
	{ 
		return false; 
	} 
} 
 
 
 
 
 
/// \brief are these matrices (and in this order) are compatible for matrix multiplication 
template< class T > 
bool Matrix::IsM2MCompatible(Matrix& m1, Matrix& m2) 
{ 
	if(m1.xDim == m2.yDim) 
	{ 
		return true; 
	} 
	else 
	{ 
		return false; 
	} 
} 
 
 
/// \brief Matrix multiplication. Also used to overload & operator. 
template< class T > 
Matrix Matrix::MMultiply(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsM2MCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!"); 
	} 
 
	Matrix temp(m1.yDim, m2.xDim); 
	temp.Init(0); 
	for(int i=0; i 
Matrix Matrix::MMultiply(Matrix& m1, Vector& m2) 
{ 
	if(m1.xDim != m2.length) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!"); 
	} 
 
	Matrix temp(m1.yDim, 1); 
	temp.Init(0); 
	for(int i=0; i 
Matrix Matrix::MMultiply(Vector& m1, Matrix& m2) 
{ 
	if(m2.Rows() != m1.length) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!"); 
	} 
 
	Matrix temp(1, m2.Columns()); 
	temp.Init(0); 
	for(int i=0; i 
Matrix Matrix::Transpose(Matrix& m) 
{ 
	Matrix temp(m.Columns(), m.Rows()); 
	for(int i=0; i 
Matrix& Matrix::Transpose() 
{ 
	if(!Matrix::IsSquare(*this)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix is not square!"); 
	} 
	 
	for(int i=0; i 
bool Matrix::IsCompatible(Matrix& m1, Matrix& m2) 
{ 
	if(m1.xDim != m2.xDim || m1.yDim != m2.yDim) 
	{ 
		return false; 
	} 
	else 
	{ 
		return true; 
	} 
 
} 
 
 
// ////////////////// 
// Boolean Operations... 
// ////////////////// 
 
/// \brief Elementwise AND operator. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::And(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Or(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Lt(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i (greater than) operator. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Gt(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i m2.data[i]) ? 1 : 0; 
	} 
	 
	return temp; 
} 
 
/// \brief Elementwise <= (less than or equal) operator. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Le(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i= (greater than or equal) operator. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Ge(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i= m2.data[i]) ? 1 : 0; 
	} 
	 
	return temp; 
} 
 
/// \brief Elementwise == operator. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Eq(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Ne(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::And(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Or(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Lt(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i (Greater than) operator between Matrix elements and a value. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Gt(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i v) ? 1 : 0; 
	} 
	 
	return temp; 
} 
 
 
/// \brief <= (Less than or equal) operator between Matrix elements and a value. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Le(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Ge(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i= v) ? 1 : 0; 
	} 
	 
	return temp; 
} 
 
 
/// \brief == operator between Matrix elements and a value. Returns 1 or 0 for each element. 
template< class T > 
Matrix Matrix::Eq(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Ne(Matrix& m1, T v) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Add(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Subtract(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Multiply(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Divide(Matrix& m1, Matrix& m2) 
{ 
	if(!Matrix::IsCompatible(m1, m2)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Add(Matrix& m1, T v2) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Subtract(Matrix& m1, T v2) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Subtract(T v2, Matrix& m1) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Multiply(Matrix& m1, T v2) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Divide(Matrix& m1, T v2) 
{ 
	if(v2 == 0) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Divide by zero in matrix by value division!"); 
	} 
 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix Matrix::Divide(T v2, Matrix& m1) 
{ 
	Matrix temp(m1.yDim, m1.xDim); 
	for(int i=0;i 
Matrix& Matrix::Add(Matrix& m) 
{ 
	if(!Matrix::IsCompatible(*this, m)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	for(int i=0;ilength;i++) 
	{ 
		this->data[i] += m.data[i]; 
	} 
	 
	return *this; 
} 
 
template< class T > 
Matrix& Matrix::Subtract(Matrix& m) 
{ 
	if(!Matrix::IsCompatible(*this, m)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	for(int i=0;ilength;i++) 
	{ 
		this->data[i] -= m.data[i]; 
	} 
	 
	return *this; 
} 
 
template< class T > 
Matrix& Matrix::Multiply(Matrix& m) 
{ 
	if(!Matrix::IsCompatible(*this, m)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	for(int i=0;ilength;i++) 
	{ 
		this->data[i] *= m.data[i]; 
	} 
	 
	return *this; 
} 
 
template< class T > 
Matrix& Matrix::Divide(Matrix& m) 
{ 
	if(!Matrix::IsCompatible(*this, m)) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Matrix sizes are not compatible!"); 
	} 
 
	for(int i=0;ilength;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& Matrix::Add(T v) 
{ 
	for(int i=0;ilength;i++) 
	{ 
		this->data[i] += v; 
	} 
	return *this; 
} 
 
template< class T > 
Matrix& Matrix::Subtract(T v) 
{ 
	for(int i=0;ilength;i++) 
	{ 
		this->data[i] -= v; 
	} 
	return *this; 
} 
 
template< class T > 
Matrix& Matrix::Multiply(T v) 
{ 
	for(int i=0;ilength;i++) 
	{ 
		this->data[i] *= v; 
	} 
	return *this; 
} 
 
template< class T > 
Matrix& Matrix::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;ilength;i++) 
	{ 
		this->data[i] /= v; 
	} 
	return *this; 
} 
 
 
 
 
 
 
// ///////// 
// OPERATORS 
// ///////// 
 
 
template< class T > 
Matrix& Matrix::operator= (Matrix& m) 
{ 
	if(&m != this) 
	{ 
		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(data, columns); 
		Utility::CheckPointer(clean); 
	}	 
	 
	return *this; 
} 
 
 
//template< class T > 
//Matrix& Matrix::operator= (Array& 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[Columns()]; 
//		Utility::CheckPointer(columns); 
//		for(int i=0; i(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[1]; 
//		Utility::CheckPointer(columns); 
//		columns[0].Set(&(data[0]),length); 
// 
//		delete clean; 
//		clean = new (std::nothrow) Cleaner(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& Matrix::operator= (Vector& m) 
{ 
	*this = (Matrix)m; 
	return *this; 
} 
 
 
template< class T > 
Matrix& Matrix::operator= (SubMatrix& 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 
Matrix& Matrix::operator= (string str) 
{ 
	Matrix temp(str); 
	*this = temp; 
	return *this; 
} 
 
 
// //// 
//Unary operators 
// //// 
 
/// \brief Unary +. Does not have any effect. Returns the same matrix.  
template< class T > 
inline Matrix Matrix::operator+ () 
{ 
	return *this; 
} 
 
 
/// \brief Unary - 
template< class T > 
inline Matrix Matrix::operator- () 
{ 
	Matrix temp(yDim, xDim); 
	for(int i=0;i 
inline Matrix Matrix::operator! () 
{ 
	Matrix temp(yDim, xDim); 
	for(int i=0;i 
inline Matrix Matrix::operator~ () 
{ 
	return Matrix::Transpose(*this); 
} 
 
 
 
//--------------- 
 
/// \brief Matrix element access (sequential: column scan). Bounds are checked. 
template< class T > 
inline T& Matrix::Elem(const int i) 
{ 
	if(i<0 || i>=length) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Index outside bounds!"); 
	} 
	return data[i]; 
} 
 
/// \brief Matrix element access (sequential: column scan). Bounds are not checked. 
template< class T > 
inline T& Matrix::ElemNC(const int i) 
{ 
	return data[i]; 
} 
 
 
/// \brief Matrix element access (row,col). Bounds are checked. 
template< class T > 
inline T& Matrix::Elem(const int j, const int i) 
{ 
	if(i<0 || i>=xDim || j<0 || j>=yDim) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Index outside bounds!"); 
	} 
 
	return columns[i].data[j]; 
 
} 
 
 
/// \brief Matrix element access (row,col). Bounds are not checked. 
template< class T > 
inline T& Matrix::ElemNC(const int j, const int i) 
{ 
	return columns[i].data[j]; 
 
} 
 
 
/// \brief Matrix element access using coordinates (x,y). Bounds are checked. 
template< class T > 
inline T& Matrix::Coor(const int i, const int j) 
{ 
	if(i<0 || i>=dims[0] || j<0 || j>=dims[1]) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Index outside bounds!"); 
	} 
 
	return columns[i].data[j]; 
 
} 
 
 
/// \brief Matrix element access using coordinates (x,y). Bounds are not checked. 
template< class T > 
inline T& Matrix::CoorNC(const int i, const int j) 
{ 
	return columns[i].data[j]; 
 
} 
 
 
//--------------- 
 
 
/// \brief Matrix element access (sequential: column scan). Bounds are checked (See ElemNC() for element access with no bounds check). 
template< class T > 
inline T& Matrix::operator() (const int i) 
{ 
	if(i<0 || i>=length) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Index outside bounds!"); 
	} 
	return data[i]; 
} 
 
/// \brief Matrix element access (row,col). Bounds are checked (See ElemNC() for element access with no bounds check). 
template< class T > 
inline T& Matrix::operator() (const int j, const int i) 
{ 
#ifdef CIMPL_BOUNDS_CHECK 
	if(i<0 || i>=xDim || j<0 || j>=yDim) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("Index outside bounds!"); 
	} 
#endif 
 
	return columns[i].data[j]; 
 
} 
 
 
/// \brief Returns a pointer to the beginning of a column. 
template< class T > 
inline Vector& Matrix::operator[] (const int i) 
{ 
	return columns[i]; 
} 
 
 
 
 
 
template< class T > 
inline SubMatrix Matrix::operator() (const int row1, const int row2, const int col1, const 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 temp; 
	temp.xDim = col2-col1+1; 
	temp.yDim = row2-row1+1; 
	temp.columns = new (std::nothrow) Vector[temp.xDim]; 
	Utility::CheckPointer(temp.columns); 
	temp.clean = new (std::nothrow) Cleaner(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 > 
SubMatrix Matrix::operator() (string str1, string str2) 
{ 
	int row1, row2, col1, col2; 
	vector bounds1 = Utility::Split(str1, ":"); 
	vector 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->operator()(row1, row2, col1, col2); 
 
} 
 
 
template< class T > 
Vector Matrix::operator() (string str) 
{ 
	Vector temp = (Vector)*this; 
	return temp(str); 
} 
 
// //// 
// Arithmetic operators 
// //// 
 
template< class T > 
Matrix& Matrix::operator+= (Matrix& m) 
{ 
	return this->Add(m); 
} 
 
template< class T > 
Matrix& Matrix::operator-= (Matrix& m) 
{ 
	return this->Subtract(m); 
} 
 
template< class T > 
Matrix& Matrix::operator*= (Matrix& m) 
{ 
	return this->Multiply(m); 
} 
 
template< class T > 
Matrix& Matrix::operator/= (Matrix& m) 
{ 
	return this->Divide(m); 
} 
 
 
template< class T > 
Matrix& Matrix::operator+= (T v) 
{ 
	return this->Add(v); 
} 
 
template< class T > 
Matrix& Matrix::operator-= (T v) 
{ 
	return this->Subtract(v); 
} 
 
template< class T > 
Matrix& Matrix::operator*= (T v) 
{ 
	return this->Multiply(v); 
} 
 
template< class T > 
Matrix& Matrix::operator/= (T v) 
{ 
	return this->Divide(v); 
} 
 
 
// TYPE CONVERSIONS 
 
//template< class T > 
//Matrix::Matrix(Array &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[Columns()]; 
//		Utility::CheckPointer(columns); 
//		for(int i=0; i(data, columns); 
//		Utility::CheckPointer(clean); 
// 
//	} 
//	else if(m.ndims > 0) 
//	{ 
//		length = m.length; 
//		data = m.data; 
//		xDim = 1; 
//		yDim = length; 
//		 
//		columns = new (std::nothrow) Vector[1]; 
//		Utility::CheckPointer(columns); 
//		columns[0].Set(&(data[0]),length); 
//		clean = new (std::nothrow) Cleaner(data, columns); 
//		Utility::CheckPointer(clean); 
//		 
//		if(ndims != 1) 
//		{ 
//			Utility::Warning("Array (which was not of dimension 2) is converted to a Nx1 (column) Matrix. Dimensionality information is lost."); 
//		} 
//	} 
// 
//} 
 
 
template< class T > 
Matrix::Matrix(Vector &v)  
{ 
	ndims = 2; 
 
	length = v.length; 
	data = v.data; 
	xDim = 1; 
	yDim = length; 
	 
	columns = new (std::nothrow) Vector[1]; 
	Utility::CheckPointer(columns); 
	columns[0].Set(&(data[0]),length); 
	clean = new (std::nothrow) Cleaner(data, columns); 
	Utility::CheckPointer(clean); 
	 
} 
 
 
template< class T > 
Matrix::Matrix(SubMatrix &m)  
{ 
	if(m.yDim < 1 || m.xDim < 1) 
	{ 
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
		Utility::RunTimeError("All Matrix dimensions should be larger than 0!"); 
	} 
 
	ndims = 2; 
	length = m.yDim*m.xDim; 
	xDim = m.xDim; 
	yDim = m.yDim; 
	data = new (std::nothrow) T[length]; 
	Utility::CheckPointer(data); 
 
	columns = new (std::nothrow) Vector[xDim]; 
	Utility::CheckPointer(columns); 
	for(int i=0; i(data, columns); 
	Utility::CheckPointer(clean); 
 
	for(int i=0;i