gauss.cpp

高斯消元

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define e 0.000000001

double liezhuyuan(double A[3509][3509],double b[3509],double x[3509])
{
	int k,i,j,i_1,n=3509;
	double max,sum,m;


	/* 消元过程 */
	for(k=0;k<n;k++)
	{
		/* 选主元*/
		max=fabs(A[k][k]);       
		i_1=k;
		for(i=k;i<n;i++)
			if( fabs(A[i][k]) > max ) 
			{ 
				max=fabs(A[i][k]); 
				i_1=i; 
			}
			if( i_1 != k )
			{
				for(j=k;j<n;j++)
				{ 
					max=A[k][j];
					A[k][j]=A[i_1][j];
					A[i_1][j]=max;
				}
				max=b[k];
				b[k]=b[i_1];
				b[i_1]=max;
			}
			if( fabs(A[k][k]) < e ) 
				return(1);
			for(i=k+1;i<n;i++)
			{
				m=A[i][k]/A[k][k];
				for(j=k;j<n;j++) 
					A[i][j]=A[i][j]-m*A[k][j];
				b[i]=b[i]-m*b[k];  
			}

	}
	/* 回带过程 */
	x[n-1]=b[n-1]/A[n-1][n-1];
	for(i=n-2;i>=0;i--)
	{   
		for(sum=0.0,j=i+1;j<n;j++) 
			sum=sum+A[i][j]*x[j];
		x[i]=b[i]-sum;
		x[i]=x[i]/A[i][i];
	}
/*	for(i=0;i<3509;i++)
	{
		for(j=0;j<3509;j++)
		{
			printf("%f ",A[i][j]);
		}
		printf("\n");
	}*/
	
}