子空间迭代法（包含雅克比法）.cpp

子空间迭代法的C++以及Fortran源程序

void CGlobalElement::ModeAnalysis()
{
	bool bLargerThanError;
	int loop,loop1;
	int nTotalDOF,nFreeDOF,nMode,nLumpMass;
	CUIntArray aiBuf;
	double dBuf,dBuf1,dError,dRatioGKtoMass,dQm;
	CArray <double,double&> adBuf,adEigenBak; 
	CMatrix matBuf,matBuf1;
	CMatrix matV,matU,matA,matB;
	CSparseMatrix smatGKBak;
	CSparseMatrix smatMass;

	if(m_apEle.GetSize()==0) return;
	m_bFlagModeAnalyzed=true;
	dError=1.0e-5;

	nTotalDOF=m_Node.GetTotalDOF();
	nFreeDOF=m_Node.GetFreeDOF();

	m_adDampRatio.SetSize(m_nMode);
	m_smatGK.ElementBufRealloc();
	m_smatGK=0.0;

	GetMassMatrix(smatMass);

	for(loop=0;loop<m_nEle;loop++){
		m_apEle[loop]->StiffAssemble(m_smatGK);
	}

	dBuf=0.0;
	for(loop=0;loop<nFreeDOF;loop++)
		dBuf+=m_smatGK(loop,loop);
	dBuf/=nFreeDOF;
	dBuf1=0.0;
	for(loop=0;loop<nFreeDOF;loop++){
		dBuf1+=smatMass(loop,loop);
	}
	if(m_iMassMatrixType==CONSISTENT_MASS){
		dBuf1/=nFreeDOF;
	}
	else{
		nLumpMass=m_LumpMass.m_apLoad.GetSize();
		dBuf1/=nLumpMass*2.0;
	}
	dRatioGKtoMass=dBuf/dBuf1;

	dBuf=1.0/dBuf;
	m_smatGK*=dBuf;

	dQm=1.0/dBuf1;
	smatMass*=dQm;

	aiBuf.SetSize(nFreeDOF);
	adBuf.SetSize(nFreeDOF);
	dBuf=0.0;
	for(loop=0;loop<nFreeDOF;loop++){
		adBuf[loop]=-1.0;
		if(smatMass(loop,loop)>0.0){
			adBuf[loop]=m_smatGK(loop,loop)/smatMass(loop,loop);
		}
		if(dBuf<adBuf[loop]) dBuf=adBuf[loop];	
	}

	for(loop=0;loop<nFreeDOF;loop++){
		if(adBuf[loop]<0.0) adBuf[loop]=dBuf+100.0;	
	}

	aiBuf[0]=0;
	for(loop=0;loop<nFreeDOF;loop++){
		for(loop1=0;loop1<loop;loop1++){
			if(adBuf[aiBuf[loop1]]>adBuf[loop]){
				aiBuf.InsertAt(loop1,loop);
				break;
			}
		}
		if(loop1==loop){
			aiBuf[loop1]=loop;
		}
	}

	nMode=m_nMode+8;
	if(2*m_nMode<nMode) nMode=2*m_nMode;
	if(nMode>nFreeDOF)	nMode=nFreeDOF;
	if(m_iMassMatrixType==LUMP_MASS&&nMode>2*nLumpMass) nMode=2*nLumpMass;
	if(m_nMode>nMode)	m_nMode=nMode;

	matV.Realloc(nFreeDOF,nMode);
	matU.Realloc(nFreeDOF,nMode);
	matA.Realloc(nMode,nMode);
	matB.Realloc(nMode,nMode);

	m_adEigen.SetSize(nMode);
	adEigenBak.SetSize(nMode);
	for(loop=0;loop<nMode;loop++)
		m_adEigen[loop]=0.0;

	matV=0.0;
	for(loop=0;loop<nFreeDOF;loop++)
		matV(loop,0)=1.0;
	for(loop=1;loop<nMode;loop++)
		matV(aiBuf[loop],loop)=1.0;

	smatGKBak=m_smatGK;
	do{
		bLargerThanError=false;
		for(loop=0;loop<nMode;loop++)
			adEigenBak[loop]=m_adEigen[loop];
		matU=smatMass*matV;

		if(!m_smatGK.LdltSolve(nFreeDOF,matU))	return;
		matBuf=smatGKBak*matU;
		matBuf1.Trans(matU);
		matA=matBuf1*matBuf;
		matBuf=smatMass*matU;
		matB=matBuf1*matBuf;

		JacobiEigen(matA,matB);
		matV=matU*matB;
		for(loop=0;loop<nMode;loop++){
			dBuf=fabs((m_adEigen[loop]-adEigenBak[loop])/m_adEigen[loop]);
			if(dBuf>dError){
				bLargerThanError=true;
				break;
			}
		}
	}while(bLargerThanError);

	dBuf=sqrt(dQm);
	m_matEigenVector.Realloc(nTotalDOF,m_nMode);
	m_matEigenVector=0.0;
	for(loop=0;loop<m_nMode;loop++)
	for(loop1=0;loop1<nFreeDOF;loop1++)
		m_matEigenVector(loop1,loop)=matV(loop1,loop)*dBuf;

	for(loop=0;loop<m_nMode;loop++){
		m_adEigen[loop]*=dRatioGKtoMass;
		m_adEigen[loop]=sqrt(m_adEigen[loop]);
	}
	m_smatGK.ElementBufEmpty();
}

