code5.cpp

this program is system of linear equations

#include "Code5.h"

BOOL CMyWinApp::InitInstance()
{
	AfxEnableControlContainer();
	CCode5 dlg;
	m_pMainWnd = &dlg;
	dlg.DoModal();
	return FALSE;
}

CMyWinApp MyApplication;

BEGIN_MESSAGE_MAP(CCode5,CDialog)
	ON_BN_CLICKED(IDC_GAUSS, OnGauss)
	ON_BN_CLICKED(IDC_CROUT, OnCrout)
	ON_BN_CLICKED(IDC_CHOLESKY, OnCholesky)
	ON_BN_CLICKED(IDC_THOMAS, OnThomas)
	ON_BN_CLICKED(IDC_GAUSSSEIDEL, OnGaussSeidel)
	ON_BN_CLICKED(IDC_RESET, OnReset)
END_MESSAGE_MAP()

CCode5::CCode5(CWnd* pParentWnd) 
: CDialog (IDD_DIALOG1, pParentWnd)
{
	x=new double [N+1];
	b=new double [N+1];
	sx=new CString [N+1];
	sb=new CString [N+1];
	A=new double *[N+1];
	L=new double *[N+1];
	U=new double *[N+1];
	sA=new CString *[N+1];
	for (int i=0;i<=N;i++)
	{
		A[i]=new double [N+1];
		L[i]=new double [N+1];
		U[i]=new double [N+1];
		sA[i]=new CString [N+1];
	}
}

CCode5::~CCode5()
{
	for (int i=0;i<=N;i++)
		delete A[i],L[i],U[i],sA[i];
	delete x,b,sx,sb,A,L,U,sA;
}

void CCode5::OnOK()
{
	CDialog::OnOK();
}

void CCode5::DisplayResults()
{
	for (int i=1;i<=N;i++)
			sx[i].Format("%lf",x[i]);
	SetDlgItemText(IDC_x1,((1<=n)?sx[1]:""));
	SetDlgItemText(IDC_x2,((2<=n)?sx[2]:""));
	SetDlgItemText(IDC_x3,((3<=n)?sx[3]:""));
	SetDlgItemText(IDC_x4,((4<=n)?sx[4]:""));
	SetDlgItemText(IDC_x5,((5<=n)?sx[5]:""));
	SetDlgItemText(IDC_MESSAGE,method+" is successful.");
}

void CCode5::DisplayError(CString condition)
{
	SetDlgItemText(IDC_x1,"");
	SetDlgItemText(IDC_x2,"");
	SetDlgItemText(IDC_x3,"");
	SetDlgItemText(IDC_x4,"");
	SetDlgItemText(IDC_x5,"");
	SetDlgItemText(IDC_MESSAGE,
		method+" fails, the matrix is not "+condition+".");
}

void CCode5::ReadInput()
{
    GetDlgItemText(IDC_A11,sA[1][1]);
	GetDlgItemText(IDC_A12,sA[1][2]);
	GetDlgItemText(IDC_A13,sA[1][3]);
	GetDlgItemText(IDC_A14,sA[1][4]);
	GetDlgItemText(IDC_A15,sA[1][5]);
	GetDlgItemText(IDC_A21,sA[2][1]);
	GetDlgItemText(IDC_A22,sA[2][2]);
	GetDlgItemText(IDC_A23,sA[2][3]);
	GetDlgItemText(IDC_A24,sA[2][4]);
	GetDlgItemText(IDC_A25,sA[2][5]);
	GetDlgItemText(IDC_A31,sA[3][1]);
	GetDlgItemText(IDC_A32,sA[3][2]);
	GetDlgItemText(IDC_A33,sA[3][3]);
	GetDlgItemText(IDC_A34,sA[3][4]);
	GetDlgItemText(IDC_A35,sA[3][5]);
	GetDlgItemText(IDC_A41,sA[4][1]);
	GetDlgItemText(IDC_A42,sA[4][2]);
	GetDlgItemText(IDC_A43,sA[4][3]);
	GetDlgItemText(IDC_A44,sA[4][4]);
	GetDlgItemText(IDC_A45,sA[4][5]);
	GetDlgItemText(IDC_A51,sA[5][1]);
	GetDlgItemText(IDC_A52,sA[5][2]);
	GetDlgItemText(IDC_A53,sA[5][3]);
	GetDlgItemText(IDC_A54,sA[5][4]);
	GetDlgItemText(IDC_A55,sA[5][5]);
	GetDlgItemText(IDC_b1,sb[1]);
	GetDlgItemText(IDC_b2,sb[2]);
	GetDlgItemText(IDC_b3,sb[3]);
	GetDlgItemText(IDC_b4,sb[4]);
	GetDlgItemText(IDC_b5,sb[5]);
	
	for (int i=1;i<=N;i++)
		if (sA[i][i]=="")
		{
			n=i-1; break;
		}
	// convert the string input to double
	for (i=1;i<=n;i++)
	{
		for (int j=1;j<=n;j++)
			A[i][j]=atof(sA[i][j]);
		b[i]=atof(sb[i]);
	}
}

