my_matrix.cpp

DEV C++ 写的一个矩阵类

#include <iostream>
#include <vector>
#include <stdexcept>
#include <iterator>
#include "my_matrix.h"
using namespace std;
 
typedef vector<vector<double> >::size_type row;
typedef vector<double>::size_type out;

matrix::matrix(int i,int j)
{
	for(int a=0;a<i;++a)
		vec_vec.push_back(vector<double>(j,0));
}

matrix matrix::operator*(matrix mat)const
{
	if(this->empty()||mat.empty())
	{
		throw runtime_error("Empty matrix cann't multiply.");
	}
	row row1=vec_vec.size();
	row row2=mat.vec_vec.size();
	out out1=vec_vec[0].size();
	out out2=mat.vec_vec[0].size();
	if(out1!=row2)
	{
		throw runtime_error("They cann't multiply.");
	}
	matrix ret;
	matrix mat_=mat.trans();
	for(int i=0;i<row1;++i)
	{
		const vector<double>& vec1=vec_vec[i];
		vector<double> vec_ret;
		for(int j=0;j<out2;++j)
		{
			const vector<double>& vec2=mat_.vec_vec[j];
			double temp=0;			
			for(int k=0;k<out1;++k)
			{
				temp+=vec1[k]*vec2[k];
			}
			vec_ret.push_back(temp);
		}
		ret.vec_vec.push_back(vec_ret);
	}
	return ret; 
}

matrix matrix::operator*(double d)const
{
	if(empty())
	{
		throw runtime_error("Empty matrix cann't add.");
	}
	row rows=vec_vec.size();
	out outs=vec_vec.size();
	
	vector<vector<double> > ret;
	for(int i=0;i<rows;++i)
	{		
		vector<double> vec;
		vec.push_back(vec_vec[i][0]*d);
		ret.push_back(vec);
	}
	for(int i=0;i<rows;++i)
	{
		vector<double>& vec=ret[i];
		for(int j=1;j<outs;++j)
			vec.push_back(vec_vec[i][j]*d);
	}
	
	return matrix(ret);
}

matrix operator*(double d,const matrix& mat)
{
	return mat*d;
}

matrix matrix::operator+(const matrix& mat)const
{
	if(empty()||mat.empty())
	{
		throw runtime_error("Empty matrix cann't add.");
	}
	row row1=vec_vec.size();
	row row2=mat.vec_vec.size();
	out out1=vec_vec[0].size();
	out out2=mat.vec_vec[0].size();
	if(row1!=row2||out1!=out2)
	{
		throw runtime_error("They cann't add.");
	}
	matrix ret;
	for(int i=0;i<row1;++i)
	{
		const vector<double>& vec1=vec_vec[i];	
		const vector<double>& vec2=mat.vec_vec[i];
		vector<double> vec_ret;	
		for(int j=0;j<out1;++j)
		{
			vec_ret.push_back(vec1[j]+vec2[j]);
		}
		ret.vec_vec.push_back(vec_ret);
	}
	return ret;
}	

matrix& matrix::operator+=(const matrix& mat)
{
	*this=*this+mat;
	return *this;
}	

matrix matrix::operator-(const matrix& mat)const
{
	if(empty()||mat.empty())
	{
		throw runtime_error("Empty matrix cann't submiss.");
	}
	row row1=vec_vec.size();
	row row2=mat.vec_vec.size();
	out out1=vec_vec[0].size();
	out out2=mat.vec_vec[0].size();
	if(row1!=row2||out1!=out2)
	{
		throw runtime_error("They cann't submiss.");
	}
	matrix ret;
	for(int i=0;i<row1;++i)
	{
		const vector<double>& vec1=vec_vec[i];	
		const vector<double>& vec2=mat.vec_vec[i];
		vector<double> vec_ret;	
		for(int j=0;j<out1;++j)
		{
			vec_ret.push_back(vec1[j]-vec2[j]);
		}
		ret.vec_vec.push_back(vec_ret);
	}
	return ret;
}	

matrix& matrix::operator-=(const matrix& mat)
{
	*this=*this-mat;
	return *this;
}

matrix matrix::operator-()const
{
	if(empty())
	{
		throw runtime_error("Empty matrix cann't have -.");
	}
	vector<vector<double> > ret;
	for(int i=0;i<vec_vec.size();++i)
	{
		vector<double> vec;
		for(int j=0;j<vec_vec[0].size();++j)
		{
			vec.push_back(-vec_vec[i][j]);
		}
		ret.push_back(vec);
	}
	return ret;
}

