ordgauss.c

来自「利用高斯法求线性方程组的解」· C语言 代码 · 共 56 行

/*#include "create1.c"*/
#include "head.c"

int ordGauss(double in_aa[C_M][C_N],double *in_bb,double *out)
{
	int n;
	int k,i,j;
	double temp;
	double m;
	double in_a[C_M][C_N];
	double in_b[Matrix_N];

	for (i=0;i<C_M ;i++ )
	{
		for (j=0;j<C_N ;j++ )
		{
			in_a[i][j]=in_aa[i][j];
		}
	}
	for (i=0;i<Matrix_N ;i++ )
	{
		in_b[i]=in_bb[i];
	}

	for (k=0;k<Matrix_N-1 ;k++ )
	{	
		if (in_a[convert_i(k,k)][k]==0.0)
		{
			printf("a[%d][%d]==%f,arithmetic failed!\n",k,k,in_a[k][k]);
			return(0);
		}
		else
		{
			for(i=k+1;i<=k+2;i++)
			{
				if (i<Matrix_N) 
				{
					m=in_a[convert_i(i,k)][k]/in_a[convert_i(k,k)][k];
					in_a[convert_i(i,(k+1))][k+1]-=m*in_a[convert_i(k,(k+1))][k+1];
					in_a[convert_i(i,(k+2))][k+2]-=m*in_a[convert_i(k,(k+2))][k+2];
					in_b[i]-=m*in_b[k];
				}
			}
		}
	}
	out[Matrix_N-1]=in_b[Matrix_N-1]/in_a[convert_i((Matrix_N-1),(Matrix_N-1))][Matrix_N-1];
	for (k=Matrix_N-2;k>=0 ;k--)
	{
		temp=0;
		for (j=k+1;j<=k+2 ;j++ )
			if(j<Matrix_N)
				temp+=in_a[convert_i(k,j)][j]*out[j];
		out[k]=(in_b[k]-temp)/in_a[convert_i(k,k)][k];
	}
	return(1);
}