matrixs.cpp

Moravec算子实现程序 Moravec算子实现程序

// Matrixs.cpp: implementation of the Matrixs class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Dan_backcom.h"
#include "Matrixs.h"
#include "iostream.h"
#include "math.h"

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

Matrixs::Matrixs(int mr,int mc)
{
	rows=mr;
    cols=mc;
	elems=new double[(rows+1)*(cols+1)];
}

Matrixs::~Matrixs()
{
	delete []elems;
}

void Matrixs::Initelems()               //为elems分配内存空间
{
	elems=new double[(rows+1)*(cols+1)];
}

double& Matrixs::operator () (int i,int j) //矩阵的表示
{
	return elems[(i-1)*cols+j-1];
}


Matrixs::Matrixs(Matrixs& ma)
{
	rows=ma.rows;
	cols=ma.cols;
	elems=new double[(rows+1)*(cols+1)];
	for(int i=1;i<=rows;i++)
	  for(int j=1;j<=cols;j++)
	   (*this)(i,j)=ma(i,j);
}



Matrixs& Matrixs::operator=(Matrixs& ma)  //矩阵的相等
{
	cols=ma.cols;
	rows=ma.rows;
	elems=new double[(rows+1)*(cols+1)];
	for(int i=1;i<=rows;i++)
	  for(int j=1;j<=cols;j++)
	   (*this)(i,j)=ma(i,j);
	return *this;
}

Matrixs Matrixs::operator+(Matrixs ma)  //矩阵的加法
{
	Matrixs c(rows,cols);
	for(int i=1;i<=rows;i++)
	  for(int j=1;j<=cols;j++)
	   c(i,j)=(*this)(i,j)+ma(i,j);
	return c;
}


Matrixs Matrixs::operator-(Matrixs ma)  //矩阵的减法
{
	Matrixs c(rows,cols);
	for(int i=1;i<=rows;i++)
	  for(int j=1;j<=cols;j++)
	   c(i,j)=(*this)(i,j)-ma(i,j);
	return c;
}


Matrixs Matrixs::operator * (Matrixs ma)   //矩阵的乘法
{
	int i,j,k;
	if(cols!=ma.rows)
	{
		cout<<"所给矩阵不能相乘！"<<endl;
		return *this;
	}
	else
	{
		Matrixs c(rows,ma.cols);
	for(i=1;i<=rows;i++)
		for(j=1;j<=ma.cols;j++)
		{
			c(i,j) = 0;
			for(k=1;k<=cols;k++)
				c(i,j) += (*this)(i,k)*ma(k,j);
		}
		return c;
	}
}


Matrixs Matrixs::operator * (double x)   //矩阵的数乘
{
	Matrixs c(rows,cols);
	for(int k=1;k<=rows;k++)
		for(int t=1;t<=cols;t++)
		c(k,t) = x * (*this)(k,t);
	return c;
}


void Matrixs::print()				//打印矩阵
{
    for(int r=1;r<=this->rows;++r)
    {
        for(int c=1;c<=this->cols;c++)
            cout<<(*this)(r,c)<<"  ";
        cout<<endl;
    }
}

Matrixs Matrixs::T()		/*矩阵转置*/
{
	Matrixs c(cols,rows);
	for(int i=1;i<=rows;i++)
		for(int j=1;j<=cols;j++)
			c(j,i) = (*this)(i,j);
		return c;
}


double Matrixs::Surplus() /*求矩阵行列式*/ 
{ 
	double x=1;
	if(rows!=cols)
	{
		cout<<"矩阵不为方阵，无行列式！"<<endl;
		return 0;
	}
	else
	{
		Matrixs c(rows,cols);
		c = (*this);
		for(int i=1;i<cols;i++)
		{
			for(int j=i+1;j<=rows;j++)
			{
				if(c(i,i))
				for(int k=cols;k>=i;k--)
				{
					c(j,k)-=(c(j,i)/c(i,i)) * c(i,k);
				}
				else if(!c(i,i) && c(j,i))
				{
					double temp;
					for(int m=i;m<=cols;m++)
					{
						temp = c(i,m);
						c(i,m) = c(j,m);
						c(j,m) = temp;
					}
				}
			}
		}
		for(int m=1;m<=rows;m++)
			x*=c(m,m);
		return x;
	}
}



