////////////////////////////////////////////////////////////////
//文件名: Matrix.h
//作者: 祁大江
//功能：矩陣類：可運行矩陣間的加、減、乘以及矩陣與數的點乘；可計算矩陣的轉置、秩、行列式、逆；
//		其中求解矩陣逆用到了經典的高斯算法；求解行列式以及秩運用的是消元法。
//說明：本程序中數組的運算，都將二維數組轉換為一維數組來進行計算的，方便傳遞參數；程序中的矩陣都是方陣。
////////////////////////////////////////////////////////////////
#define N 9  //定義矩陣元素的個數
#include<iostream.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include"matrix.h"
Matrix::Matrix()
{
	list=new float[N];
	nElem=0;
}
void Matrix::Destroy()
{
	delete []list;
}
int Matrix::Elem(void)
{
	return nElem;
}
int Matrix::Elem(float elem)  //給矩陣元素賦值
{
	if(nElem<N)
	{
		list[nElem++]=elem;
		return nElem;
	}
	else
	{
		return 0;
	}
}
void Matrix::Print(void)
{
	for(int i=0;i<N;i++) 
	{
		cout<<list[i]<<'\t';
		if(((i+1)%int(sqrt(N)))==0)
		{
			cout<<endl;
		}
	}
	cout<<endl;
}
float Matrix::GetElem(int i) //返回矩陣中所要求的元素
{
	if((i>=0)&&(i<nElem))
	{
		return list[i];
	}
	else
	{
		return 0;
	}
}
/*int Matrix::zhi()
{
	int i,j,k;
	int t=0;
	float m;
	float temp[N*N]={0};
	int s=int(sqrt(N));
	for(i=0;i<sqrt(N);i++)
	{
		for(k=i+1;k<sqrt(N);k++)
		{
			for(j=0;j<sqrt(N);j++)
			{
				temp[t*int(sqrt(N))+j]=list[i*int(sqrt(N))+j]+list[k*int(sqrt(N))+j];
			}
			t++;
		}
	}
	cout<<"等比矩陣為:"<<endl;
	for(i=0;i<(N*(sqrt(N)-1))/2;i++) 
	{
		cout<<temp[i]<<'\t';
		if(((i+1)%int(sqrt(N)))==0)
		{
			cout<<endl;
		}
	}
	cout<<endl;
	for(i=0;i<(sqrt(N)*(sqrt(N)-1))/2;i++)
	{
		for(k=0;k<sqrt(N);k++)
		{
			m=temp[i*int(sqrt(N))]/list[k*int(sqrt(N))];
			for(j=1;j<sqrt(N);j++)
			{
				if((temp[i*int(sqrt(N))+j]/list[k*int(sqrt(N))+j])!=m)
				{
					break;
				}
			}
			if(j==sqrt(N))
			{
				s--;
			}
		}
	}
	return s;
}*/
int Matrix::zhi()//求矩陣的秩，沒有找到經典的算法
{
	int i,j,k,t;
	float temp;
	int flg=0;
	int f;
	int s=0;
	if(list[0]==0)//確保迭代的第一個元素不為0
	{
		for(i=1;i<sqrt(N);i++) //行
		{
			for(k=0;k<sqrt(N);k++) //列
			{
				if(list[i*int(sqrt(N))+k]!=0) //找到一個不為0的元素
				{
					flg=1; //標志找到了不為0的元素
					//先交換行
					for(j=0;j<sqrt(N);j++)
					{
						temp=list[j];
						list[j]=list[i*int(sqrt(N))+j];
						list[i*int(sqrt(N))+j]=temp;
					}
					//再交換列
					for(j=0;j<sqrt(N);j++)
					{
						temp=list[j*int(sqrt(N))];
						list[j*int(sqrt(N))]=list[j*int(sqrt(N))+k];
						list[j*int(sqrt(N))+k]=temp;
					}
					break;
				}
			}
			if(flg==1)
			{
				break;
			}
		}
		if(flg==0)
		{
			//說明該矩陣是0矩陣，其秩為0
			return 0;
		}
	}
	/*cout<<"等價矩陣為:"<<endl;
	for(i=0;i<N;i++) 
	{
		cout<<list[i]<<'\t';
		if(((i+1)%int(sqrt(N)))==0)
		{
			cout<<endl;
		}
	}*/
	for(i=0;i<sqrt(N);i++)
	{
		f=0;
		if(i==0)
		{
			temp=list[i*int(sqrt(N))+i];
			for(j=0;j<sqrt(N);j++) //使對角線上的元素為1
			{
				list[i*int(sqrt(N))+j]/=temp;
				k=i;
			}
		}
		else
		{
			if(i>1)
			{
				temp=list[(i-1)*int(sqrt(N))+(i-1)];
				for(j=0;j<sqrt(N);j++)//使對角線上的元素為1
				{
					list[(i-1)*int(sqrt(N))+j]/=temp;
					k=i-1;
				}
			}
			for(t=i;t<sqrt(N);t++) //消元，把第i-1行乘一個系數依次往下加
			{
				temp=list[t*int(sqrt(N))+k];
				for(j=0;j<sqrt(N);j++)
				{
					list[t*int(sqrt(N))+j]-=list[(i-1)*int(sqrt(N))+j]*temp;
				}
			}
		}
		if(list[i*int(sqrt(N))+i]==0)//確保以后迭代的對角線元素不為0
		{
			for(k=i;k<sqrt(N);k++) //行
			{
				for(t=i;t<sqrt(N);t++)//列
				{
					if(list[k*int(sqrt(N))+t]!=0)//找到一個不為零的元素
					{
						f=1;
						//先換行
						for(j=i;j<sqrt(N);j++)
						{
							temp=list[i*int(sqrt(N))+j];
							list[i*int(sqrt(N))+j]=list[k*int(sqrt(N))+j];
							list[k*int(sqrt(N))+j]=temp;
						}
						//再換列
						for(j=i;j<sqrt(N);j++)
						{
							temp=list[j*int(sqrt(N))+i];
							list[j*int(sqrt(N))+i]=list[j*int(sqrt(N))+t];
							list[j*int(sqrt(N))+t]=temp;
						}
						break;
					}
				}
				if(f==1)//找到了不為零的元素，交換后就跳出
				{
					break;
				}
			}
			if (f==0)
			{
				//說明剩下的元素全為0,直接跳出整個for循環
				break;
			}
		}
	}
	for(i=0;i<sqrt(N);i++)//最后秩代為上三角陣后，每一行只要就一個元素不為0，該矩陣的秩就加1
	{
		for(j=0;j<sqrt(N);j++)
		{
			if(list[i*int(sqrt(N))+j]!=0)
			{
				s++;
				break;
			}

		}

	}
	/*cout<<"迭代為上三角矩陣為:"<<endl;
	for(i=0;i<N;i++) 
	{
		cout<<list[i]<<'\t';
		if(((i+1)%int(sqrt(N)))==0)
		{
			cout<<endl;
		}
	}*/
	return s;
}
float Matrix::hang()//求行列式，該算法還可以
{
	int i,j,k,t;
	float s=1.0;
	float temp=0;
	int flg=0;
	if(list[0]==0)//確保迭代的第一個元素不為0
	{
		for(i=1;i<sqrt(N);i++)
		{
			if(list[i*int(sqrt(N))]!=0)
			{
				flg=1;
				for(j=0;j<sqrt(N);j++)
				{
					temp=list[j];
					list[j]=list[i*int(sqrt(N))+j];
					list[i*int(sqrt(N))+j]=temp;
				}
				break;
			}
		}
		if(flg==0)//說明行列式為0
		{
			return 0;
		}
	}
	for(i=0;i<sqrt(N);i++)
	{
		flg=0;
		if(i==0)
		{
			temp=list[i*int(sqrt(N))+i];
			for(j=0;j<sqrt(N);j++)
			{
				list[i*int(sqrt(N))+j]/=temp;
				k=i;
			}
		}
		else
		{
			if(i>1)
			{
				temp=list[(i-1)*int(sqrt(N))+(i-1)];
				for(j=0;j<sqrt(N);j++)
				{
					list[(i-1)*int(sqrt(N))+j]/=temp;
					k=i-1;
				}
			}
			for(t=i;t<sqrt(N);t++)//消元
			{
				temp=list[t*int(sqrt(N))+k];
				for(j=0;j<sqrt(N);j++)
				{
					list[t*int(sqrt(N))+j]-=list[(i-1)*int(sqrt(N))+j]*temp;
				}
			}
		}
		if(list[i*int(sqrt(N))+i]==0)//確保以后迭代的對角線元素元素不為0
		{
			for(k=i+1;k<sqrt(N);k++)
			{
				if(list[k*int(sqrt(N))+i]!=0)
				{
					flg=1;
					for(j=0;j<sqrt(N);j++)
					{
						temp=list[i*int(sqrt(N))+j];
						list[i*int(sqrt(N))+j]=list[k*int(sqrt(N))+j];
						list[k*int(sqrt(N))+j]=temp;
					}
					break;
				}
			}
			if (flg==0)
			{
				return 0;
			}
		}
		s*=list[i*int(sqrt(N))+i];
	}
	return s;
}
void Matrix::ni()//求逆，經典的高斯算法
{
	int i,j,k,t;
	float ru[2*N]={0};
	float temp=0;
	int flg=0;
	//在該矩陣右邊填充一個同階的單位陣，為以下的求解作好準備
	for (i=0;i<sqrt(N);i++)
	{
		for(j=0;j<2*sqrt(N);j++)
		{
			if(j<sqrt(N))
			{
				ru[i*int(2*sqrt(N))+j]=list[i*int(sqrt(N))+j];
			}
			else
			{
				if(((j-i)-int(sqrt(N)))==0)
				{
					ru[i*int(2*sqrt(N))+j]=1;
				}
				else
				{
					ru[i*int(2*sqrt(N))+j]=0;
				}
			}
		}
	}
	/*cout<<"初始矩陣為："<<endl;
	for(i=0;i<int(sqrt(N));i++)
	{
		for(j=0;j<2*int(sqrt(N));j++)
		{
			cout<<ru[i*int(2*sqrt(N))+j]<<'\t';
		}
		cout<<endl;
	}*/
	if(ru[0]==0)//確保迭代的第一個元素不為0
	{
		for(i=1;i<int(sqrt(N));i++)
		{
			if(ru[i*int(sqrt(N))]!=0)
			{
				flg=1;
				for(j=0;j<2*int(sqrt(N));j++)
				{
					temp=ru[j];
					ru[j]=ru[i*int(2*sqrt(N))+j];
					ru[i*int(2*sqrt(N))+j]=temp;
				}
				break;
			}
		}
		if (flg==0)//說明該矩陣不可逆，直接跳出 
		{
			return;
		}
	}
	//先迭代為一個上三角形
	for(i=0;i<int(sqrt(N));i++)
	{
		flg=0;
		if(i==0)
		{
			temp=ru[i*int(2*sqrt(N))+i];
			for(j=0;j<2*int(sqrt(N));j++)
			{
				ru[i*int(2*sqrt(N))+j]/=temp;
				k=i;
			}
		}
		else
		{
			if(i>1)
			{
				temp=ru[(i-1)*int(2*sqrt(N))+(i-1)];
				for(j=0;j<2*int(sqrt(N));j++)
				{
					ru[(i-1)*int(2*sqrt(N))+j]/=temp;
					k=i-1;
				}
			}
			for(t=i;t<int(sqrt(N));t++) //消元
			{
				temp=ru[t*int(2*sqrt(N))+k];
				for(j=0;j<2*int(sqrt(N));j++)
				{
					ru[t*int(2*sqrt(N))+j]-=ru[(i-1)*int(2*sqrt(N))+j]*temp;
				}
			}
		}
		if(ru[i*int(2*sqrt(N))+i]==0)//確保迭代后對角線上元素不為0
		{
			for(k=i+1;k<int(sqrt(N));k++)
			{
				if(ru[k*int(2*sqrt(N))+i]!=0)
				{
					flg=1;
					for(j=0;j<2*int(sqrt(N));j++)
					{
						temp=ru[i*int(2*sqrt(N))+j];
						ru[i*int(2*sqrt(N))+j]=ru[k*int(2*sqrt(N))+j];
						ru[k*int(2*sqrt(N))+j]=temp;
					}
					break;
				}
			}
			if (flg==0)//說明該矩陣不可逆，直接跳出
			{
				return ;
			}
		}
		if (i==(int(sqrt(N))-1))//把最后一行元素對角線上的元素化為1
		{
			temp=ru[i*int(2*sqrt(N))+i];
			for(j=0;j<2*int(sqrt(N));j++)
			{
				ru[i*int(2*sqrt(N))+j]/=temp;
			}
		}
	}
	/*cout<<"矩陣迭代為上三角形后為："<<endl;
	for(i=0;i<int(sqrt(N));i++)
	{
		for(j=0;j<2*int(sqrt(N));j++)
		{
			cout<<ru[i*int(2*sqrt(N))+j]<<'\t';
		}
		cout<<endl;
	}*/
	//再轉換為對角形
	for(i=int(sqrt(N))-1;i>=0;i--)
	{
		for(t=i-1;t>=0;t--)
		{
			temp=ru[t*int(2*sqrt(N))+i];
			for(j=0;j<2*int(sqrt(N));j++)
			{
				ru[t*int(2*sqrt(N))+j]-=ru[i*int(2*sqrt(N))+j]*temp;
			}
		}

	}
	/*cout<<"矩陣迭代為對角形后為："<<endl;
	for(i=0;i<int(sqrt(N));i++)
	{
		for(j=0;j<2*int(sqrt(N));j++)
		{
			cout<<ru[i*int(2*sqrt(N))+j]<<'\t';
		}
		cout<<endl;
	}*/
	//從填充的矩陣中找出所求的逆陣
	for(i=0;i<int(sqrt(N));i++)
	{
		for(j=int(sqrt(N));j<2*int(sqrt(N));j++)
		{
			list[i*int(sqrt(N))+(j-int(sqrt(N)))]=ru[i*int(2*sqrt(N))+j];
		}
	}
}
void Matrix::dicheng(float alpha)
{
	for(int i=0;i<N;i++)
	{
		list[i]=alpha*list[i];
	}
}
void Matrix:: T() //轉置
{
	int i,j;
	Matrix temp;
	for(i=0;i<sqrt(N);i++)
	{
		for(j=0;j<sqrt(N);j++)
		{
			temp[j*int(sqrt(N))+i]=list[i*int(sqrt(N))+j];
		}
	}
	for(i=0;i<sqrt(N);i++)
	{
		for(j=0;j<sqrt(N);j++)
		{
			list[i*int(sqrt(N))+j]=temp[i*int(sqrt(N))+j];
		}
	}
}
Matrix operator+(Matrix a,Matrix b)
{
	Matrix temp;
	for(int i=0;i<N;i++)
	{
		temp[i]=a[i]+b[i];
	}
	return temp;
}
Matrix operator-(Matrix a,Matrix b)
{
	Matrix temp;
	for(int i=0;i<N;i++)
	{
		temp[i]=a[i]-b[i];
	}
	return temp;
}
Matrix operator*(Matrix a,Matrix b)
{
	int i,j,k;
	Matrix temp;
	for(i=0;i<N;i++) //初始化為零
	{
		temp.Elem(0);
	}
	for(i=0;i<sqrt(N);i++)
	{
		for(j=0;j<sqrt(N);j++)
		{
			for(k=0;k<sqrt(N);k++)
			{
				temp[i*int(sqrt(N))+j]+=a[i*int(sqrt(N))+k]*b[k*int(sqrt(N))+j];
			}
		}
	}
	return temp;
}
float & Matrix::operator[](int index)
{
	if(index>=N||index<0)
	{
		cout<<"\nerror:下標"<<index<<"越界"<<endl;
		exit(2);
	}
	return list[index];
}
Matrix Matrix::operator =(Matrix &a)
{
	for(int i=0;i<N;i++)
	{
		list[i]=a[i];
	}
	return *this;
}
Matrix Matrix::operator *(float alpha)//點乘
{
	for(int i=0;i<N;i++)
	{
		list[i]=alpha*list[i];
	}
	return *this;
}
void main()
{
	Matrix list1,list2,list3;
	double time1=0;	//開始時間                
	double time2=0;	//結束時間	
	float ha=0;
	int z=0;
	float tmp1[9]={1,2,0,1,2,0,1,1,3};
	float tmp2[9]={1,2,3,4};
	time1=(double)clock();
	for(int i=0;i<N;i++) 
	{
		list1.Elem(tmp1[i]);
		list2.Elem(tmp2[i]);
	}
	list1.Print();
	list2.Print();
	z=list1.zhi();
	cout<<"矩陣秩為:"<<z<<endl;
	list3=list1*2;
	cout<<"矩陣點乘后為:"<<endl;
    list3.Print();
	list2.T();
	cout<<"矩陣轉置后為:"<<endl;
	list2.Print();
	//list3=list1*list2;
	cout<<"矩陣相乘后為:"<<endl;
	list3.Print();
	ha=list2.hang();
	cout<<"矩陣行列式為:"<<ha<<endl;
	if(list1.hang()!=0)
	{
		list1.ni();
		cout<<"矩陣逆矩陣為:"<<endl;
		list1.Print();
	}
	else
	{
		cout<<"矩陣不可逆!不存在逆矩陣!"<<endl;
	}
	//cout<<list1.GetElem(list1.Elem()-1)<<endl;
	//list1.Destroy();
	time2=(double)clock(); 
	cout<<"運行所需時間為："<<(time2-time1)/1000<<" 秒"<<endl; 
}