// Matrix.cpp: implementation of the CMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ProGram2.h"
#include "Matrix.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

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

CMatrix::CMatrix()
{
	m_NumCols=1;
	m_NumRows=1;
	m_pData=NULL;
	Init(m_NumRows,m_NumCols);
}

CMatrix::CMatrix(int nRows, int nCols)
{////////////指定行列構造函數
	m_NumCols=nCols;
	m_NumRows=nRows;
	m_pData=NULL;
	Init(m_NumRows,m_NumCols);
}
CMatrix::CMatrix(int nRows,int nCols, double value[])
{//、、、、指定值構造函數
	m_NumCols=nCols;
	m_NumRows=nRows;
	m_pData=NULL;
	Init(m_NumRows,m_NumCols);
	SetData(value);
	
}
CMatrix::CMatrix(int nSize)
{//方陣構造函數
	m_NumCols=nSize;
	m_NumRows=nSize;
	m_pData=NULL;
	Init(nSize,nSize);
}
CMatrix::CMatrix(int nSize,double value[])
{//指定值的方陣構造函數
	m_NumCols=nSize;
	m_NumRows=nSize;
	m_pData=NULL;
	Init(nSize,nSize);
	SetData(value);
}
CMatrix::CMatrix(const CMatrix &other)
{
	m_NumCols=other.m_NumCols;
	m_NumRows=other.m_NumRows;
	m_pData=NULL;
	Init(m_NumRows,m_NumCols);
	memcpy(m_pData,other.m_pData,sizeof(double)*m_NumRows*m_NumCols);
}
////
////
////初始化//生成0矩陣
bool CMatrix::Init(int nRows,int nCols)
{
	if(m_pData)
	{
		delete[] m_pData;
		m_pData=NULL;
	}
	m_NumRows=nRows;
	m_NumCols=nCols;
	int nSize=nRows*nCols;
	if(nSize<0)
		return false;

	//分配內存
	m_pData=new double[nSize];
	if(m_pData==NULL)
		return false;
	if(IsBadReadPtr(m_pData,sizeof(double)*nSize))////////////
		return false;
	memset(m_pData,0,sizeof(double)*nSize);
	return true;
}
bool CMatrix::MakeUnitMatrix(int nSize)
{
	if(!Init(nSize,nSize))//初始化為0矩陣
		return false;
	for(int i=0;i<nSize;i++)
		SetElement(i,i,1);
	return true;
}
///
///
///析構函數
CMatrix::~CMatrix()
{
	if(m_pData)
	{
		delete[] m_pData;
		m_pData=NULL;
	}
}
//////////////////////////////////////////////////////////
///
///
///元素與值操作////////////
///
///
///設置矩陣各個元素的值
void CMatrix::SetData(double value[])
{
	memset(m_pData,0,sizeof(double)*m_NumRows*m_NumCols);//設置空值
	memcpy(m_pData,value,sizeof(double)*m_NumRows*m_NumCols);
}
///
///
///設置指定元素的值
bool CMatrix::SetElement(int nRow,int nCol,double value)
{
	if(nCol<0 ||nRow<0 ||nCol>m_NumCols || nRow>m_NumRows)
		return false;
	if(m_pData==NULL)
		return false;
	m_pData[nCol+nRow*m_NumCols]=value;
	return true;
}
///
///
///獲取指定元素的值
double CMatrix::GetElement(int nRow,int nCol)const
{
	ASSERT(nCol>=0 && nRow>=0 && nCol<m_NumCols && nRow<m_NumRows);
	ASSERT(m_pData);
	return m_pData[nCol+nRow*m_NumCols];
}
///
///
///獲取行列數
int CMatrix::GetNumCols()const
{
	return m_NumCols;
}
int CMatrix::GetNumRows()const
{
	return m_NumRows;
}
///
///
///獲取矩陣的數據
double *CMatrix::GetData()const
{
	return m_pData;
}    
///
///
///獲取指定行的向量
int CMatrix::GetRowVector(int nRow,double *pVector)const
{
	if(m_pData==NULL)
		delete pVector;
	pVector = new double[m_NumCols];
	ASSERT(pVector!=NULL);
	for(int j=0;j<m_NumCols;j++)
		pVector[j]=GetElement(nRow,j);
	return m_NumCols;
}                                          
///
///
///獲取指定列的向量
int CMatrix::GetColVector(int nCol,double *pVector)const
{
	if(m_pData==NULL)
		delete pVector;
	pVector = new double[m_NumRows];
	ASSERT(pVector!=NULL);
	for(int j=0;j<m_NumRows;j++)
		pVector[j]=GetElement(j,nCol);
	return m_NumCols;	
}
////////////////////////////////////////////
///
///
///數學運算符重載
///
///
///'='號，賦值
CMatrix &CMatrix::operator=(const CMatrix &other)
{
	if(&other!=this)
	{
		m_NumRows=other.GetNumRows();
		m_NumCols=other.GetNumCols();
		Init(m_NumRows,m_NumCols);
		memcpy(m_pData,other.GetData(),sizeof(double)*m_NumRows*m_NumCols);////////改動
	}
	return *this;

}
///
///
///'=='號，判斷相等
bool CMatrix::operator==(const CMatrix &other)const
{
	if(m_NumRows!=other.GetNumRows() || m_NumCols!=other.GetNumCols())
		return false;
	for(int i=0;i<m_NumRows;i++)
	{
		for(int j=0;j<m_NumCols;j++)
		{
			if(GetElement(i,j)!=other.GetElement(i,j))
				return false;
		}
	}
	return true;
}
///
///
///'!='號，判斷不等
bool CMatrix::operator!=(const CMatrix &other)const
{
	return !(*this==other);
}
///
///
///'+'號，矩陣相加
CMatrix CMatrix::operator+(const CMatrix &other)const
{
	ASSERT(m_NumRows==other.GetNumRows() || m_NumCols==other.GetNumCols());
	CMatrix result(*this);
	for(int i=0;i<m_NumRows;i++)
	{
		for(int j=0;j<m_NumCols;j++)
		{
			result.SetElement(i,j,result.GetElement(i,j)+other.GetElement(i,j));
		}
	}
	return result;
}
///
///
///'-'號，矩陣相減
CMatrix CMatrix::operator-(const CMatrix &other)const
{
	ASSERT(m_NumRows==other.GetNumRows() || m_NumCols==other.GetNumCols());
	CMatrix result(*this);

	for(int i=0;i<m_NumRows;i++)
	{
		for(int j=0;j<m_NumCols;j++)
		{
			result.SetElement(i,j,result.GetElement(i,j)-other.GetElement(i,j));
		}
	}
	return result;
}
///
///
///'*'號，矩陣乘以double數值
CMatrix CMatrix::operator*(double value)const
{
	ASSERT(m_NumRows==this->GetNumRows() || m_NumCols==this->GetNumCols());
	CMatrix result(*this);

	for(int i=0;i<m_NumRows;i++)
	{
		for(int j=0;j<m_NumCols;j++)
		{
			result.SetElement(i,j,result.GetElement(i,j)*value);
		}
	}
	return result;
}
///
///
///'*'號，矩陣相乘
CMatrix CMatrix::operator*(const CMatrix &other)const
{
	ASSERT(m_NumCols==other.GetNumRows());
	CMatrix result(m_NumRows,other.GetNumCols());
	double value;
	for(int i=0;i<result.GetNumRows();i++)
	{
		for(int j=0;j<result.GetNumCols();j++)
		{
			value=0.0;
			for(int k=0;k<m_NumCols;k++)
			{
				value+=GetElement(i,k)*other.GetElement(k,j);
			}
			result.SetElement(i,j,value);
		}
	}
	return result;
}

