jacobi

雅克比算法

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]);
		}
	}
}