matrixoperation1.cpp

该程序包实现了几个常用的模式识别分类器算法

// 实现矩阵操作
// 矩阵运算.cpp : Defines the entry point for the console application.
//
//#include "stdafx.h"
#include "stdlib.h"
#include <iostream>
using namespace std;
#include "global.h"

//构造函数
CMatrix::CMatrix()
{
	pdata=NULL;
	nline=0;
	ncol=0;
}
//析构函数
CMatrix::~CMatrix()
{
	if(pdata!=NULL)
	{
	   delete[] pdata;
	}
}

//输出矩阵
void CMatrix::shownumber()
{
	int i,j;
	for(i=0;i<nline;i++)
	{
	   for(j=0;j<ncol;j++)
	   {
		cout<<pdata[i*ncol+j]<<"   ";
	   }
	   cout<<endl;
	}
}
//方阵的值与方阵的幂还有逆的判断，行和列是否相等
states CMatrix::issamesize1()
{
	states flag;
	if(this->nline==this->ncol)
	{
	   flag=noerror;
	}
	else
	{
	   flag=nosamesize;
	}
	return flag;
}
//矩阵加法与矩阵减法的判断，第一个矩阵的行与第二个矩阵的行，第一个矩阵的列与第二个矩阵的列是否相等
states CMatrix::issamesize2(CMatrix &b)
{    
	states flag;
	if(this->nline==b.nline&&this->ncol==b.ncol)
	{   
	   flag=noerror;
	}
	else
	{
	   flag=nosamesize;
	}
	return flag;
}
//矩阵乘法的判断，第一个矩阵的列是否与第二个矩阵的行相等
states CMatrix::issamesize3(CMatrix &b)
{
	states flag;
	if(this->ncol==b.nline)
	{ 
	   flag=noerror;
	}
	else
	{
	   flag=nosamesize;
	}
	return flag;
}
//矩阵除法的判断，第一个矩阵的列与第二个矩阵的行，第二个矩阵的行与第二个矩阵的列是否相等 
states CMatrix:: issamesize4(CMatrix &b)
{   
	states flag;
	if(this->ncol==b.nline&&b.nline==b.ncol)
	{
	   flag=noerror;
	}
	else
	{
	   flag=nosamesize;
	}
	return flag;
}
//指针的判断
states CMatrix::isempty()
{
	states flag;
	if(this->pdata!=NULL)
	{
	   flag=noerror;
	}
	else
	{
	   flag=nopdataspace;
	}
	return flag;
}
//求方阵的值
DOUBLE CMatrix::det()
{  
	DOUBLE result;
	//如果是方阵
	if(this->issamesize1()==noerror&&this->isempty()==noerror)
	{
		this->flag_det = 1;
		if(this->nline==1)
		{
			result=this->pdata[0];
		}
		else if(this->nline==2)
		{
			result=this->pdata[0]*this->pdata[3]-this->pdata[1]*this->pdata[2];
		}
		else
		{
			int aa, bb;
			CMatrix mb;
			mb.nline=this->nline;
			mb.ncol=this->ncol;
			mb.pdata=new DOUBLE[mb.nline*mb.ncol];
			int totalNum = (this->nline)*(this->ncol);
			DOUBLE max = 0;
			DOUBLE x;
			int line = 0;
			int collumn = 0;
			DOUBLE sign = 1;
			//寻找全局绝对值最大的元素，记录行号和列号，然后交换到第0行，第0列
			for (aa=0; aa<this->nline; aa++)
			{
				for (bb=0; bb<this->ncol; bb++)
				{
					x = this->pdata[aa*this->ncol+bb];
					mb.pdata[aa*this->ncol+bb] = x;
					if (x>0 && x>max)
					{
						line = aa;
						collumn = bb;
						max = x;
					}
					else if (x<0 && (x+max)<0)
					{
						line = aa;
						collumn = bb;
						max = -x;
					}
				}
			}
			if (line!=0)
			{
				sign = sign * -1;
				for (bb=0; bb<this->ncol; bb++)
				{
					x=mb.pdata[line*mb.ncol+bb];
					mb.pdata[line*mb.ncol+bb] = mb.pdata[bb];
					mb.pdata[bb] = x;
				}
			}
			if (collumn!=0)
			{
				sign = sign * -1;
				for (aa=0; aa<this->nline; aa++)
				{
					x=mb.pdata[aa*mb.ncol+collumn];
					mb.pdata[aa*mb.ncol+collumn] = mb.pdata[aa*mb.ncol];
					mb.pdata[aa*mb.ncol] = x;
				}
			}
			//消零第一列元素，除了（0，0）位置的元素
			for (aa=1; aa<this->nline; aa++)
			{
				for (bb=1; bb<this->ncol; bb++)
				{
					mb.pdata[aa*this->ncol+bb] -= mb.pdata[bb] * mb.pdata[aa*this->ncol] / mb.pdata[0];
				}
			}
			CMatrix newM;
			newM.nline = this->nline - 1;
			newM.ncol = this->ncol - 1;
			newM.pdata = new DOUBLE[newM.nline * newM.ncol];
			for (aa=1; aa<this->nline; aa++)
			{
				for (bb=1; bb<this->ncol; bb++)
				{
					newM.pdata[(aa-1)*(newM.ncol)+bb-1] = mb.pdata[aa*(this->ncol)+bb];
				}
			}
			//迭代计算行列式
			result = sign * mb.pdata[0] * newM.det();
		}
	}
	//如果不是方阵
	else
	{ 
	   this->flag_det=0;
	   result=0;
	}
	return result;
}
//方阵的幂 
states CMatrix::pow(CMatrix &c,int k)
{
	int i, j, a, aa;
	states flag;
	//如果是方阵
	if(this->issamesize1()==noerror&&this->isempty()==noerror)
	{
	   flag=noerror;
	   c.nline=this->nline;
	   c.ncol=this->ncol;
	   c.pdata=new DOUBLE[c.nline*c.ncol];
	   CMatrix pp;
	   pp.pdata=new   DOUBLE[c.nline*c.ncol];
	   for(i=0;i<this->nline;i++)
	   {
		for(j=0;j<this->ncol;j++)
		{ 
		 pp.pdata[i*this->ncol+j]=this->pdata[i*this->ncol+j];
		}
	   }
	   //k大于1的时候参与循环
	   if(k>1)
	   {
		for(a=1;a<k;a++)
		{  
		 for(i=0;i<this->nline;i++)
		 {
		  for(j=0;j<this->ncol;j++)
		  {    
		   c.pdata[i*c.ncol+j]=0;
		   for(k=0;k<this->ncol;k++)
		   {  
			c.pdata[i*c.ncol+j]+=pp.pdata[i*this->ncol+k]*this->pdata[k*this->ncol+j];
		   }
		  }
		 }
		 for(i=0;i<this->nline;i++)
		 {
		  for(j=0;j<this->ncol;j++)
		  {
		   pp.pdata[i*this->ncol+j]=c.pdata[i*c.ncol+j];
		  }
		 }
		}
	   }
	   //k为1的时候直接赋值
	   else
	   {
		for(i=0;i<this->nline;i++)
		{
		 for(j=0;j<this->ncol;j++)
		 {
		  c.pdata[i*this->ncol+j]=this->pdata[i*this->ncol+j];
		 }
		}
	   }
	}
	//如果不是方阵
	else
	{
	   flag=error;
	}
	return flag;
}
//方阵的逆 
states CMatrix::inv(CMatrix &c) 
{
	int n, i, j, k;
	int indexI, indexJ;
	int nline = this->nline;
	int ncol = this->ncol;
	int size = nline*ncol;
	DOUBLE max, temp;
	states flag;
	flag=this->issamesize1();
	if(flag==noerror&&this->isempty()==noerror)
	{                                                
		DOUBLE x;
		x=this->det();
		//det不为0的时候参与计算
		if(x!=0)
		{   
			flag=noerror;
			//nline大于等于2的时候参与循环
			if(this->nline>=2)
			{
				CMatrix ma, ml, mr;
				ma.nline=nline;
				ml.nline=nline;
				mr.nline=nline;
				ma.ncol=ncol;
				ml.ncol=ncol;
				mr.ncol=ncol;
				ma.pdata=new DOUBLE[size];
				ml.pdata=new DOUBLE[size];
				mr.pdata=new DOUBLE[size];
				for(i=0;i<size;i++)
				{
					ma.pdata[i] = this->pdata[i];
					ml.pdata[i] = 0;
					mr.pdata[i] = 0;
				}
				for(i=0;i<nline;i++)
				{
					ml.pdata[i*nline+i] = 1;
					mr.pdata[i*nline+i] = 1;
				}
				for(n=0; n<nline-1; n++)
				{
					if (ma.pdata[n*nline+n]<0)
						max = -ma.pdata[n*nline+n];
					else
						max = ma.pdata[n*nline+n];
					indexI = n;
					indexJ = n;
					for(i=n; i<nline; i++)
					{
						for(j=n; j<nline; j++)
						{
							if (ma.pdata[i*nline+j]>max)
							{
								max = ma.pdata[i*nline+j];
								indexI = i;
								indexJ = j;
							}
							else if (ma.pdata[i*nline+j]<(-max))
							{
								max = -ma.pdata[i*nline+j];
								indexI = i;
								indexJ = j;
							}
						}
					}
					if (indexI != n)//行交换
					{
						for (j=n; j<nline; j++)
						{
							temp = ma.pdata[indexI*nline+j];
							ma.pdata[indexI*nline+j] = ma.pdata[n*nline+j];
							ma.pdata[n*nline+j] = temp;
						}
						for (j=0; j<nline; j++)
						{
							temp = ml.pdata[indexI*nline+j];
							ml.pdata[indexI*nline+j] = ml.pdata[n*nline+j];
							ml.pdata[n*nline+j] = temp;
						}
					}
					if (indexJ != n)//列交换
					{
						for (i=n; i<nline; i++)
						{
							temp = ma.pdata[i*nline+indexJ];
							ma.pdata[i*nline+indexJ] = ma.pdata[i*nline+n];
							ma.pdata[i*nline+n] = temp;
						}
						for (i=0; i<nline; i++)
						{
							temp = mr.pdata[i*nline+indexJ];
							mr.pdata[i*nline+indexJ] = mr.pdata[i*nline+n];
							mr.pdata[i*nline+n] = temp;
						}
					}
					for(i=n+1; i<nline; i++)//行变换
					{
						temp = ma.pdata[i*nline+n]/ma.pdata[n*nline+n];
						for(j=n+1; j<nline; j++)
						{
							ma.pdata[i*nline+j] -= ma.pdata[n*nline+j]*temp;
						}
						for(j=0; j<nline; j++)
						{
							ml.pdata[i*nline+j] -= ml.pdata[n*nline+j]*temp;
						}
					}
					for(j=n+1; j<nline; j++)//列变换
					{
						temp = ma.pdata[n*nline+j]/ma.pdata[n*nline+n];
						for(i=0; i<nline; i++)
						{
							mr.pdata[i*nline+j] -= mr.pdata[i*nline+n]*temp;
						}
					}
					for(i=0; i<nline; i++)//列变换
					{
						mr.pdata[i*nline+n] = mr.pdata[i*nline+n]/ma.pdata[n*nline+n];
					}
				}
				for(i=0; i<nline; i++)//列变换
				{
					mr.pdata[i*nline+nline-1] = mr.pdata[i*nline+nline-1]/ma.pdata[size-1];
				}
				mr.mul(ml, c);
			}
			else //nline为1的时候
			{  
				c.nline=1;
				c.ncol=1;
				c.pdata=new DOUBLE[1];
				c.pdata[0]=1/x;
			}
		}
		else//det为0的时候直接输出错误信息
		{
			flag=error;
			cout<<endl;
			cout<<"the matrix'det is zero"<<endl;
		}
	}
	else//如果不是方阵
	{
		flag=error;
	}
	return flag;
}
//矩阵的转置 
states CMatrix::tran(CMatrix &c)
{
	int i,j;
	states flag;
	//如果矩阵不为空
	if(this->isempty()==noerror)
	{
	   flag=noerror;
	   c.nline=this->ncol;
	   c.ncol=this->nline;
	   c.pdata=new DOUBLE[c.nline*c.ncol];
	   for(i=0;i<c.nline;i++)
	   {
		for(j=0;j<c.ncol;j++)
		{
		 c.pdata[i*c.ncol+j]=this->pdata[j*this->ncol+i];
		}
	   }
	}
	//如果矩阵为空
	else
	{
	   flag=error;
	}
	return flag;
}                                                         
//矩阵的加法
states CMatrix::add(CMatrix &b, CMatrix &c)
{
	int i,j;
	states flag;
	//如果两矩阵可以参加加法运算
	if(this->issamesize2(b)==noerror&&this->isempty()==noerror&&b.isempty()==noerror)
	{
	   flag=noerror;
	   if(c.pdata!=NULL)
	   {
		delete[] c.pdata;
	   }
	   c.nline=this->nline;
	   c.ncol=b.ncol;
	   c.pdata=new DOUBLE[c.nline*c.ncol];
	   for(i=0;i<c.nline*c.ncol;i++)
	   {
		c.pdata[i]=this->pdata[i]+b.pdata[i];
	   }
	}
	//如果两矩阵不能参加加法运算
	else
	{
	   flag=error;
	}
	return flag;
}
//矩阵的减法
states CMatrix::sub(CMatrix &b, CMatrix &c)
{
	int i,j;
	states flag;
	//如果两矩阵可以参加减法运算
	if(this->issamesize2(b)==noerror&&this->isempty()==noerror&&b.isempty()==noerror)
	{
	   flag=noerror;
	   if(c.pdata!=NULL)
	   {
		delete[] c.pdata;
	   }
	   c.nline=this->nline;
	   c.ncol=b.ncol;
	   c.pdata=new DOUBLE[c.nline*c.ncol];
	   for(i=0;i<c.nline*c.ncol;i++)
	   {
		c.pdata[i]=this->pdata[i]-b.pdata[i];
	   }
	}
	//如果两矩阵不能参加减法运算
	else
	{
	   flag=error;
	}
	return flag;
}
//矩阵的乘法
states   CMatrix::mul(CMatrix &b, CMatrix &c)
{
	int i,j,k;
	states flag;
	//如果两矩阵可以参加乘法运算
	if(this->issamesize3(b)==noerror && this->isempty()==noerror && b.isempty()==noerror)
	{
	   flag=noerror;
	   if(c.pdata!=NULL)
	   {
		delete[] c.pdata;
	   }
	   c.nline=this->nline;
	   c.ncol=b.ncol;
	   c.pdata=new DOUBLE[c.nline*c.ncol];
	   for(i=0;i<c.nline;i++)
	   {
		for(j=0;j<c.ncol;j++)
		{   
		 c.pdata[i*c.ncol+j]=0;
		 for(k=0;k<this->ncol;k++)
		 {
		  c.pdata[i*c.ncol+j]+=this->pdata[i*this->ncol+k]*b.pdata[k*b.ncol+j];
		 }
		}
	   }
	}
	//如果两矩阵不能参加减法运算
	else
	{
	   flag=error;
	}
	return flag;
}     
//矩阵的除法 
states CMatrix::divi(CMatrix &b, CMatrix &c)
{
	int i,j;
	states flag;
	//如果两矩阵可以参加除法运算
	if(this->issamesize4(b)==noerror&&this->isempty()==noerror&&b.isempty()==noerror)
	{   
	   flag=noerror;
	   c.nline=b.nline;
	   c.ncol=b.ncol;
	   c.pdata=new DOUBLE[c.nline*c.ncol];
	   for(i=0;i<b.nline;i++)
	   {
		for(j=0;j<b.ncol;j++)
		{  
		 c.pdata[i*c.ncol+j]=b.pdata[i*c.ncol+j];
		}
	   }
	   states flag1;
	   //调用求逆运算
	   flag1=b.inv(c);
	   if(flag==noerror)
	   {
		for(i=0;i<b.nline;i++)
		{
		 for(j=0;j<b.ncol;j++)
		 {  
		  b.pdata[i*b.ncol+j]=c.pdata[i*b.ncol+j];
		 }
		}
	   }
	   //调用方阵的逆
	   flag1=this->mul(b,c);
	}
	//如果两矩阵不能参加减法运算
		 else
	{
	   flag=error;
	}
	return flag;
}                                                              