void CGlobalElement::JacobiEigen(CMatrix &matA, CMatrix &matB)
{
	CMatrix matBBak,matR,matRT,matV;
	double dAlfa,dBeta,dBuf,dBuf1;
	double dA,dB,dC;
	double dError;
	int nRow;
	int iRow,iCol;
	int loop,loop1;
	CArray<double,double&> adBuf;
	CUIntArray aiBuf;

	nRow=matA.GetRow();
	matBBak.Realloc(nRow,nRow);
	matR.Realloc(nRow,nRow);
	matRT.Realloc(nRow,nRow);
	matV.Realloc(nRow,nRow);
	matBBak=matB;

	matV=0.0;
	for(loop=0;loop<nRow;loop++)
		matV(loop,loop)=1.0;

	dError=1.0e-24;

	dBuf=0.0;
	for(loop=0;loop<nRow;loop++){
		for(loop1=loop+1;loop1<nRow;loop1++){
			dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1));
			if(dBuf<fabs(dBuf1)){
				dBuf=fabs(dBuf1);
				iRow=loop;iCol=loop1;
			}
			dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1));
			if(dBuf<fabs(dBuf1)){
				dBuf=fabs(dBuf1);
				iRow=loop;iCol=loop1;
			}
		}
	}

	while(dBuf>dError){
		matR=0.0;
		for(loop=0;loop<nRow;loop++) matR(loop,loop)=1.0;
		dA=matA(iRow,iRow)*matB(iRow,iCol)-matB(iRow,iRow)*matA(iRow,iCol);
		dC=-matA(iCol,iCol)*matB(iRow,iCol)+matB(iCol,iCol)*matA(iRow,iCol);
		if(dA==0.0&&dC==0.0){
			dAlfa=0.0;
			dBeta=-matA(iRow,iCol)/matA(iCol,iCol);
		}
		else{
			dB=matA(iRow,iRow)*matB(iCol,iCol)-matA(iCol,iCol)*matB(iRow,iRow);
			if(dB>=0.0){
				dBuf=dB/2.0;
				dAlfa=-dC/(dBuf+sqrt(dBuf*dBuf-dA*dC));
			}
			else{
				dBuf=dB/2.0;
				dAlfa=-dC/(dBuf-sqrt(dBuf*dBuf-dA*dC));
			}
			dBeta=dA/dC*dAlfa;
		}
		matR(iRow,iCol)=dAlfa;
		matR(iCol,iRow)=dBeta;
		matRT.Trans(matR);

		matA=matRT*matA*matR;
		matB=matRT*matB*matR;
		matV*=matR;

		dBuf=0.0;
		for(loop=0;loop<nRow;loop++){
			for(loop1=loop+1;loop1<nRow;loop1++){
				dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1));
				if(dBuf<fabs(dBuf1)){
					dBuf=fabs(dBuf1);
					iRow=loop;iCol=loop1;
				}
				dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1));
				if(dBuf<fabs(dBuf1)){
					dBuf=fabs(dBuf1);
					iRow=loop;iCol=loop1;
				}
			}
		}
	}

	matRT.Trans(matV);
	matR=matRT*matBBak;
	for(loop=0;loop<nRow;loop++){
		dBuf=0.0;
		for(loop1=0;loop1<nRow;loop1++){
			dBuf+=matR(loop,loop1)*matV(loop1,loop);
		}
		if(dBuf<0.0) dBuf=-dBuf;
		dBuf=sqrt(dBuf);
		for(loop1=0;loop1<nRow;loop1++)
			matV(loop1,loop)/=dBuf;
	}

	adBuf.SetSize(nRow);
	for(loop=0;loop<nRow;loop++){
		adBuf[loop]=matA(loop,loop)/matB(loop,loop);
	}
	aiBuf.SetSize(nRow);
	aiBuf[0]=0;
	for(loop=0;loop<nRow;loop++){
		for(loop1=0;loop1<loop;loop1++){
			if(fabs(adBuf[aiBuf[loop1]])>fabs(adBuf[loop])){
				aiBuf.InsertAt(loop1,loop);
				break;
			}
		}
		if(loop1==loop){
			aiBuf[loop1]=loop;
		}
	}

	for(loop=0;loop<nRow;loop++){
		m_adEigen[loop]=adBuf[aiBuf[loop]];
	}
	for(loop=0;loop<nRow;loop++){
		for(loop1=0;loop1<nRow;loop1++){
			matB(loop1,loop)=matV(loop1,aiBuf[loop]);
		}
	}
}