Matrixs Matrixs::MatrixOpp() /*矩阵求逆*/ 
{ 
	if(rows!=cols || (*this).Surplus()==0)
	{
		cout<<"求逆矩阵错误！"<<endl;
	}
    int i,j,m,n; 
    double X; 
	Matrixs C(rows,cols);
    Matrixs SP(rows,cols); 
	Matrixs AB(rows,cols);
	Matrixs A(rows,cols);
      
     
    X=1/Surplus(); 
         
    for(i=1;i<=rows;i++) 
    for(j=1;j<=cols;j++) 
    {
		A = (*this); 		
		for(m=1;m<=cols;m++) 
		A(i,m)=0; 
		for(n=1;n<=rows;n++) 
		A(n,j)=0; 
		A(i,j)=1; 
		SP(j,i)=(1 - 2*((i+j)%2)) * A.Surplus(); 
		
    } 
     C=SP * X; 
     
     return C; 
}


int Matrixs::invers_matrix(int n)//矩阵求逆
{
	double *m1;
	m1 = this->elems;
	int *is,*js;
	int i,j,k,l,u,v;
	double temp,max_v;
	is=(int *)malloc(n*sizeof(int));
	js=(int *)malloc(n*sizeof(int));
	if(is==NULL||js==NULL){
		printf("out of memory!\n");
		return(0);
	}
	for(k=0;k<n;k++){
		max_v=0.0;
		for(i=k;i<n;i++)
			for(j=k;j<n;j++){
				temp=fabs(m1[i*n+j]);
				if(temp>max_v){
					max_v=temp; is[k]=i; js[k]=j;
				}
			}
			if(max_v==0.0){
				free(is); free(js);
				printf("invers is not availble!\n");
				return(0);
			}
			if(is[k]!=k)
				for(j=0;j<n;j++){
					u=k*n+j; v=is[k]*n+j;
					temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
				}
				if(js[k]!=k)
					for(i=0;i<n;i++){
						u=i*n+k; v=i*n+js[k];
						temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
					}
					l=k*n+k;
					m1[l]=1.0/m1[l];
					for(j=0;j<n;j++)
						if(j!=k){
							u=k*n+j;
							m1[u]*=m1[l];
						}
						for(i=0;i<n;i++)
							if(i!=k)
								for(j=0;j<n;j++)
									if(j!=k){
										u=i*n+j;
										m1[u]-=m1[i*n+k]*m1[k*n+j];
									}
									for(i=0;i<n;i++)
										if(i!=k){
											u=i*n+k;
											m1[u]*=-m1[l];
										}
	}
	for(k=n-1;k>=0;k--){
		if(js[k]!=k)
			for(j=0;j<n;j++){
				u=k*n+j; v=js[k]*n+j;
				temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
			}
			if(is[k]!=k)
				for(i=0;i<n;i++){
					u=i*n+k; v=i*n+is[k];
					temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
				}
	}
	free(is); free(js);
	return(1);
}


void Matrixs::mult(double *m2,double *result,int i_1,int j_12,int j_2)//矩阵相乘
{ 
	double *m1;
	m1 = this->elems;
	int i,j,k;
	for(i=0;i<i_1;i++)
		for(j=0;j<j_2;j++)
		{
			result[i*j_2+j]=0.0;
			for(k=0;k<j_12;k++)
				result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];
		}
		return;
}


void Matrixs::transpose(double *m2,int m,int n)  //矩阵转置          
{
	double *m1;
	m1 = this->elems;
	int i,j;                                        
	for(i=0;i<m;i++)                                
		for(j=0;j<n;j++)                        
			m2[j*m+i]=m1[i*n+j];            
		return;                                 		
}