mainangle.cpp

列主元gauss解方程。方程的大小可以对源代码稍稍修改即可做到。

#include <iostream.h>
#include <math.h>
#include <windows.h>


int FindMaxLine(float **A,int column,int Height);
void ExchangeLines(float **A,int Width,int lineA,int lineB,float *b);
void Angel(float **A,int column,int Width,float *b);
float SolveEquation(float **A,float *b,int xi,int Width);

void main()
{

	int N=0;
	float **A=NULL;
	float *b=NULL;
	int i=0;
	int ai=0,aj=0;
	int j=0;
	int temp;

//apply memory

	N = 9;
	A = new float *[N];
	for (i = 0;i < N; i++) 
	{
		A[i]=new float[N];
	}
	b=new float [N];
	for (ai=0;ai<N;ai++)
	{
		for (aj=0;aj<N;aj++)
		{
			A[ai][aj]=0;
		}
	}


	A[0][0]=31;    A[0][1]=-13;   A[0][5]=-10;
	A[1][0]=-13;   A[1][1]=35;    A[1][2]=-9;       A[1][4]=-11;
	A[2][1]=-9;    A[2][2]=31;    A[2][3]=-10;
	A[3][2]=-10;   A[3][3]=79;    A[3][4]=-30;      A[3][8]=-9;
	A[4][3]=-30;   A[4][4]=57;    A[4][5]=-7;       A[4][7]=-5;
	A[5][4]=-7;    A[5][5]=47;    A[5][6]=-30;
	A[6][5]=-30;   A[6][6]=41;
	A[7][4]=-5;    A[7][7]=27;    A[7][8]=-2;
	A[8][3]=-9;    A[8][7]=-2;    A[8][8]=29;

	b[0]=-15;b[1]=27;b[2]=-23;b[3]=0;b[4]=-20;b[5]=12;b[6]=-7;b[7]=7;b[8]=10;


/*

	A[0][0]=2;    A[0][1]=-4;    A[0][2]=6;
	A[1][0]=4;    A[1][1]=-9;    A[1][2]=2; 
	A[2][0]=1;    A[2][1]=-1;    A[2][2]=3;
	
	b[0]=3;       b[1]=5;         b[2]=4;
*/	
//do the first line and column because of the unknown parameters
	temp=fabs(A[0][0]);
	j=0;                   //the index of the max

	for (i=0;i<N;i++)
	{
		if(temp<fabs(A[i][0]))
		{
			temp=fabs(A[i][0]);
			j=i;
		}
	}
	if (j!=0)                  //exchange  two lines
	{
		for (i=0;i<N;i++)
		{
			temp=A[0][i];
			A[0][i]=A[j][i];
			A[j][i]=temp;
		}
		temp=b[0];
		b[0]=b[j];
		b[j]=temp;
	}

	//the first line and column action
	for (i=0;i<N;i++)
	{
		A[0][i]=A[0][i];
	}
	for (i=1;i<N;i++)
	{
		A[i][0]=A[i][0]/A[0][0];
	}


	for (ai=1,aj=1;aj<N;)
	{
		i=FindMaxLine(A,aj,N);
		if (ai!=i)
		{
			ExchangeLines(A,N,ai,i,b);
		}
		Angel(A,aj,N,b);
		ai++;
		aj++;		
	}
	
	for (i=N-1;i>=0;i--)
	{
		b[i]=SolveEquation(A,b,i,N);
	}

	cout<<"the answer is "<<endl;


	for (i=0;i<N;i++)
	{
		cout<<b[i]<<"    ";
	}
	cout<<endl;
	
/*	
	for (ai=0,aj=0;ai<N;)
	{
		cout<<A[ai][aj]<<"     ";
		aj++;
		if (aj==N)
		{
			cout<<b[ai]<<endl;
			ai++;
			aj=0;
		}
	}
*/



//release the memory
	delete []b;

	for(i=0;i<N;i++)
	{
		delete [] A[i];
	}

 	delete [] A;
		
	A = NULL;	
	
}


int FindMaxLine(float **A,int column,int Width)
{
	/*
	A coefficient matrix
	colume  the colume that need compared
	Height the Height   A
	the function return the line that 
	*/
	int maxline;
	float temp;
	int i;
	int j=0;
	float temp2;

//calculate sk    
	temp=A[column][column];
	for (j=0;j<=column-1;j++)
	{
		temp-=A[column][j]*A[j][column];
	}
	temp=fabs(temp);


	maxline=column;
	for ( i=column+1;i<Width;i++)
	{

		temp2=A[i][column];
		for (j=0;j<=column-1;j++)
		{
			temp2-=A[i][j]*A[j][column];
		}
		temp2=fabs(temp2);

		if (temp<temp2)
		{
			maxline=i;
			temp=temp2;
		}
	}
	return maxline;	
}


void ExchangeLines(float **A,int Width,int lineA,int lineB,float *b)
{
	float temp;
	for(int i=0;i<Width;i++)
	{
		temp=A[lineA][i];
		A[lineA][i]=A[lineB][i];
		A[lineB][i]=temp;
	}
	temp=b[lineA];
	b[lineA]=b[lineB];
	b[lineB]=temp;
}

void Angel(float **A,int column,int Width,float *b)
{
	//the startline and column is equal;
	float temp;
	int i;
	int j;
//calculate the U[column][j]

	for (j=column;j<Width;j++)
	{
		temp=A[column][j];
		for (i=0;i<=column-1;i++)
		{
			temp-=A[column][i]*A[i][j];
		}
		A[column][j]=temp;
	}
//calculate the b[column]
	temp=b[column];
	for (i=0;i<=column-1;i++)
	{
		temp-=A[column][i]*b[i];
	}
	b[column]=temp;


//calculate the l[j][column]
	for (j=column+1;j<Width;j++)
	{
		temp=A[j][column];
		for (i=0;i<=column-1;i++)
		{
			temp-=A[j][i]*A[i][column];
		}

		temp/=A[column][column];

		A[j][column]=temp;

	}
}


float SolveEquation(float **A,float *b,int xi,int Width)
{
/* solve the  xi*/
	int i;
	float temp;
	temp=0;
	for (i=xi+1;i<Width;i++)
	{
		float aa=A[xi][i]*b[i];
		temp+=A[xi][i]*b[i];
	}
	temp=(b[xi]-temp)/A[xi][xi];
	return temp;
}