code9.cpp

Eigenvalues and eigenvectors

//Code9.cpp
#include "Code9.h"

CMyWinApp MyApplication;

BOOL CMyWinApp::InitInstance()
{
    CCode9* pFrame = new CCode9;
    m_pMainWnd = pFrame;
    pFrame->ShowWindow(SW_SHOW);
    pFrame->UpdateWindow();
    return TRUE;
}

BEGIN_MESSAGE_MAP(CCode9,CFrameWnd)
	ON_WM_PAINT()
    ON_BN_CLICKED(IDC_BUTTON,OnButton)
END_MESSAGE_MAP()

CCode9::CCode9()
{
	int i,j;
	Create(NULL,"Code9: Eigenvalues and Eigenvectors",WS_OVERLAPPEDWINDOW,CRect(0,0,800,635),NULL);
	Arial80.CreatePointFont(80,"Arial");
	idc=301; fStatus=0; fMenu=0;
	table[1].hm=CPoint(20,310); 
	table[2].hm=CPoint(400,310);
	input[0][0].hm=CPoint(50,70);
	for (i=1;i<=M;i++)
		for (j=1;j<=M;j++)
			input[i][j].hm=CPoint(input[0][0].hm.x+(j-1)*80,input[0][0].hm.y+(i-1)*25);
	input[nInputItems][1].hm=CPoint(input[0][0].hm.x+100,input[0][0].hm.y+nInputItems*25);
	for (i=1;i<=M;i++)
		for (j=1;j<=M;j++)
			input[i][j].ed.Create(WS_CHILD | WS_VISIBLE | WS_BORDER, 
					CRect(CPoint(input[i][j].hm.x,input[i][j].hm.y),
					CSize(70,20)),this,idc++);
	input[nInputItems][1].ed.Create(WS_CHILD | WS_VISIBLE | WS_BORDER, 
		CRect(CPoint(input[0][0].hm.x+450,input[0][0].hm.y-50),
		CSize(100,20)),this,idc++);
	btn.Create("Compute",WS_CHILD | WS_VISIBLE | BS_DEFPUSHBUTTON,
		CRect(CPoint(input[0][0].hm.x,input[0][0].hm.y-50),CSize(100,30)),this,IDC_BUTTON);
}

CCode9::~CCode9()
{	
}

void CCode9::OnPaint()
{
	CPaintDC dc(this);
	dc.SelectObject(Arial80);
	dc.TextOut(input[0][0].hm.x+400,input[0][0].hm.y-50,"Epsilon");
	if (fStatus)
	{
		dc.TextOut(table[1].hm.x,table[1].hm.y-20,
			((fMenu==1)?"Power Method for the Most Dominant Eigenvalue":
			"QR Method for a symmetric matrix"));
		dc.TextOut(table[2].hm.x,table[2].hm.y-20,
			((fMenu==1)?"Shifted Power Method for the Least Dominant Eigenvalue":
			"Eigenvector in QR Method"));
	}
}

void CCode9::OnButton()
{
	int i,j,k;
	double **A,**B,**I;
	A=new double *[M+1];
	B=new double *[M+1];
	I=new double *[M+1];
	for (i=1;i<=M;i++)
	{
		A[i]=new double [M+1];
		B[i]=new double [M+1];
		I[i]=new double [M+1];
	}
	for (i=1;i<=M;i++)
		for (j=1;j<=M;j++)
			input[i][j].ed.GetWindowText(input[i][j].item);
	input[nInputItems][1].ed.GetWindowText(input[nInputItems][1].item);
	
	// determine the actual size of matrix
	for (i=1;i<=M;i++)
		if (input[i][i].item=="")
		{
			m=i-1; break;
		}
	
	// read the contents of A
	for (i=1;i<=m;i++)
		for (j=1;j<=m;j++)
		{
			I[i][j]=0;
			if (i==j)
				I[i][j]=1;
			A[i][j]=atof(input[i][j].item);
		}

	// test for matrix symmetry
	k=0;
	for (i=1;i<=m;i++)
		for (j=1;j<=m;j++)
			if (A[i][j]==A[j][i])
				k++;
	fMenu=((k==m*m)?2:1);
	if (fMenu==1)
	{
		PowerMtd(1,A);
		PowerTable(1); 
		if (StopA!=-1)
		{
			for (i=1;i<=m;i++)
				for (j=1;j<=m;j++)
					B[i][j]=A[i][j]-eigen[StopA].lambda*I[i][j];
			PowerMtd(2,B);
			PowerTable(2);
		}
	}
	if (fMenu==2)
	{
		QRMtd(A);
		QRTable(1);
		QRTable(2);
	}
}

