sparsematrix.cpp

针对有限元分析的稀疏矩阵采用一维存储时的矩阵运算编写的

// SparseMatrix.cpp: implementation of the CSparseMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Matrix.h"
#include "SparseMatrix.h"

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

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

CSparseMatrix::CSparseMatrix()
{
	m_bDecomposed=false;
	m_bElementBufEmpty=true;
	m_nRow=1;
	m_aiDiag=new unsigned long[1];
	m_aiDiag[0]=0;
	m_adValue=new double[1]; 
}

CSparseMatrix::CSparseMatrix(int nRow,unsigned long* aiDiag)
{
	Realloc(nRow,aiDiag);
}

CSparseMatrix::~CSparseMatrix()
{
	delete m_aiDiag;
	delete m_adValue;
}

CSparseMatrix& CSparseMatrix::operator =(const double dBuf)
{
	if(m_bElementBufEmpty){
		ElementBufRealloc();
	}
	unsigned long nElement=m_aiDiag[m_nRow-1]+1;
	for(unsigned long loop=0;loop<nElement;loop++)
		m_adValue[loop]=dBuf;
	return *this;
}

CSparseMatrix& CSparseMatrix::operator*=(const double dBuf)
{
	if(m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while multiplying";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		return *this;
	}
	unsigned long nElement=m_aiDiag[m_nRow-1]+1;
	for(unsigned long loop=0;loop<nElement;loop++)
		m_adValue[loop]*=dBuf;
	return *this;
}

CSparseMatrix& CSparseMatrix::operator =(CSparseMatrix& sm)
{
	unsigned long  nElement,nElement1;
	nElement=m_aiDiag[m_nRow-1]+1;
	nElement1=sm.m_aiDiag[sm.m_nRow-1]+1;
	if(m_nRow != sm.m_nRow || nElement != nElement1){
		Realloc(sm.m_nRow,sm.m_aiDiag);
	}
	memcpy(m_aiDiag,sm.m_aiDiag,m_nRow*sizeof(long));
	m_bElementBufEmpty=sm.m_bElementBufEmpty;
	if(!sm.m_bElementBufEmpty){
		ElementBufRealloc();
		memcpy(m_adValue,sm.m_adValue,nElement1*sizeof(double));
	}
	m_bDecomposed=sm.m_bDecomposed;
	return *this;
}

CSparseMatrix& CSparseMatrix::operator+=(CSparseMatrix& sm)
{
	bool bBandWidthEqual;
	unsigned long loop,loop1,nElement,nBandWidth,nRow;
	if(m_bElementBufEmpty||sm.m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while adding";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		return *this;
	}
	nRow=m_nRow;
	if(nRow>sm.m_nRow) nRow=sm.m_nRow;
	bBandWidthEqual=true;
	for(loop=0;loop<nRow;loop++){
		if(m_aiDiag[loop]!=sm.m_aiDiag[loop]){
			bBandWidthEqual=false;
			break;
		}
	}

	if(bBandWidthEqual){
		nElement=m_aiDiag[nRow-1]+1;
		for(loop=0;loop<nElement;loop++)
			m_adValue[loop]+=sm.m_adValue[loop];
	}
	else{
		m_adValue[0]+=sm.m_adValue[0];
		for(loop=1;loop<nRow;loop++){
			nBandWidth=m_aiDiag[loop]-m_aiDiag[loop-1];
			if(nBandWidth>sm.m_aiDiag[loop]-sm.m_aiDiag[loop-1])
				nBandWidth=sm.m_aiDiag[loop]-sm.m_aiDiag[loop-1];
			for(loop1=0;loop1<nBandWidth;loop1++)
				m_adValue[m_aiDiag[loop]-loop1]+=
					sm.m_adValue[sm.m_aiDiag[loop]-loop1];
		}
	}
	m_bDecomposed=false;
	return *this;
}

