mat.cpp

矩阵类,自己编写,通过较多考验,符合matlab编程习惯,可处理无穷维矩阵和向量

// Mat.cpp: implementation of the Mat class.
// 矩阵运算 2006年1月由xiaofcen编制
// 所有涉及矩阵运算的首先必须初始化，即要调用 zero函数
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "LVB.h"
#include "mat.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

mat::mat()
{
	row=1;
	col=1;
	p=NULL;
}

mat::~mat()
{
	Freethis();
}

mat::mat(const mat& other)
{
	row = other.row;
	col = other.col;
	p = NULL;
	zero(row,col);
	memcpy(p, other.p, sizeof(double)*row*col);
}


void mat::zero(int m, int n)
{
	Freethis();
	row=m;	col=n;
	int nSize = row*col;
	if (nSize <= 0)
		ASSERT(true);
	p = new double[nSize];	
	if (p == NULL)
		exit(0);				        	// 内存分配失败
	if (IsBadReadPtr(p, sizeof(double) * nSize))
		exit(0);	
	memset(p,0, sizeof(double) * nSize);                         // 将各元素值置0
}

void mat::Freethis()
{
	if (p)
	{
		delete[] p;
		p = NULL;
	}
}

//------------显示某个元素-------------
double mat::r(int nRow, int nCol) const
{
	if (nCol <=0 || nCol > col || nRow <= 0 || nRow > row)
	{
		return 0.0;
		ASSERT(true);	
	}
	else
	{
		if (p == NULL)
		{
			return 0.0;
			ASSERT(true);
		}
		else
		{
			return p[ (nCol-1) +(nRow-1)* col];
		}
	}
}

//------------设置某个元素-------------
void mat::set(int nRow, int nCol, double value)
{
	if (nCol <= 0 || nCol > col || nRow <= 0 || nRow > row)
		exit(0);					      
	if (p == NULL)
		exit(0);						   
	p[(nCol-1) +(nRow-1)* col] = value;
}

//------------显示某一行-------------
mat mat::cpr(int m)
{
	mat out;out.zero(1,col);
	if( (m <= 0) || (m > row) )
		exit(0);
	for (int i=1;i<=col;i++)
	{
		out.set(1,i,this->r(m,i));
	}
	return out;
}

//------------显示某一列-------------
mat mat::cpc(int n)
{
	mat out;out.zero(row,1);
	if( (n <= 0) || ( n> col) )
		exit(0);
	for (int i=1;i<=row;i++)
	{
		out.set(i,1,this->r(i,n));
	}
	return out;
}

//------------设置某行元素-------------
void mat::setr(int nRow,const mat& b)
{
	if (nRow <= 0 || nRow > row)
		exit(0);					         
	if (p == NULL)
		exit(0);
	for (int i=1;i<=col;i++)
	{
		set(nRow,i,b.r(1,i));
	}
}

//------------设置某列元素-------------
void mat::setc(int m,const mat& b)
{
	if (m <=0 || m > col)
		exit(0);					         
	if (p == NULL)
		exit(0);
	for (int i=1;i<=row;i++)
	{
		set(i,m,b.r(i,1));
	}
}

//-------------转置-------------------
mat mat::T() const
{
	mat out; out.zero(col,row);
	for (int i=1;i<=row;i++)
	{
		for (int j=1;j<=col;j++)
		{
			out.set(i,j,this->r(j,i));			
		}
	}
	return out;
}

//-------------单位阵----------------
void mat::eyemat(int n)
{
	zero(n,n);
	for (int i=1;i<=n;i++)
	{ 
		for (int j=1;j<=n;j++)
		{
			if (i==j)
				set(i,j,1.0);				
			else
				set(i,j,0.0);
		}
	}
}

mat& mat::operator=(const mat& other)
{
	if (&other != this)
	{
		zero(other.row,other.col);
		memcpy(p, other.p, sizeof(double)*row*col);
	}
	return *this;
}

mat mat::operator*(double value) const
{
	mat	result(*this) ;	
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			 result.set(i,j,result.r(i,j)*value);
	}
	return result ;
}

mat mat::operator*(int value) const
{
	mat	result(*this) ;	
	for (int i = 1 ; i <= row ; i++)
	{
		for (int j = 1 ; j <= col; j++)
			 result.set(i,j,result.r(i,j)*value);
	}
	return result ;
}

mat	mat::operator*(const mat& other) const
{
	// 首先检查行列数是否符合要求
	ASSERT (col == other.row);
	mat result;result.zero(row,other.col);
	double	value ;
	for (int i = 1 ; i <=row ; i++)
	{
		for (int j = 1 ; j <= other.col ; j++)
		{
			value = 0.0 ;
			for (int k = 1 ; k <= col ; k++)
			{
				value += this->r(i, k) * other.r(k, j);
			}
			result.set(i,j,value) ;
		}
	}
	return result ;
}