void CCode9::PowerMtd(int fPower,double **c)
{
	double u[M+1],w[M+1],lamb[maxIter+1],error[maxIter+1];
	int i,j,k;
	for (i=1;i<=m;i++)
	{
		w[i]=0;
		if (i==2)
			w[i]=1;
		eigen[0].v[i]=w[i];
		eigen[0].vB[i]=w[i];
	}
	epsilon=atof(input[nInputItems][1].item);
	for (k=0;k<=maxIter;k++)
	{
		for (i=1;i<=m;i++)
		{
			u[i]=0;
			for (j=1;j<=m;j++)
				u[i] += c[i][j]*w[j];
		}
		lamb[k]=u[1];
		for (i=1;i<=m;i++)
			if (fabs(lamb[k])<=fabs(u[i]))
			{
				lamb[k]=u[i];
				if (fPower==1)
					eigen[k+1].lambda=lamb[k];
				if (fPower==2)
					eigen[k+1].lambdaB=lamb[k];
			}
		for (i=1;i<=m;i++)
		{
			w[i]=u[i]/lamb[k];
			if (fPower==1)
				eigen[k+1].v[i]=w[i];
			if (fPower==2)
				eigen[k+1].vB[i]=w[i];
		}
		if (k>0)
		{
			error[k]=lamb[k]-lamb[k-1];
			if (fPower==1)
				eigen[k].error=error[k];
			if (fPower==2)
				eigen[k].errorB=error[k];
			if (fabs(error[k])<epsilon)
			{
				if (fPower==1)
				{
					StopA=k;
					eigen[StopA].lambda=lamb[StopA];
					fStatus=1;
					Invalidate();
					break;
				}
				if (fPower==2)
				{
					StopB=k;
					eigen[StopB].lambdaB=lamb[StopB];
					break;
				}
			}
			else
			{
				if (k==maxIter && fPower==1)
					StopA=-1;
				if (k==maxIter && fPower==2)
					StopB=-1;
			}
		}
	}			
}

void CCode9::QRMtd(double **c)
{
	int i,j,k,u,p;
	double a[maxIter+1][M+1][M+1],r[M+1][M+1],q[M+1][M+1],s[M+1][M+1],w[M+1][M+1];
	double t[M+1],P[M+1][M+1],I[M+1][M+1],B[M+1][M+1];
	double Alpha,R,Sum;
	epsilon=atof(input[nInputItems][1].item);
	for (i=1;i<=m;i++)
		for (j=1;j<=m;j++)
		{
			r[i][j]=0;
			a[0][i][j]=c[i][j];
			I[i][j]=0;
			if (i==j)
				I[i][j]=1;
		}
	//Householder Method
	for (k=1;k<=m-2;k++)
	{
		Sum=0;
		for (i=k+1;i<=m;i++)
			Sum+=pow(a[0][i][k],2);
		if (a[0][k+1][k]<0)
			Alpha=sqrt(Sum);
		else
			Alpha=(-1)*sqrt(Sum);
		R=sqrt((pow(Alpha,2)*1/2)-(Alpha*a[0][k+1][k]*1/2));
		for (i=1;i<=k;i++)
			t[i]=0;
		t[k+1]=(a[0][k+1][k]-Alpha)/(2*R);
		for (i=k+2;i<=m;i++)
			t[i]=a[0][i][k]/(2*R);
		for (i=1;i<=m;i++)
			for (j=1;j<=m;j++)
				P[i][j]=I[i][j]-(2*t[i]*t[j]);
		for (i=1;i<=m;i++)
			for (j=1;j<=m;j++)
			{
				Sum=0;
				for (u=1;u<=m;u++)
					Sum+=P[i][u]*a[0][u][j];
				B[i][j]=Sum;
			}
		for (i=1;i<=m;i++)
			for (j=1;j<=m;j++)
			{
				Sum=0;
				for (u=1;u<=m;u++)
					Sum += B[i][u]*P[u][j];
				a[0][i][j]=Sum;
			}
	}
	for (i=1;i<=m;i++)
		for (j=1;j<=m;j++)
			q[i][j]=a[0][i][j];
	//QR algorithm
	StopA=-1;
	for (p=0;p<=maxIter;p++)
	{
		for (i=1;i<=m;i++)
			for (j=i;j<=m;j++)
			{
				if (i==j)
				{
					Sum=0;
					for (k=1;k<=m;k++)
						Sum+=pow(q[k][i],2);
					r[i][j]=pow(Sum,0.5);
					for (k=1;k<=m;k++)
						q[k][j]=(q[k][j]/r[i][j]);
				}
				if (i<j)
				{
					Sum=0;
					for (k=1;k<=m;k++)
						Sum+=(q[k][i]*q[k][j]);
					r[i][j]=Sum;
					for (k=1;k<=m;k++)
						q[k][j]=q[k][j]-(r[i][j]*q[k][i]);
				}
			}
		if (p==0)
			for (i=1;i<=m;i++)
				for (j=1;j<=m;j++)
				{
					s[i][j]=q[i][j];
					eigen[p].vQR[i][j]=q[j][i];
				}
		else
		{
			for (i=1;i<=m;i++)
				for (j=1;j<=m;j++)
				{
					Sum=0;
					for (k=1;k<=m;k++)
						Sum+=s[i][k]*q[k][j];
					w[i][j]=Sum;
				}
			for (i=1;i<=m;i++)
				for (j=1;j<=m;j++)
				{
					s[i][j]=w[i][j];
					eigen[p].vQR[i][j]=w[j][i];
				}
		}
		for (i=1;i<=m;i++)
			for (j=1;j<=m;j++)
			{
				Sum=0;
				for (k=1;k<=m;k++)
					Sum+=r[i][k]*q[k][j];
				a[p][i][j]=Sum;
				q[i][j]=a[p][i][j];
				eigen[p].lambdaQR[i]=a[p][i][i];
			}
		if (p>0)
		{
			j=0;
			for (i=1;i<=m;i++)
			{
				eigen[p].errorQR[i]=eigen[p].lambdaQR[i]-eigen[p-1].lambdaQR[i];
				if (fabs(eigen[p].errorQR[i])<epsilon)
					j++;
			}
			if (j==m)
			{
				fStatus=1;
				Invalidate();
				StopA=p;
				break;
			}
		}
	}

}

