gauss.cpp

本算法为全选主元高斯消去法

// Gauss.cpp : Defines the entry point for the console application.
//
#include <iostream.h>
#include "stdafx.h"
#include "stdlib.h"//增加
#include "math.h"//增加
int Gauss(double a[],double b[],long n)
{
	long *js,is,p,q;
	double d,t;
	js=(long *)malloc(n* sizeof(long));//malloc 内存块  sizeof 存储容量 
	long bChoice=1;
	for(long k=0;k<=n-2;k++)
	{
		d=0.0;
		for(long i=k;i<=n-1;i++)
			for(long j=k;j<=n-1;j++)
			{
				t=fabs(a[i*n+j]);
				if(t>d)
				{
					d=t;
					js[k]=j;
					is=i;
				}
			}
			if(d+1.0==1.0) 
				bChoice=0;
			else 
			{
				if(js[k] !=k)
					for(i=0;i<=n-1;i++)//交换位置
					{
						p=i*n+k;
						q=i*n+js[k];
						t=a[p];
						a[p]=a[q];
						a[q]=t;
					}
					if(is!=k)
					{
						for(long j=k;j<=n-1;j++)
						{
							p=k*n+j;
							q=is*n+j;
							t=a[p];
							a[p]=a[q];
							a[q]=t;
						}
						t=b[k];
						b[k]=b[is];
						b[is]=t;
					}
			}
			if(bChoice==0)
			{
				free(js);
				cout<<"fail"<<endl;
				//printf("fail\n");
				return 0;
			}
			d=a[k*n+k];
			for(long j=k+1;j<=n-1;j++)//归一化
			{
				p=k*n+j;
				a[p]=a[p]/d;
			}
			b[k]=b[k]/d;
			for(i=k+1;i<=n-1;i++)
			{
				for(j=k+1;j<=n-1;j++)
				{
					p=i*n+j;
					a[p]=a[p]-a[i*n+k]*a[k*n+j];//消去
				}
				b[i]=b[i]-a[i*n+k]*b[k];
			}
	}
	  d=a[(n-1)*n+n-1];
	  if(fabs(d)+1.0==1.0)
	  {
		  free(js);
		  cout<<"fail"<<endl;
		  //printf("fail\n");
		  return 0;
	  }
	  b[n-1]=b[n-1]/d;//回代过程
	  for(long i=n-2;i>=0;i--)
	  {
		  t=0.0;
		  for(long j=i+1;j<=n-1;j++)
			  t=t+a[i*n+j]*b[j];
		  b[i]=b[i]-t;
	  }
	  js[n-1]=n-1;
	  for(k=n-1;k>=0;k--)
	  if(js[k] !=k)
	  {
			  t=b[k];
			  b[k]=b[js[k]];
			  b[js[k]]=t;
	  }
		  free(js);
		  return 1;
}
int main(int argc, char* argv[])
{
	static double a[16]={0.2368,0.2471,0.2568,1.2671,0.1968,0.2071,1.2168,0.2271,0.1581,1.1675,0.1871,1.1161,0.1254,0.1397,0.1490};
	static double b[4]={1.8471,1.7471,1.6471,1.5471};
	if(Gauss(a,b,4) !=0)
		for(long i=0;i<=3;i++)
	       cout<<"x("<<i<<")="<<b[i]<<endl;
		//printf("x(%2d)=%14e\n",i,b[i]);
	return 0;
}