mat	mat::operator+(const mat& other) const
{
	ASSERT (row == other.row && col == other.col); // 首先检查行列数是否相等
	mat	result;result.zero(row,col);
	for (int i = 1; i <= row ; i++)
	{
		for (int j = 1; j <= col; j++)
		{
			result.set(i,j,r(i,j)+other.r(i,j));
		}
	}
	return result ;
}

mat	mat::operator-(const mat& other) const
{
	ASSERT (row == other.row && col == other.col); // 首先检查行列数是否相等
	mat	result;result.zero(row,col);
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			result.set(i,j,r(i,j)-other.r(i,j));
	}
	return result;
}

mat mat::operator/(double a) const
{
	mat	result;result.zero(row,col);	
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			 result.set(i,j,this->r(i,j)/a);
	}
	return result ;
}

mat mat::operator/(int a) const
{
	mat	result;result.zero(row,col);
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			 result.set(i,j,this->r(i,j)/a);
	}
	return result ;
}

//----------------矢量除-----------------
mat mat::operator/(const mat& other) const
{
	mat out;
	int m=other.row;
	int n=other.col;
	if (m !=n)
	{
		exit(0);
	}
	else
	{
		mat thismat;thismat.zero(row,col);
		for (int i=1;i<=row;i++)
		{
			for (int j=1;j<=col;j++)
			{
				thismat.set(i,j,r(i,j));
			}
		}
		mat inva=other.invmat();
		out=thismat*inva;
	}
	return out;
}

//-----------------------矢量差乘------------------
mat mat::operator^(const mat& c) const
{
	mat out;	
	int m=c.row;
	int n=c.col;
	if ( (m == row ) && ( n == col)  && (m*n == 3) )
	{
		if (m =3)
		{
			out.zero(3,1);
			out.set(1,1,r(2,1)*c.r(3,1)-r(3,1)*c.r(2,1));
			out.set(2,1,r(3,1)*c.r(1,1)-r(1,1)*c.r(3,1));
			out.set(3,1,r(1,1)*c.r(2,1)-r(2,1)*c.r(1,1));
		}
		else
		{
			out.zero(1,3);
			out.set(1,1,r(1,2)*c.r(1,3)-r(1,3)*c.r(1,2));
			out.set(2,1,r(1,3)*c.r(1,1)-r(1,1)*c.r(1,3));
			out.set(3,1,r(1,1)*c.r(1,2)-r(1,2)*c.r(1,1));
		}
	}
	else
	{
		exit(0);
	}
	return out;
}

//---------------矢量点乘--------------------
double mat::operator|(const mat& a) const
{
	//无论是矩阵还是向量，只对第一列有效/列
	double out=0.0;
	int m=a.row;
	int n=a.col;
	if (( m==row ) && (n ==col) )
	{
		double temp=0.0;
		if (m==1)
		{
			for (int i=1;i<=n;i++)
			{
				temp+=a.r(1,i)*r(1,i);
			}
			out=temp;
		}
		else
		{
			if (n==1)
			{
				for (int i=1;i<=m;i++)
				{
					temp+=a.r(i,1)*r(i,1);
				}
				out=temp;
			}	
		}
	}
	return out;
}

//--------------- 2范数 -----------------------
double mat::normmat() const
{
	double out=0.0;
	int m=row;
	int n=col;
	if ((m !=1) && (n !=1))
	{
		exit(0);
	}
	if (m==1)
	{
		double temp=0.0;
		for (int i=1;i<=n;i++)
		{
			temp+=r(1,i)*r(1,i);
		}
		out=sqrt(temp);
	}
	else
	{
		if (n==1)
		{
			double temp=0.0;
			for (int i=1;i<=m;i++)
			{
				temp+=r(i,1)*r(i,1);
			}
			out=sqrt(temp);
		}
	}
	return out;
}


//---------------伴随阵-------------------------
mat mat::adj() const
{
	mat out;out.zero(row,col);
	for(int i=1;i<=row;i++)
	{
		for(int j=1;j<=col;j++)
		{
			if (((i+j)%2)!=0)
			{
				out.set(i,j,(-1.0)*(sub_mat(i,j).detmat()));
			}
			else
			{
				out.set(i,j,       (sub_mat(i,j).detmat()));
			}
		}
	}
	mat reuslt=out.T();
	return reuslt;
}

//---------------行列式的值-----------------------
double mat::detmat() const
{
	double out=0.0;
	if ( (col==1) &&(row==1) )
	{
		out=r(1,1);
	}
	else
	{
		if ((col==2)&&(row==2))
		{
			out=r(1,1)*r(2,2)-r(1,2)*r(2,1);
		}
		else
		{
			int i=1;
			for(int j=1;j<=col;j++)
			{
				if ( ((i+j)%2 )!=0)
				{
					out+=(-1.0)*sub_mat(i,j).detmat()*r(i,j);
				}
				else
				{
					out+=       sub_mat(i,j).detmat()*r(i,j);
				}
			}
		}
	}
	return out;
}

//--------------------余子式------------------------
mat mat::sub_mat(int rr,int cc) const
{
	mat out;out.zero(row-1,col-1);