void CCode5::OnGaussSeidel()
{
	const int MAX=10;
	double Sum,error;
	double *xOld;
	int i,j,k;
	bool fStatus=1;
	CString condition="diagonally-dominant";
	method="Gauss-Seidel";
	xOld=new double [N+1];
	ReadInput();

	for (i=1;i<=n;i++)
	{
		Sum=0;
		for (j=1;j<=n;j++)
			if (i!=j)
				Sum += fabs(A[i][j]);
		if (fabs(A[i][i])<=Sum)
			fStatus=0;
	}

	for (i=1;i<=n;i++)
		x[i]=0;
	
	// perform iterations & update x
	for (k=0;k<=MAX;k++)
	{
		error=0;
		for (i=1;i<=n;i++)
		{
			Sum=0;
			for (j=1;j<=n;j++)
				if (i!=j)
					Sum += A[i][j]*x[j];
			xOld[i]=x[i];
			x[i]=(b[i]-Sum)/A[i][i];
			error=((error>fabs(x[i]-xOld[i]))?error:fabs(x[i]-xOld[i]));
		}
	}
	delete xOld;
	if (fStatus)
		DisplayResults();
	else
		DisplayError(condition);
}

void CCode5::OnGauss()
{
	int i,j,k;
	double Sum,p;
	bool fStatus=1;
	CString condition="invertible";
	method="Gauss";
	ReadInput();

	// perform row operations
	for (k=1;k<=n-1;k++)
		for (i=k+1;i<=n;i++)
		{
			if (A[k][k]==0)
			{
				fStatus=0;
				break;
			}
			else
				p=A[i][k]/A[k][k];
			for (j=1;j<=n;j++)
				A[i][j] -= p*A[k][j];
			b[i] -= p*b[k];
		}

	// perform backward substitutions
	for (i=n;i>=1;i--)
	{
    	Sum=0;
		x[i]=0;
		for (j=i;j<=n;j++)
         	Sum += A[i][j]*x[j];
		x[i]=(b[i]-Sum)/A[i][i];
	}
	if (fStatus)
		DisplayResults();
	else
		DisplayError(condition);
}

void CCode5::OnCrout()
{
	int i,j,k,r;
	double Sum,*w;
	bool fStatus=1;
	CString condition="invertible";
	method="Crout";
	w=new double [N+1];
	ReadInput();

	// initialize L and U
	for (i=1;i<=n;i++)
		for (j=1;j<=n;j++)
		{
			L[i][j]=U[i][j]=0;
			if (i==j)
				U[i][j]=1;
		}

	// compute L and U
	for (i=1;i<=n;i++)
		L[i][1]=A[i][1];
	for (j=2;j<=n;j++)
		if (L[1][1]==0)
		{
			fStatus=0;
			break;
		}
		else
			U[1][j]=A[1][j]/L[1][1];
	for (j=2; j<=n-1;j++)
	{
		for (i=j; i<=n; i++)
		{
			Sum=0;
			for (k=1;k<=j-1;k++)
				Sum += L[i][k]*U[k][j];
			L[i][j]=A[i][j]- Sum;
		}
		for (k=j+1;k<=n;k++)
		{
			Sum=0;
			for (r=1;r<=j-1;r++)
				Sum += L[j][r]*U[r][k];
			if (L[j][j]==0)
			{
				fStatus=0;
				break;
			}
			else
				U[j][k]=(A[j][k]-Sum)/L[j][j];
		}	
	}
	Sum=0;
	for (k=1;k<=n-1;k++)
		Sum += L[n][k]*U[k][n];
	L[n][n]=A[n][n]-Sum;

	// forward substitutions for finding w
	w[1]=b[1]/L[1][1];
	for (j=1;j<=n;j++)
	{
		Sum=0;
		for (k=1; k<=j-1;k++)
			Sum += L[j][k]*w[k];
		if (L[j][j]==0)
		{
			fStatus=0;
			break;
		}
		else
			w[j] =(b[j]-Sum)/L[j][j];
	}
	
	// backward substitutions for finding x
	for (i=n;i>=1;i--)
	{
		Sum=0;
		for (k=i+1; k<=n; k++)
			Sum += U[i][k]*x[k];
		x[i]=w[i]-Sum;
	}
	delete w;
	if (fStatus)
		DisplayResults();
	else
		DisplayError(condition);
}