void CSparseMatrix::Realloc(int nRow, unsigned long* aiDiag)
{
	ElementBufEmpty();
	delete m_aiDiag;
	m_nRow=nRow;
	m_aiDiag=new unsigned long [m_nRow];
	for(int loop=0;loop<nRow;loop++)
		m_aiDiag[loop]=aiDiag[loop];
}

void CSparseMatrix::ElementBufEmpty()
{
	m_bDecomposed=false;
	m_bElementBufEmpty=true;
	delete m_adValue;
	m_adValue=new double[1]; 
}

void CSparseMatrix::ElementBufRealloc()
{
	unsigned long  nElement=m_aiDiag[m_nRow-1]+1;
	delete m_adValue;
	m_adValue=new double [nElement];
	m_bElementBufEmpty=false;
}

CMatrix CSparseMatrix::operator*(CMatrix& m)
{
	int loop,loop1,loop2,iBuf;
	int nRow,nMatCol,nMatRow,nBandWidth;
	double dBuf;

	if(m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while multiplying";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		CMatrix mTmp;
		return mTmp;
	}

	nRow=nMatRow=m.GetRow();
	nMatCol=m.GetCol();

	CMatrix mTmp(nRow,nMatCol);
	mTmp=0.0;
	for(loop=0;loop<nMatCol;loop++){
		for(loop1=0;loop1<nRow;loop1++){
			if(loop1==0)
				nBandWidth=1;
			else
				nBandWidth=m_aiDiag[loop1]-m_aiDiag[loop1-1];
			iBuf=m_aiDiag[loop1];
			for(loop2=0;loop2<nBandWidth-1;loop2++){
				dBuf=m_adValue[iBuf-nBandWidth+loop2+1];
				mTmp(loop1,loop)+=dBuf*m(loop1-nBandWidth+loop2+1,loop);
				mTmp(loop1-nBandWidth+loop2+1,loop)+=dBuf*m(loop1,loop);
			}
			mTmp(loop1,loop)+=m_adValue[iBuf]*m(loop1,loop);
		}
	}
	return mTmp;
}

bool CSparseMatrix::LdltSolve(int nRow,double *adB)
{
	long loop;
	unsigned long loop1,loop2;
	unsigned long  iBegin,jBegin,iAdd_ii,iAdd_jj,iAdd_kk;
	unsigned long  iRow,nBandWidth;
	double 	dBuf;

	if(m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while ldlt solving";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		return false;
	}

	if(m_bDecomposed==false){
		for(loop=1;loop<nRow;loop++)
		{
			iBegin=loop-m_aiDiag[loop]+m_aiDiag[loop-1]+1;
			iAdd_ii=m_aiDiag[loop];
			for(loop1= iBegin +1;loop1<loop+1;loop1++)
			{
				jBegin =loop1-m_aiDiag[loop1]+m_aiDiag[loop1-1]+1;
				if(jBegin < iBegin) jBegin = iBegin;
				iAdd_jj=m_aiDiag[loop1];
				dBuf=0.0;	
				for(loop2=jBegin;loop2<loop1;loop2++)
				{
					if(m_adValue[m_aiDiag[loop2]]==0.0){
						AfxMessageBox("Error in GK", MB_OK, 0 );
						return false;
					}
					dBuf+=m_adValue[iAdd_ii -loop+loop2]*m_adValue[iAdd_jj-loop1+loop2]
							/m_adValue[m_aiDiag[loop2]];
				}
				m_adValue[iAdd_ii-loop+loop1]-=dBuf;
			}
		}
		m_bDecomposed=true;
	}

	if(m_adValue[0]==0.0){
		AfxMessageBox("Error in GK", MB_OK, 0 );
		return false;
	}
	adB[0]=adB[0]/m_adValue[0];
	for(loop=1;loop<nRow;loop++)
	{
		dBuf=0.0;
		iBegin=loop-m_aiDiag[loop]+m_aiDiag[loop-1]+1;	
		iAdd_ii=m_aiDiag[loop];
		for(loop1=iBegin;loop1<loop;loop1++)
		{
			dBuf+=adB[loop1]*m_adValue[iAdd_ii-loop+loop1];
		}
		adB[loop]-=dBuf;
		if(m_adValue[iAdd_ii]==0.0){
			AfxMessageBox("Error in GK", MB_OK, 0 );
			return false;
		}
		adB[loop]/=m_adValue[iAdd_ii];
	}

	nBandWidth=1;	iRow=nRow-1;
	for(loop=nRow-2;loop>=0;loop--)
	{
		nBandWidth++;	
		while(nBandWidth>(m_aiDiag[iRow]-m_aiDiag[iRow-1])){
			nBandWidth--;
			iRow--;
		}
		dBuf=0.0;
		iAdd_ii=m_aiDiag[loop];
		for(loop1=1;loop1<nBandWidth;loop1++)
		{
			if(loop1<(m_aiDiag[loop1+loop]-m_aiDiag[loop1+loop-1]))
			{
				if(m_adValue[iAdd_ii]==0.0){
					AfxMessageBox("Error in GK", MB_OK, 0 );
					return false;
				}
				iAdd_kk=m_aiDiag[loop1+loop];
				dBuf+=m_adValue[iAdd_kk-loop1]*adB[loop1+loop]/m_adValue[iAdd_ii];
			}
		}
		adB[loop]-=dBuf;
	}
	return true;
}

