cpp1.h

高斯消去法解线性方程组

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

int gauss(int n,double a[],double b[])
{
	int *js,flag,k,i,j,is,p,q;
	double max,t;

	js=(int *)malloc(n*sizeof(int));
	flag=1;//置非奇异标置
	for(k=0;k<=n-2;k++)
	{
		max=0.0;
		for(i=k;i<n;i++)
			for(j=k;j<n;j++)
			{
				t=fabs(a[i*n+j]);
				if(t>max)//记忆行列的交换信息
				{ max=t; is=i; js[k]=j; }
			}

			if(max+1.0==1.0)//如果最大元为0，则矩阵为奇异矩阵，置奇异标志
				flag=0;
			else
			{
				if(is!=k)//进行行交换
				{
					for(j=k;j<n;j++)//系数矩阵行交换
					{
						p=is*n+j;
						q=k*n+j;
						t=a[p];
						a[p]=a[q];
						a[q]=t;
					}
					//常数矩阵的行交换,千万不能遗漏!!
					t=b[k];
					b[k]=b[is];
					b[is]=t;
				}
				if(js[k]!=k)//进行列交换
					for(i=0;i<n;i++)
					{
						p=i*n+k;
						q=i*n+js[k];
						t=a[p];
						a[p]=a[q];
						a[q]=t;
					}
			}

			if(flag==0)
			{ free(js); printf("fail\n"); return(0);}

			max=a[k*n+k];
			for(j=k+1;j<n;j++)//第k步进行第k行的归一化,即从该行的第k+1个元素开始
			{
				p=k*n+j;
				a[p]=a[p]/max;
			}
			b[k]=b[k]/max;

			for(i=k+1;i<n;i++)//从第k+1行对k+1及以后的元素进行系数矩阵的消元
			{
				for(j=k+1;j<n;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];
			}
	}

	max=a[n*(n-1)+n-1];
	if(fabs(max)+1.0==1.0)
	{ free(js); printf("fail\n"); return(0);}

	b[n-1]=b[n-1]/max;
	for(i=n-2;i>=0;i--)
	{
		t=0.0;
		for(j=i+1;j<n;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);
}