realmatrix.cpp

这是一个实数矩阵计算的源代码

		}
		return newMatrix;
	}
	if(pHead!=NULL)
	{
		newMatrix.Create(pHead->nRow,pHead->nCol+1);
		if(col>pHead->nCol) col=pHead->nCol+1;
		if(col<1) col=1;
		
		//copy data to newMatrix

		CNode* newpData=new CNode;
		CNode* newpCurrent=new CNode;
		newpCurrent=newMatrix.pHead->pNext;
		newpData=newpCurrent;

		src.pData=src.pHead->pNext;
		pData=pHead->pNext;
		pCurrent=pData;
		while(pCurrent!=pTail)
		{
			while(pData->pRight!=NULL&&pData->pRight->nCol<col)
			{

				CNode* newpData1=new CNode;
				pData=pData->pRight;
				newpData1->nRow=pData->nRow;
				newpData1->nCol=pData->nCol;
				newpData1->data=pData->data;
				newpData1->pDown=NULL;
				newpData1->pRight=NULL;
				newpData->pRight=newpData1;
				newpData=newpData->pRight;

				CNode *pTemp=newMatrix.pHead;
				for(int i=0;i<newpData->nCol;++i)
					pTemp=pTemp->pNext;
				while(pTemp->pDown!=NULL)
					pTemp=pTemp->pDown;
				pTemp->pDown=newpData;
			}

			if(src.pData!=src.pTail&&pData->pRight->nCol==col)
			{
				if(src.pData->pRight!=NULL&&src.pData->pRight->nCol==1)
				{
					CNode* pTemp1=new CNode;
					pTemp1->nRow=src.pData->pRight->nRow;
					pTemp1->nCol=pData->pRight->nCol;
					pTemp1->data=src.pData->pRight->data;
					pTemp1->pRight=NULL;
					pTemp1->pDown=NULL;
					newpData->pRight=pTemp1;
					newpData=newpData->pRight;
	
					CNode* pTemp=newMatrix.pHead;
					for(int i=0;i<newpData->nCol;++i)
						pTemp=pTemp->pNext;
					while(pTemp->pDown!=NULL)
						pTemp=pTemp->pDown;
					pTemp->pDown=newpData;
				}
				src.pData=src.pData->pNext;
			}
			while(pData->pRight!=NULL&&pData->pRight->nCol>=col)
			{

				CNode* newpData1=new CNode;
				pData=pData->pRight;
				newpData1->nRow=pData->nRow;
				newpData1->nCol=pData->nCol+1;
				newpData1->data=pData->data;
				newpData1->pDown=NULL;
				newpData1->pRight=NULL;
				newpData->pRight=newpData1;
				newpData=newpData->pRight;

				CNode *pTemp=newMatrix.pHead;
				for(int i=0;i<newpData->nCol;++i)
					pTemp=pTemp->pNext;
				while(pTemp->pDown!=NULL)
					pTemp=pTemp->pDown;
				pTemp->pDown=newpData;
			}
			pCurrent=pCurrent->pNext;
			pData=pCurrent;
			newpCurrent=newpCurrent->pNext;
			newpData=newpCurrent;
		}
	}
	return newMatrix;
}