bool CSparseMatrix::Decomposition(int nRow)
{
	int loop,loop1,loop2;
	int iBegin,jBegin,iAdd_ii,iAdd_jj;
	double dBuf;

	if(m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while decomposing";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		return false;
	}

	if(m_bDecomposed==false){
		for(loop=1;loop<nRow;loop++)
		{
			iBegin=loop-m_aiDiag[loop]+m_aiDiag[loop-1]+1;
			iAdd_ii=m_aiDiag[loop];
			for(loop1= iBegin +1;loop1<loop+1;loop1++)
			{
				jBegin =loop1-m_aiDiag[loop1]+m_aiDiag[loop1-1]+1;
				if(jBegin < iBegin) jBegin = iBegin;
				iAdd_jj=m_aiDiag[loop1];
				dBuf=0.0;	
				for(loop2=jBegin;loop2<loop1;loop2++)
				{
					if(m_adValue[m_aiDiag[loop2]]==0.0){
						AfxMessageBox("Error in GK", MB_OK, 0 );
						return false;
					}
					dBuf+=m_adValue[iAdd_ii -loop+loop2]*m_adValue[iAdd_jj-loop1+loop2]
							/m_adValue[m_aiDiag[loop2]];
				}
				m_adValue[iAdd_ii-loop+loop1]-=dBuf;
			}
		}
		m_bDecomposed=true;
	}
	return true;
}

double& CSparseMatrix::operator()(int iRow, int iCol)
{
	unsigned long iBuf;
	unsigned long nBandWidth;

	if(m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while while doing operator()";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		m_dBuf=0.0;
		return m_dBuf;
	}

	if(iRow>m_nRow){
		CString sText;
		sText="The row is out of range, while doing operator()";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		m_dBuf=0.0;
		return m_dBuf;
	}
	if(iCol>m_nRow){
		CString sText;
		sText="The col is out of range, while operator()";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		m_dBuf=0.0;
		return m_dBuf;
	}

	if(iRow>=iCol){
		if(iRow==0){
			nBandWidth=1;
		}
		else{
			nBandWidth=m_aiDiag[iRow]-m_aiDiag[iRow-1];
		}
		iBuf=iRow-iCol;
		if(nBandWidth<=iBuf){
			m_dBuf=0.0;
			return m_dBuf;
		}
		else{
			return m_adValue[m_aiDiag[iRow]-iBuf];
		}
	}
	else{
		nBandWidth=m_aiDiag[iCol]-m_aiDiag[iCol-1];
		iBuf=iCol-iRow;
		if(nBandWidth<=iBuf){
			m_dBuf=0.0;
			return m_dBuf;
		}
		else{
			return m_adValue[m_aiDiag[iCol]-iBuf];
		}
	}
}