matrix matrix::elim(int& count)//rows==outs;
{
	if(empty())
	{
		throw runtime_error("Empty matrix cann't eliminate.");
	}
	row rows=vec_vec.size();
	out outs=vec_vec[0].size();		
	for(int i=0;i<rows;++i)
	{		
		if(maxrow(i)!=i)
		{
			changeline(maxrow(i),i);
			++count;
		}
		const vector<double>& vec=vec_vec[i];
		for(int j=i+1;j<rows;++j)
		{
			vector<double>& vec_=vec_vec[j];
			double a=vec_vec[j][i];
			for(int k=i;k<outs;++k)
			{
				if(vec[i]==0)
					vec_[k]=0;
				else
					vec_[k]=vec_[k]-a*vec[k]/vec[i];
			}
		}
	}
	
	return *this;
}

double fab(double d)
{
	if(d<0)
		return -d; 
	else
		return d;
}

int matrix::maxrow(int i)
{
	vector<vector<double> >::size_type rows=vec_vec.size();
	int a=i;
	for(int j=i+1;j<rows;++j)
	{
		if(fab(vec_vec[j][a])>fab(vec_vec[i][a]))
			i=j;
	}
	return i;
}	
	
double matrix::mold()
{	
	if(empty())
	{
		throw runtime_error("Empty matrix cann't eliminate.");
	}
	row rows=vec_vec.size();
	out outs=vec_vec[0].size();
	if(rows!=outs)
	{
		throw runtime_error("The matrix has no mold.");
	}
	
	if(rows==1&&outs==1)
		return vec_vec[0][0];
		
	int count=0;	
	matrix ret(elim(count));
	double temp=1;		
	for(int i=0;i<rows;++i)
		temp=temp*ret.vec_vec[i][i];
	if(fab(temp)<1e-12)
		temp=0;
	if(count%2==0)
		return temp;
	else 
		return -temp;
}

matrix matrix::inver()
{
	if(empty())
	{
		throw runtime_error("Empty matrix cann't inverse.");
	}
	row rows=vec_vec.size();
	out outs=vec_vec[0].size();	
	matrix mat(*this);
	double m=mat.mold();
	if(rows!=outs||m==0)
	{
		throw runtime_error("The matrix cann't inverse.");
	}	
	if(rows==1&&outs==1)
	{
		mat.vec_vec[0][0]=1/mat.vec_vec[0][0];
		return mat;
	}
	for(int i=0;i<rows;++i)
	{
		matrix ret(*this);
		vector<vector<double> >::iterator vec_vecit=ret.vec_vec.begin();		
		ret.vec_vec.erase(vec_vecit+i);
		matrix ret_=ret;
		for(int j=0;j<outs;++j)
		{		
			ret=ret_;
			for(int k=0;k<rows-1;++k)
			{
				vector<double>::iterator vecit=ret.vec_vec[k].begin();
				ret.vec_vec[k].erase(vecit+j);						
			}
			double m_=ret.mold();
			int a;
			if((i+j)%2==0)
				a=1;
			else 
				a=-1;				
			mat.vec_vec[j][i]=a*m_/m;
		}
	}
	return mat;
}

matrix matrix::trans()
{
	if(this->empty())
	{
		throw runtime_error("Empty matrix cann't transform.");
	}
	row row1=vec_vec.size();
	out out1=vec_vec[0].size();
	vector<vector<double> > vec_v;
	const vector<double>& vec_ret=vec_vec[0];
	for(int j=0;j<out1;++j)
	{
		vector<double> vec1;
		vec1.push_back(vec_ret[j]);
		vec_v.push_back(vec1);
	}	
	for(int i=1;i<row1;++i)
	{	
		const vector<double>& vec_ret=vec_vec[i];
		for(int j=0;j<out1;++j)
		{
			vector<double>& vec1=vec_v[j];
			vec1.push_back(vec_ret[j]);
		}
	}
	return matrix(vec_v);
}

istream& operator>>(istream& in,matrix& mat)
{
	if(!mat.empty())
		mat=matrix();
	in.clear();
	double temp;
	int row;
	int out;
	in>>row>>out;
	for(int i=0;i<row;++i)
	{
		vector<double> vec;
		for(int j=0;j<out;++j)
		{
			in>>temp;
			vec.push_back(temp);
		}
		mat.vec_vec.push_back(vec);
	}
	if(!in)
	{
		mat=matrix(); 
		throw runtime_error("Failed to assign the matrix.");
	}	
	return in;
}

ostream& operator<<(ostream& os,const matrix& mat)
{
	if(mat.empty())
	{
		throw runtime_error("Empty matrix cann't output.");
	}
	vector<vector<double> >::size_type row=mat.vec_vec.size();
	vector<double>::size_type out=mat.vec_vec[0].size();
	for(int i=0;i<row;++i)
	{
		const vector<double>& vec=mat.vec_vec[i];
		for(int j=0;j<out;++j)
		{
			os.width(10);
			if(fab(vec[j])<1e-12)
				os<<0<<"   ";
			else
				os<<vec[j]<<"   ";
		}
		os<<endl;
	}
	return os;
}