void CCode9::PowerTable(int c)
{
	int Stop,i,j,k;
	CString str;
	CRect rcTable[3];
	rcTable[c]=CRect(table[c].hm.x,table[c].hm.y,table[c].hm.x+370,table[c].hm.y+290);
	table[c].list.DestroyWindow();
    table[c].list.Create(WS_VISIBLE | WS_CHILD | WS_DLGFRAME | LVS_REPORT
        | LVS_NOSORTHEADER,rcTable[c],this,idc++);
	table[c].list.InsertColumn(0,"i",LVCFMT_CENTER,25);
	table[c].list.InsertColumn(1,((c==1)?"lambda hi":"lambda lo"),LVCFMT_CENTER,70);
	for (i=2;i<=m+1;i++)
	{
		str.Format("v[%d]",i-1);
		table[c].list.InsertColumn(i,str,LVCFMT_CENTER,70);
	}
	table[c].list.InsertColumn(m+2,"error",LVCFMT_CENTER,70);
	Stop=((c==1)?StopA:StopB);
	for (k=0;k<=Stop;k++)
	{
		str.Format("%d",k); 
		table[c].list.InsertItem(k,str,0);
		if (k>0)
		{
			str.Format("%lf",((c==1)?eigen[k].lambda:eigen[k].lambdaB+eigen[StopA].lambda));
			table[c].list.SetItemText(k,1,str);
			str.Format("%lf",((c==1)?eigen[k].error:eigen[k].errorB));
			table[c].list.SetItemText(k,m+2,str);
		}
		for (j=1;j<=m;j++)
		{
			str.Format("%lf",((c==1)?eigen[k].v[j]:eigen[k].vB[j])); 
			table[c].list.SetItemText(k,j+1,str);
		}
	}
}

void CCode9::QRTable(int c)
{
	int Stop,i,j,k,p;
	CString str;
	CRect rcTable[4];
	rcTable[c]=CRect(table[c].hm.x,table[c].hm.y,table[c].hm.x+370,table[c].hm.y+290);
	table[c].list.DestroyWindow();
    table[c].list.Create(WS_VISIBLE | WS_CHILD | WS_DLGFRAME | LVS_REPORT
        | LVS_NOSORTHEADER,rcTable[c],this,idc++);
	table[c].list.InsertColumn(0,"i",LVCFMT_CENTER,25);
	if (c==1)
		for (i=1;i<=m;i++)
		{
			str.Format("lambda[%d]",i);
			table[c].list.InsertColumn(i,str,LVCFMT_CENTER,70);
			str.Format("error[%d]",i);
			table[c].list.InsertColumn(m+i,str,LVCFMT_CENTER,70);
		}
	if (c==2)
	{
		k=0;
		for (i=1;i<=m;i++)
			for(j=1;j<=m;j++)
			{
				str.Format("v[%d][%d]",i,j);
				table[c].list.InsertColumn(i+k,str,LVCFMT_CENTER,70);
				k++;
			}
	}
	Stop=StopA;
	for (k=0;k<=Stop;k++)
	{
		str.Format("%d",k); 
		table[c].list.InsertItem(k,str,0);
		if (k>0)
		{
			if (c==1)
				for (i=1;i<=m;i++)
				{
					str.Format("%lf",eigen[k].lambdaQR[i]);
					table[c].list.SetItemText(k,i,str);
					str.Format("%lf",eigen[k].errorQR[i]);
					table[c].list.SetItemText(k,m+i,str);
				}
			if (c==2)
			{
				p=0;
				for (i=1;i<=m;i++)
				{
					for (j=1;j<=m;j++)
					{
						str.Format("%lf",eigen[k].vQR[i][j]);
						table[c].list.SetItemText(k,j+p,str);
					}
					p += m;
				}
			}
		}
	}
}