gauss.cpp

常用数值算法的C++源码

#include "stdio.h"
#include "math.h"
double *x;

void Gauss(double *a,int n)
{
	int    m,N=n+1;
	double max,temp,sum;
	for(int k=0;k<n-1;k++)
	{
		max=fabs(a[k*N+k]);
		m=k;
		//寻找列主元
		for(int i=k+1;i<n;i++)
		{
			temp=fabs(a[i*N+k]);
			if(temp>max)
			{
				max=temp;
				m=i;
			}
		}
		//换行
		if(m!=k)
		{
			for(i=k;i<N;i++)
			{
				temp=a[m*N+i];
				a[m*N+i]=a[k*N+i];
				a[k*N+i]=temp;
			}
		}
		//消元
		for(i=k+1;i<n;i++)
		{
			if(a[i*N+k]!=0)
			{
				for(int j=k;j<N;j++)
					x[j]=a[i*N+j]-a[k*N+j]*a[i*N+k]/a[k*N+k];
				for(j=k;j<N;j++)
					a[i*N+j]=x[j];
			}
		}
	}
    //回代
	x[n-1]=a[(n-1)*N+n]/a[(n-1)*N+(n-1)];
	for(int i=n-2;i>=0;i--)
	{
		sum=0;
		for(int j=i+1;j<n;j++)
			sum=sum+a[i*N+j]*x[j];
		x[i]=(a[i*N+n]-sum)/a[i*N+i];
	}
}

void main()
{
	double r[90]={31,-13,0,0,0,-10,0,0,0,  -15,
		          -13,35,-9,0,-11,0,0,0,0,  27,
				  0,-9,31,-10,0,0,0,0,0,   -23,
				  0,0,-10,79,-30,0,0,0,-9,   0,
				  0,0,0,-30,57,-7,0,-5,0,  -20,
				  0,0,0,0,-7,47,-30,0,0,    12,
				  0,0,0,0,0,-30,41,0,0,     -7,
				  0,0,0,0,-5,0,0,27,-2,      7,
				  0,0,0,-9,0,0,0,-2,29,     10};
	double i[10];
	int n=9;
	x=i;
	Gauss(r,n);
	printf("解向量为：\n");
	for(int j=0;j<9;j++)
		printf("%8.5f\n",x[j]);
}