bool CSparseMatrix::LdltSolve(int nRow,CMatrix& m)
{
	int loop1a;
	unsigned long loop,loop1,loop2;
	unsigned long  iBegin,jBegin,iAdd_ii,iAdd_jj,iAdd_kk;
	unsigned long  iRow,nBandWidth,nMatRow,nMatCol;
	double 	dBuf;

	if(m_bElementBufEmpty){
		CString sText;
		sText="The element buffer of SparseMatrix has not been allocated, while ldlt solving";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		return false;
	}
	nMatRow=m.GetRow();
	nMatCol=m.GetCol();
	if(nRow>nMatRow){
		CString sText;
		sText="The row of Matrix is not enough, while ldlt solving";
		AfxMessageBox( sText, MB_OK, 0 );
		//cout<<sText<<endl;
		return false;
	}

	if(m_bDecomposed==false){
		for(loop=1;loop<nRow;loop++)
		{
			iBegin=loop-m_aiDiag[loop]+m_aiDiag[loop-1]+1;
			iAdd_ii=m_aiDiag[loop];
			for(loop1= iBegin +1;loop1<loop+1;loop1++)
			{
				jBegin =loop1-m_aiDiag[loop1]+m_aiDiag[loop1-1]+1;
				if(jBegin < iBegin) jBegin = iBegin;
				iAdd_jj=m_aiDiag[loop1];
				dBuf=0.0;	
				for(loop2=jBegin;loop2<loop1;loop2++)
				{
					if(m_adValue[m_aiDiag[loop2]]==0.0){
						AfxMessageBox("Error in GK", MB_OK, 0 );
						return false;
					}
					dBuf+=m_adValue[iAdd_ii -loop+loop2]*m_adValue[iAdd_jj-loop1+loop2]
							/m_adValue[m_aiDiag[loop2]];
				}
				m_adValue[iAdd_ii-loop+loop1]-=dBuf;
			}
		}
		m_bDecomposed=true;
	}

	if(m_adValue[0]==0.0){
		AfxMessageBox("Error in GK", MB_OK, 0 );
		return false;
	}

	for(loop=0;loop<nMatCol;loop++){
		m(0,loop)=m(0,loop)/m_adValue[0];
		for(loop1=1;loop1<nRow;loop1++)
		{
			dBuf=0.0;
			iBegin=loop1-m_aiDiag[loop1]+m_aiDiag[loop1-1]+1;	
			iAdd_ii=m_aiDiag[loop1];
			for(loop2=iBegin;loop2<loop1;loop2++)
			{
				dBuf+=m(loop2,loop)*m_adValue[iAdd_ii-loop1+loop2];
			}
			m(loop1,loop)-=dBuf;
			if(m_adValue[iAdd_ii]==0.0){
				AfxMessageBox("Error in GK", MB_OK, 0 );
				return false;
			}
			m(loop1,loop)/=m_adValue[iAdd_ii];
		}

		nBandWidth=1;	iRow=nRow-1;
		for(loop1a=nRow-2;loop1a>=0;loop1a--)
		{
			nBandWidth++;	
			while(nBandWidth>(m_aiDiag[iRow]-m_aiDiag[iRow-1])){
				nBandWidth--;
				iRow--;
			}
			dBuf=0.0;
			iAdd_ii=m_aiDiag[loop1a];
			for(loop2=1;loop2<nBandWidth;loop2++)
			{
				if(loop2<(m_aiDiag[loop2+loop1a]-m_aiDiag[loop2+loop1a-1]))
				{
					if(m_adValue[iAdd_ii]==0.0){
						AfxMessageBox("Error in GK", MB_OK, 0 );
						return false;
					}
					iAdd_kk=m_aiDiag[loop2+loop1a];
					dBuf+=m_adValue[iAdd_kk-loop2]*m(loop2+loop1a,loop)/m_adValue[iAdd_ii];
				}
			}
			m(loop1a,loop)-=dBuf;
		}
	}
	return true;
}