列主元高斯消元法.txt

来自「拉格朗日插值」· 文本 代码 · 共 83 行

#include<stdio.h>
#include<math.h>
#define H 20
#define L 30
main()
{
  int n,i,j,k,p;
  float a[H][L],m[H][L],b[H],x[L];
  float max,temp1,temp2,E;
  printf("Input n:");
  scanf("%d",&n);
  printf("Input a:\n");
  for(i=0;i<=n-1;i++)
  {
    for(j=0;j<=n-1;j++)
	{
	  scanf("%f,",&a[i][j]);
	}
  }
  printf("Input b:\n");
  for(i=0;i<=n-1;i++)
  {
    scanf("%f,",&b[i]);
  }
  for(k=0;k<=n-2;k++)
  {
    max=(float)fabs(a[k][k]);
	for(i=k;i<=n-1;i++)
	{
	  if((float)fabs(a[i][k])>(float)fabs(a[k][k]))
	  {
	    p=i;
		max=a[p][k];
	  }
	}
    if(max==0)
	{
	  printf("qi yi!\n");
	}
    if(p!=k)
	{
	  temp2=b[p];
	  b[p]=b[k];
	  b[k]=temp2;
	  for(j=k;j<=n-1;j++)
	  {
	    temp1=a[p][j];
        a[p][j]=a[k][j];
	    a[k][j]=temp1;
	  }
	} 
	for(i=k+1;i<=n-1;i++)
	{
	  m[i][k]=a[i][k]/a[k][k];
	  b[i]=b[i]-b[k]*m[i][k];
	  for(j=k+1;j<=n-1;j++)
	  {
	    a[i][j]=a[i][j]-a[k][j]*m[i][k];
	  }
	}
  }
  if(a[n-1][n-1]==0)
  {
    printf("qi yi!\n");
  }
  else
  {
	x[n-1]=b[n-1]/a[n-1][n-1];
	for(k=n-2;k>=0;k--)
	{
	  E=0.0;
	  for(j=k+1;j<=n-1;j++)
	  {
	    E=E+a[k][j]*x[j];
	  }
      x[k]=(b[k]-E)/a[k][k];
	}
    for(i=0;i<=n-1;i++)
	{
      printf("x[%d]=%f\n",i+1,x[i]);
	}
  }
}