///
///
///公共函數，求矩陣的逆，返回逆矩陣
CMatrix CMatrix::Invert()
{
	CMatrix result(*this);
	if(result.InvertGaussJordan())
		return result;
}
///
///
///矩陣轉置
CMatrix CMatrix::Transpose()const
{
	CMatrix trans(m_NumCols,m_NumRows);
	for(int i=0;i<m_NumRows;i++)
	{
		for(int j=0;j<m_NumCols;j++)
			trans.SetElement(j,i,GetElement(i,j));
	}
	return trans;
}
///
///
///私有函數，高斯列主元求矩陣的逆
bool CMatrix::InvertGaussJordan()
{
	int *pnRow,*pnCol,i,j,k,l,u,v;
	double d=0.0, p=0.0;

	//分配內存
	pnRow=new int[m_NumCols];
	pnCol=new int[m_NumCols];
	if(pnRow==NULL || pnCol==NULL)
		return false;

	//消元
	for(k=0;k<m_NumCols;k++)
	{
		d=0.0;
		for(i=k;i<m_NumCols;i++)
		{
			for(j=k;j<m_NumCols;j++)
			{
				l=i*m_NumCols+j;
				p=fabs(m_pData[l]);
				if(p>d)
				{
					d=p;
					pnRow[k]=i;
					pnCol[k]=j;
				}
			}
		}

		//失敗
		if(d==0.0)
		{
			delete[] pnRow;
			delete[] pnCol;
			return false;
		}
		if(pnRow[k]!=k)
		{
			for(j=0;j<m_NumCols;j++)
			{
				u=k*m_NumCols+j;
				v=pnRow[k]*m_NumCols+j;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}
		if(pnCol[k]!=k)
		{
			for(i=0;i<m_NumCols;i++)
			{
				u=i*m_NumCols+k;
				v=i*m_NumCols+pnCol[k];
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}			
		}
		l=k*m_NumCols+k;
		m_pData[l]=1.0/m_pData[l];
		for(j=0;j<m_NumCols;j++)
		{
			if(j!=k)
			{
				u=k*m_NumCols+j;
				m_pData[u]=m_pData[u]*m_pData[l];
			}
		}
		for(i=0;i<m_NumCols;i++)
		{
			if(i!=k)
			{
				for(j=0;j<m_NumCols;j++)
				{
					if(j!=k)
					{
						u=i*m_NumCols+j;
						m_pData[u]=m_pData[u]-m_pData[i*m_NumCols+k]*m_pData[k*m_NumCols+j];
					}
				}
			}
		}
		for(i=0;i<m_NumCols;i++)
		{
			if(i!=k)
			{
				u=i*m_NumCols+k;
				m_pData[u]=-m_pData[u]*m_pData[l];
			}
		}
	}
	//調整恢復行列次序
	for(k=m_NumCols-1;k>=0;k--)
	{
		if(pnCol[k]!=k)
		{
			for(j=0;j<m_NumCols;j++)
			{
				u=k*m_NumCols+j;
				v=pnCol[k]*m_NumCols+j;
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}
		}
		if(pnRow[k]!=k)
		{
			for(i=0;i<m_NumCols;i++)
			{
				u=i*m_NumCols+k;
				v=i*m_NumCols+pnRow[k];
				p=m_pData[u];
				m_pData[u]=m_pData[v];
				m_pData[v]=p;
			}			
		}
	}
	//清理內存
	delete[] pnCol;
	delete[] pnRow;
	//成功返回
	return true;
}
///
///
///私有函數，求對稱矩陣的特征值與特征根
bool CMatrix::JacobiEigenv(double EigenValue[],CMatrix &EigenVector,int nMaxIt,double eps)
{

	int i,j,p,q,u,w,t,s,l;
	double fm, cn, sn,omega,x,y,d;
	if(!EigenVector.Init(m_NumCols,m_NumCols))
		return false;
	l=1;
	for(i=0;i<m_NumCols;i++)
	{
		EigenVector.m_pData[i*m_NumCols+i]=1.0;
		for(j=0;j<m_NumCols;j++)
			if(i!=j)
				EigenVector.m_pData[i*m_NumCols+j]=0.0;
	}
	while(1)
	{
		fm=0.0;
		for(i=1;i<m_NumCols;i++)
		{
			for(j=0;j<i;j++)
			{
				d=fabs(m_pData[i*m_NumCols+j]);
				if((i!=j) && (d>fm))
				{
					fm=d;
					p=i;
					q=j;
				}
			}
		}
		if(fm<eps)
		{
			for(i=0;i<m_NumCols;++i)
				EigenValue[i]=GetElement(i,i);
			return true;
		}
		if(l>nMaxIt)
			return false;
		l=l+1;
		u=p*m_NumCols+q;
		w=p*m_NumCols+p;
		t=q*m_NumCols+p;
		s=q*m_NumCols+q;
		x=-m_pData[u];
		y=(m_pData[s]-m_pData[w])/2.0;
		omega=x/sqrt(x*x+y*y);

		if(y<0.0)
			omega=-omega;
		sn=1.0+sqrt(1.0-omega*omega);
		sn=omega/sqrt(2.0*sn);
		cn=sqrt(1.0-sn*sn);
		fm=m_pData[w];
		m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData[u]*omega;
		m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega;
		m_pData[u]=0.0;
		m_pData[t]=0.0;
		for(j=0;j<m_NumCols;j++)
		{
			if((j!=p) && (j!=q))
			{
				u=p*m_NumCols+j;
				w=q*m_NumCols+j;
				fm=m_pData[u];
				m_pData[u]=fm*cn+m_pData[w]*sn;
				m_pData[w]=-fm*sn+m_pData[w]*cn;
			}
		}
		for(i=0;i<m_NumCols;i++)
		{
			if((i!=p) && (i!=q))
			{
				u=i*m_NumCols+p;
				w=i*m_NumCols+q;
				fm=m_pData[u];
				m_pData[u]=fm*cn+m_pData[w]*sn;
				m_pData[w]=-fm*sn+m_pData[w]*cn;
			}
		}
		for(i=0;i<m_NumCols;i++)
		{
			u=i*m_NumCols+p;
			w=i*m_NumCols+q;
			fm=EigenVector.m_pData[u];
			EigenVector.m_pData[u]=fm*cn+EigenVector.m_pData[w]*sn;
			EigenVector.m_pData[w]=-fm*sn+EigenVector.m_pData[w]*cn;
		}
	}
	for(i=0;i<m_NumCols;++i)
		EigenValue[i]=GetElement(i,i);
	return true;
}
///
///
///私有函數，求矩陣特征值和特征向量的第二種方法
bool CMatrix::JacobiEigenv2(double EigenValue[],CMatrix &EigenVector,int nMaxIt,double eps)
{
	int i,j,p,q,u,w,t,s;
	double ff,fm,cn,sn,omega,x,y,d;
	if(!EigenVector.Init(m_NumCols,m_NumCols))
		return false;
	for(i=0;i<m_NumCols;i++)
	{
		EigenVector.m_pData[i*m_NumCols+i]=1.0;
		for(j=0;j<m_NumCols;j++)
			if(i!=j)
				EigenVector.m_pData[i*m_NumCols+j]=0.0;
	}
	ff=0.0;
	for(i=1;i<m_NumCols;i++)
	{
		for(j=0;j<i-1;j++)
		{
			d=m_pData[i*m_NumCols+j];
			ff=ff+d*d;
		}
	}
	ff=sqrt(2.0*ff);
Loop_0:
	ff=ff/(1.0*m_NumCols);
Loop_1:
	for(i=1;i<m_NumCols;i++)
	{
		for(j=0;j<i;j++)
		{
			d=fabs(m_pData[i*m_NumCols+j]);
			if(d>ff)
			{
				p=i;
				q=j;
				goto Loop_2;
			}
		}
	}
	if(ff<eps)
	{
		for(i=0;i<m_NumCols;i++)
			EigenValue[i]=GetElement(i,i);
		return true;
	}
	goto Loop_0;
Loop_2:
	u=p*m_NumCols+q;
	w=p*m_NumCols+p;
	t=q*m_NumCols+p;
	s=q*m_NumCols+q;
	x=-m_pData[u];
	y=(m_pData[s]-m_pData[w])/2.0;
	omega=x/sqrt(x*x+y*y);
	if(y<0.0)
		omega=-omega;
	sn=1.0+sqrt(1.0-omega*omega);
	sn=omega/sqrt(2.0*sn);
	cn=sqrt(1.0-sn*sn);
	fm=m_pData[w];
	m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData[u]*omega;
	m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega;
	m_pData[u]=0.0;
	m_pData[t]=0.0;
	for(j=0;j<m_NumCols;j++)
	{
		if((j!=p)&&(j!=q))
		{
			u=p*m_NumCols+j;
			w=q*m_NumCols+j;
			fm=m_pData[u];
			m_pData[u]=fm*cn+m_pData[w]*sn;
			m_pData[w]=-fm*sn+m_pData[w]*cn;
		}
	}
	for(i=0;i<m_NumCols;i++)
	{
		if((i!=p)&&(i!=q))
		{
			u=i*m_NumCols+p;
			w=i*m_NumCols+q;
			fm=m_pData[u];
			m_pData[u]=fm*cn+m_pData[w]*sn;
			m_pData[w]=-fm*sn+m_pData[w]*cn;
		}
	}
	for(i=0;i<m_NumCols;i++)
	{
		u=i*m_NumCols+p;
		w=i*m_NumCols+q;
		fm=EigenVector.m_pData[u];
		EigenVector.m_pData[u]=fm*cn+EigenVector.m_pData[w]*sn;
		EigenVector.m_pData[w]=-fm*sn+EigenVector.m_pData[w]*cn;
	}
	goto Loop_1;
}
CMatrix CMatrix::Eigenv()
{
	ASSERT(m_NumCols==m_NumRows);
	CMatrix result(*this);
	double *EigenValue=new double[m_NumCols];
	CMatrix EigenVector;
	if(result.JacobiEigenv(EigenValue,EigenVector))
	{
		result.Init(m_NumRows+1,m_NumCols);
		for(int i=0;i<m_NumCols;i++)
		{
			result.SetElement(0,i,EigenValue[i]);
			for(int j=1;j<=m_NumRows;j++)
			{
				result.SetElement(j,i,EigenVector.GetElement(j-1,i));
			}
		}
		 result.Swape();
		 return result;
	}

}
CMatrix CMatrix::Eigenv2()
{
	ASSERT(m_NumCols==m_NumRows);
	CMatrix result(*this);
	double *EigenValue=new double[m_NumCols];
	CMatrix EigenVector;
	if(result.JacobiEigenv2(EigenValue,EigenVector))
	{
		result.Init(m_NumRows+1,m_NumCols);
		for(int i=0;i<m_NumCols;i++)
		{
			result.SetElement(0,i,EigenValue[i]);
			for(int j=1;j<=m_NumRows;j++)
			{
				result.SetElement(j,i,EigenVector.GetElement(j-1,i));
			}
		}
		 result.Swape();
		 return result;
	}

}
void CMatrix::Swape(void)
{
	ASSERT(*this!=NULL);
	double ch;
	for(int i=0;i<m_NumCols-1;i++)
	{
		for(int j=i+1;j<m_NumCols;j++)
		{
			if(GetElement(0,i)<GetElement(0,j))
			{
				for(int k=0;k<m_NumRows;k++)
				{
					ch=GetElement(k,i);
					SetElement(k,i,GetElement(k,j));
					SetElement(k,j,ch);
				}
			}
		}
	}
}