CRealMatrix CRealMatrix::PoseBlock(int& number,int* Ip,CRealMatrix& MatrixP)
{
	int n=pHead->nRow;
	CRealMatrix newMatrix(n,n);
	if(pHead!=NULL&&pHead->nCol==pHead->nRow)
	{
		int* Iq=new int[n+1];		//存放换行的序号
		int* MaxCol=new int[n+1];
		int* MinCol=new int[n+1];
		if(!MatrixP.IsEmpty())
			MatrixP.Empty();
		MatrixP.Create(n,n);	//存放变换矩阵



		//统计每行最后一个非零元素的列号

		pData=pHead->pNext;
		pCurrent=pData;
		while(pCurrent!=pTail)
		{
			while(pData->pRight!=NULL)
			{
				pData=pData->pRight;
			}
			if(pCurrent->pRight!=NULL)
				MinCol[pCurrent->nRow]=pCurrent->pRight->nCol;
			else
				MinCol[pCurrent->nRow]=0;
			if(pData!=pCurrent)
				MaxCol[pCurrent->nRow]=pData->nCol;
			else
				MaxCol[pCurrent->nRow]=0;
			pCurrent=pCurrent->pNext;
			pData=pCurrent;
		}

		//按最大列号从小到大进行排序

		for(int i=1;i<=n;++i)
		{
			int m_wMin=0;
			for(int j=1;j<=n;++j)
			{
				bool lp1=0;

				for(int m=1;m<i;++m)
				{
					if(j==Iq[m]) 
					{
						lp1=1;
						break;
					}
				}
				if(!lp1) break;

			}
			int p=j;
			m_wMin=MaxCol[p];
			for(j=1;j<=n;++j)
			{
				bool lp=0;
				for(int m=1;m<i;++m)
				{
					if(j==Iq[m]) 
					{
						lp=1;
						break;
					}
				}
				if(lp) continue;

				if(MaxCol[j]<m_wMin)
				{
					m_wMin=MaxCol[j];
					p=j;
				}
			}
			Iq[i]=p;
		}


		//形成排序后的新矩阵
		CNode* newpData=new CNode;
		CNode* newpCurrent=new CNode;

		newpData=newMatrix.pHead->pNext;
		newpCurrent=newpData;

		for(i=1;i<=n;++i)
		{
			pData=pHead;
			for(int j=0;j<Iq[i];++j)
				pData=pData->pNext;

			//copy data
			while(pData->pRight!=NULL)
			{
				CNode* newpData1=new CNode;
				pData=pData->pRight;
				newpData1->nRow=i;
				newpData1->nCol=pData->nCol;
				newpData1->data=pData->data;
				newpData1->pDown=NULL;
				newpData1->pRight=NULL;
				newpData->pRight=newpData1;
				newpData=newpData->pRight;

				//construct col pointer
				CNode *pTemp=newMatrix.pHead;
				for(int k=0;k<newpData->nCol;++k)
					pTemp=pTemp->pNext;
				while(pTemp->pDown!=NULL)
					pTemp=pTemp->pDown;
				pTemp->pDown=newpData;
			}

			//next row
			newpCurrent=newpCurrent->pNext;
			newpData=newpCurrent;
		}
		

		//寻找最小块
		number=1;
		for(i=1;i<=n;++i)
		{
			if(i==MaxCol[Iq[i]]&&MaxCol[Iq[i]]+1==MinCol[Iq[i]+1])
			{
				Ip[number]=i;
				++number;
			}
		}
		Ip[number]=n;
		for(i=1;i<=n;++i)
			MatrixP.Store(i,Iq[i],1);
		delete[] Iq;
		delete[] MinCol;
		delete[] MaxCol;
	}
	return newMatrix;
}

BOOL CRealMatrix::Decompose1()
{
	int n=pHead->nRow;
	if(n==pHead->nCol&&n!=0)
	{
		CNode* psrcCurrent=pHead->pNext;
		CNode* psrcData=psrcCurrent;
		while(psrcCurrent!=pTail)
		{
			pData=psrcCurrent;
			while(pData->pRight!=NULL&&pData->pRight->nCol<=psrcCurrent->nRow)
				pData=pData->pRight;
			psrcData=pData;
			//计算Dii
			double value=GetDiaEle0(pData->nRow,pData->nRow);
			if(value==0)      return	1;//对角元为0

			//计算Uij
			for(INT j=psrcCurrent->nRow+1;j<=n;++j)
				double valueU=GetUpEle0(psrcCurrent->nRow,j);
			CNode* pNode=psrcData;
			if(pNode->nCol>psrcCurrent->nRow)
				pNode->data/=value;
			while(pNode->pRight!=NULL)
			{
				pNode=pNode->pRight;
				pNode->data/=value;
			}
			psrcCurrent=psrcCurrent->pNext;
			psrcData=psrcCurrent;
		}
		return 2;
	}
	return 0;
}