void CCode5::OnCholesky()
{
	int i,j,k;
	double Sum,*w;
	bool fStatus=1;
	CString condition="positive-definite";
	method="Cholesky";
	w=new double [N+1];
	ReadInput();

	// check for positive-definite
	for (i=1;i<=n;i++)
		for (j=i;j<=n;j++)
			if ((A[i][j]!=A[j][i]) 
					|| A[i][i]<=0
					|| (i!=j && A[i][i]<=fabs(A[i][j])))
			{
				fStatus=0;
				break;
			}

	// initialize L and U
	for (i=1;i<=n;i++)
		for (j=1;j<=n;j++)
			L[i][j]=U[i][j]=0;

	// compute L and U
	for (k=1;k<=n;k++)
	{
		Sum=0;
		for (j=1;j<=k-1;j++)
			Sum += pow(L[k][j],2);
		L[k][k]=sqrt(A[k][k]-Sum);
		for (i=1;i<=k-1;i++)
		{
			Sum=0;
			for (j=1;j<=i-1;j++)
				Sum += L[i][j]*L[k][j];
			L[k][i]=(A[k][i]-Sum)/L[i][i];
		}
	}

	// forward substitutions for finding w
	//w[1]=b[1]/L[1][1];
	for (j=1;j<=n;j++)
	{
		Sum=0;
		for (k=1; k<=j-1;k++)
			Sum += L[j][k]*w[k];
		w[j] =(b[j]-Sum)/L[j][j];
	}
	
	// backward substitutions for finding x
	for (i=n;i>=1;i--)
	{
		Sum=0;
		for (k=i+1; k<=n; k++)
			Sum += U[i][k]*x[k];
		x[i]=w[i]-Sum;
	}
	delete w;
	if (fStatus)
		DisplayResults();
	else
		DisplayError(condition);
}

void CCode5::OnThomas()
{
	int i,j;
	bool fStatus=1;
	double *e,*f,*g;
	CString condition="tridiagonal";
	method="Thomas";
	e=new double [N+1];
	f=new double [N+1];
	g=new double [N+1];
	ReadInput();

	// check for tridiagonality
	for (i=1;i<=n;i++)
	{
		if (i==1)
			for (j=3;j<=n;j++)
				if (A[i][j]!=0)
					fStatus=0;
		if (i==n)
			for (j=1;j<=n-2;j++)
				if (A[i][j]!=0)
					fStatus=0;
		if (i>1 && i<n)
			if (A[i][i+2]!=0 && A[i][i-2]!=0)
				fStatus=0;
	}

	// compute e,f,g
	for (i=1;i<=n;i++)
	{
		f[i]=A[i][i];
		if (i<=n-1)
			g[i]=A[i][i+1];
		if (i>=2)
			e[i]=A[i][i-1];
	}
    for (i=2;i<=n;i++)
    {
		e[i]/=f[i-1];
		f[i] -= e[i]*g[i-1];
	}
    for (i=2;i<=n;i++)
		b[i] -= e[i]*b[i-1];
    x[n]=b[n]/f[n];
    for (i=n-1;i>=1;i--)
		x[i]=(b[i]-g[i]*x[i+1])/f[i];
	delete e,f,g;
	if (fStatus)
		DisplayResults();
	else
		DisplayError(condition);
}

void CCode5::OnReset()
{
	int i,j;
	for (i=1;i<=N;i++)
	{
		for (j=1;j<=N;j++)
			sA[i][j]="";
		sb[i]="";
		sx[i]="";
	}
    SetDlgItemText(IDC_A11,sA[1][1]);
	SetDlgItemText(IDC_A12,sA[1][2]);
	SetDlgItemText(IDC_A13,sA[1][3]);
	SetDlgItemText(IDC_A14,sA[1][4]);
	SetDlgItemText(IDC_A15,sA[1][5]);
	SetDlgItemText(IDC_A21,sA[2][1]);
	SetDlgItemText(IDC_A22,sA[2][2]);
	SetDlgItemText(IDC_A23,sA[2][3]);
	SetDlgItemText(IDC_A24,sA[2][4]);
	SetDlgItemText(IDC_A25,sA[2][5]);
	SetDlgItemText(IDC_A31,sA[3][1]);
	SetDlgItemText(IDC_A32,sA[3][2]);
	SetDlgItemText(IDC_A33,sA[3][3]);
	SetDlgItemText(IDC_A34,sA[3][4]);
	SetDlgItemText(IDC_A35,sA[3][5]);
	SetDlgItemText(IDC_A41,sA[4][1]);
	SetDlgItemText(IDC_A42,sA[4][2]);
	SetDlgItemText(IDC_A43,sA[4][3]);
	SetDlgItemText(IDC_A44,sA[4][4]);
	SetDlgItemText(IDC_A45,sA[4][5]);
	SetDlgItemText(IDC_A51,sA[5][1]);
	SetDlgItemText(IDC_A52,sA[5][2]);
	SetDlgItemText(IDC_A53,sA[5][3]);
	SetDlgItemText(IDC_A54,sA[5][4]);
	SetDlgItemText(IDC_A55,sA[5][5]);
	SetDlgItemText(IDC_b1,sb[1]);
	SetDlgItemText(IDC_b2,sb[2]);
	SetDlgItemText(IDC_b3,sb[3]);
	SetDlgItemText(IDC_b4,sb[4]);
	SetDlgItemText(IDC_b5,sb[5]);
	SetDlgItemText(IDC_x1,sx[1]);
	SetDlgItemText(IDC_x2,sx[2]);
	SetDlgItemText(IDC_x3,sx[3]);
	SetDlgItemText(IDC_x4,sx[4]);
	SetDlgItemText(IDC_x5,sx[5]);
}