/*
matrix operation examples:
		//方阵的行列式
	   case '1':
		result=mymatrix1.det();
		if(mymatrix1.flag_det==1)
		{   
		 cout<<endl;
		 cout<<"   det(mymatrix1)=   "<<result<<endl;; 
		}
		else
		{  
		   cout<<endl;
		   cout<<"   error in det."<<endl;
		   mymatrix1.flag_det=1;
		}
		break;
		//方阵的幂 
	   case '2':
		cout<<"input   your kk:"<<endl;
		cin>>kk;
		flag=mymatrix1.pow(mymatrix3,kk);
		break;
		//方阵的逆
	   case '3':
		flag=mymatrix1.inv(mymatrix3);
		break;
		//矩阵的转置
	   case '4':
		flag=mymatrix1.tran(mymatrix3);
		break;
		//矩阵的加法
	   case '5':
		mymatrix2.getnumber();
		flag=mymatrix1.add(mymatrix2,mymatrix3);
		break;
		//矩阵的减法
	   case '6':
		mymatrix2.getnumber();
		flag=mymatrix1.sub(mymatrix2,mymatrix3);
		break;
		//矩阵的乘法
	   case '7':
		mymatrix2.getnumber();
		flag=mymatrix1.mul(mymatrix2,mymatrix3);
		break;
		//矩阵的除法，为了方便在这里我们只讨论被除矩阵是方阵的情况
	   case '8':
		mymatrix2.getnumber();
		flag=mymatrix1.divi(mymatrix2,mymatrix3);
		break;
*/