double CRealMatrix::GetUpEle0(INT row,INT col)
{
	double sum=0;
	double valueD,valueL,valueU;
	CNode* pRow=pCurrent;
	CNode* pCol=pHead;
	CNode* pTemp=pHead;

	for(INT k=1;k<=row-1;++k)
	{
		pTemp=pTemp->pNext;
		while(pTemp->nRow<k)
			pTemp=pTemp->pNext;
		pCol=pTemp;
		valueU=GetUpEle(col,pCol);
		if(valueU==0)	continue;
		pCol=pTemp;
		valueL=GetUpEle(row,pCol);
		if(valueL==0)	continue;
		valueD=GetDiaEle(k,pTemp);
		sum+=valueL*valueD*valueU;
	}
	CNode* pNode=pHead->pNext;
	while(pNode->nRow<row)
		pNode=pNode->pNext;
	while(pNode->pRight!=NULL&&pNode->pRight->nCol<=col)
		pNode=pNode->pRight;
	if(pNode->nCol==col)
	{
		pNode->data-=sum;
		if(pNode->data==0)
		{
			DeleteNode(row,col);
			return 0;
		}
		return pNode->data;
	}
	Add(row,col,-sum);
	return -sum;
}
double CRealMatrix::GetDiaEle0(INT row,INT col)
{
	double sum=0;
	double valueD,valueU;
	CNode* pCol=pHead;
	CNode* pTemp=pHead;

	for(INT k=1;k<=row-1;++k)
	{
		pTemp=pTemp->pNext;
		while(pTemp->nRow<k)
			pTemp=pTemp->pNext;
		pCol=pTemp;
		valueU=GetUpEle(row,pCol);
		if(valueU==0) continue;
		valueD=GetDiaEle(k,pTemp);
		sum+=valueU*valueD*valueU;
	}//#end for
	if(pData->nCol==row)
	{
		pData->data-=sum;
		return pData->data;
	}
	Add(row,row,-sum);
	return -sum;
}
void CRealMatrix::Normalize0()
{
	if(pHead==NULL)	return;
	if(pHead->nRow==pHead->nCol&&pHead->nRow>=2)
	{
		int n=pHead->nRow;

		double* dbTemp=new double[n+1];

		//set 0
		for(int i=1;i<=n;++i)
			dbTemp[i]=0;

		CNode* pTemp=new CNode;
		pData=pHead->pNext;
		pCurrent=pData;
		pTemp=pCurrent;
		while(pCurrent!=pTail)
		{
			while(pData->pRight!=NULL)
			{
				pData=pData->pRight;
				if(pData->nCol==pData->nRow)
				{
					dbTemp[pData->nCol]=0;
					continue;
				}
				dbTemp[pData->nCol]=pData->data;
			}

			for(i=1;i<=n;++i)
			{
				if(dbTemp[i]!=0)
				{
					pTemp=pHead->pNext;
					while(pTemp->pNext!=pTail&&pTemp->pNext->nRow<=i)
						pTemp=pTemp->pNext;

					//add a node
					while(pTemp->pRight!=NULL&&pTemp->pRight->nCol<=pCurrent->nRow)
						pTemp=pTemp->pRight;

					if(pTemp->nCol==pCurrent->nRow&&pTemp->nRow==i)	continue;
					
					CNode* pNode=new CNode;
					pNode->pDown=NULL;
					pNode->pNext=NULL;
					pNode->pRight=NULL;
					pNode->pPrev=NULL;

					pNode->nRow=i;
					pNode->nCol=pCurrent->nRow;
					pNode->data=dbTemp[i];

					CNode* pNew=new CNode;
					if(pTemp->pRight!=NULL)
					{
						pNew=pTemp->pRight;
						pNode->pRight=pNew;
					}
					pTemp->pRight=pNode;


					//col table

					CNode* pDTemp=new CNode;
					pDTemp=pCurrent;
					while(pDTemp->pDown!=NULL)
						pDTemp=pDTemp->pDown;
					pDTemp->pDown=pNode;
				}
			}
			
			pCurrent=pCurrent->pNext;
			pData=pCurrent;
			for(i=1;i<=n;++i)
				dbTemp[i]=0;
		}//# end while
		delete[] dbTemp;
	}//# edn if
}//# end


void CRealMatrix::Normalize1()
{
	if(pHead->nRow==0||(pHead->nRow!=pHead->nCol))		return;

	pCurrent=pHead->pNext;
	pData=pCurrent;

	while(pCurrent!=pTail)
	{
		if(pData->pRight==NULL) continue;
		CNode* pTemp=pData;
		pData=pData->pRight;
		while(pData->pRight!=NULL&&pData->pRight->nCol<=pCurrent->nRow)
			DeleteNode(pCurrent->nRow,pData->nCol,pTemp);
		pCurrent=pCurrent->pNext;
		pData=pCurrent;
	}
}

void CRealMatrix::Normalize2()
{
	if(pHead->nRow==0||(pHead->nRow!=pHead->nCol))		return;

	pCurrent=pHead->pNext;
	pData=pCurrent;

	while(pCurrent!=pTail)
	{
		if(pData->pRight==NULL)
		{
			pCurrent=pCurrent->pNext;
			pData=pCurrent;
			continue;
		}
		while(pData->pRight!=NULL&&pData->pRight->nCol<=pCurrent->nRow)
			pData=pData->pRight;

		CNode* pTemp=pData;

		if(pData->nCol<pCurrent->nRow||pData->pRight==NULL)
		{
			pCurrent=pCurrent->pNext;
			pData=pCurrent;
			continue;
		}
		pData=pData->pRight;

		while(pData!=NULL&&pData->nCol>pCurrent->nRow)
			DeleteNode(pCurrent->nRow,pData->nCol,pTemp);
		pCurrent=pCurrent->pNext;
		pData=pCurrent;
	}
}
 

/////////////////////////////////////////////////////////////////
//other operation
///////////////////////////////////////////////////////

CRealMatrix CRealMatrix::CreateG0(CRealMatrix &srcR,CRealMatrix &srcHR)
{

	CRealMatrix srcHRH;

	if(pHead==NULL||srcR.pHead==NULL) return srcHRH;

	if(pHead->nRow==srcR.pHead->nRow)
	{
		double* dbMatrixR=new double[pHead->nRow+1];
		int nRow=pHead->nRow;
		int nCol=pHead->nCol;

		//given n*m dimension of srcHR;
		if(!srcHR.IsEmpty())
			srcHR.Empty();

		//give n*n dimension of srcHRH
		srcHRH.Create(nCol,nCol);

		//take the value of srcR into dbmatrix 

		srcR.pData=srcR.pHead->pNext;
		srcR.pCurrent=srcR.pData;
		while(srcR.pCurrent!=srcR.pTail)
		{
			if(srcR.pData->pRight!=NULL)
			{
				srcR.pData=srcR.pData->pRight;
				dbMatrixR[srcR.pData->nRow]=srcR.pData->data;
			}
			else
				dbMatrixR[srcR.pData->nRow]=srcR.pData->data;
			srcR.pCurrent=srcR.pCurrent->pNext;
			srcR.pData=srcR.pCurrent;
		}

		//concalute srcHR----H(T)R and form the triangle of H(T)RH
		srcHR=*this;

		srcHR.pData=srcHR.pHead->pNext;
		srcHR.pCurrent=srcHR.pData;
		while(srcHR.pCurrent!=srcHR.pTail)
		{
			while(srcHR.pData->pRight!=NULL)
			{
				srcHR.pData=srcHR.pData->pRight;
				srcHR.pData->data*=dbMatrixR[srcHR.pData->nRow];
			}
			srcHR.pCurrent=srcHR.pCurrent->pNext;
			srcHR.pData=srcHR.pCurrent;
		}
		srcHR=srcHR.Transpose();

		//caluation HRH
		srcHR.pCurrent=srcHR.pHead->pNext;
		srcHR.pData=srcHR.pCurrent;
		CNode* pTemp;
		while(srcHR.pCurrent!=srcHR.pTail&&srcHR.pCurrent->nRow!=0)
		{
			double sum=0;
			while(srcHR.pData->pRight!=NULL)
			{
				srcHR.pData=srcHR.pData->pRight;
				if(dbMatrixR[srcHR.pData->nCol]==0) continue;
				sum+=srcHR.pData->data*
					srcHR.pData->data/dbMatrixR[srcHR.pData->nCol];
			}
			srcHRH.Store(srcHR.pCurrent->nRow,srcHR.pCurrent->nRow,sum);

			pTemp=srcHR.pCurrent->pNext;
			CNode* pTemp1;
			while(pTemp->nRow!=0&&pTemp!=srcHR.pTail)
			{
				sum=0;
				srcHR.pData=srcHR.pCurrent;
				pTemp1=pTemp;
				while(srcHR.pData->pRight!=NULL)
				{
					srcHR.pData=srcHR.pData->pRight;
					while(pTemp->pRight!=NULL&&
						pTemp->pRight->nCol<srcHR.pData->nCol)
						pTemp=pTemp->pRight;
					if(pTemp->pRight!=NULL&&pTemp->pRight->nCol==srcHR.pData->nCol)
					{
						pTemp=pTemp->pRight;
						sum+=srcHR.pData->data*pTemp->data/dbMatrixR[srcHR.pData->nCol];
					}
					else if(pTemp->pRight==NULL)
						break;
				}
				if(sum!=0)
					srcHRH.Store(srcHR.pCurrent->nRow,pTemp1->nRow,sum);
				pTemp=pTemp1->pNext;
			}
			srcHR.pCurrent=srcHR.pCurrent->pNext;
			srcHR.pData=srcHR.pCurrent;
		}
		delete[] dbMatrixR;
	} //#end if
	return srcHRH;
}

CRealMatrix CRealMatrix::CreateG1(CRealMatrix &srcR,CRealMatrix &srcHR)
{

	CRealMatrix srcHRH;

	if(pHead==NULL||srcR.pHead==NULL) return srcHRH;

	if(pHead->nRow==srcR.pHead->nRow)
	{
		double* dbMatrixR=new double[pHead->nRow+1];
		int nRow=pHead->nRow;
		int nCol=pHead->nCol;

		//given n*m dimension of srcHR;
		if(!srcHR.IsEmpty())
			srcHR.Empty();

		//give n*n dimension of srcHRH
		srcHRH.Create(nCol,nCol);

		//take the value of srcR into dbmatrix 

		srcR.pData=srcR.pHead->pNext;
		srcR.pCurrent=srcR.pData;
		while(srcR.pCurrent!=srcR.pTail)
		{
			if(srcR.pData->pRight!=NULL)
			{