高斯消元法

#include <iostream.h>
#include <math.h>
const int n=4;

/*下面是求列中最大值函数，考察范围是：a[k][k]、a[k+1][k]、……、a[n-1][k],返回结果是主元的行号，为选择主元作准备。*/
int Max(double a[n][n+1],int k)
{
 int t=k;
 double max=fabs(a[k][k]);
 for(int i=k+1;i<=n-1;i++)
 {
  if(fabs(a[i][k])>max)
  {
   t=i;
   max=fabs(a[i][k]);
  }
 }
 return t;
}

/*下面是一个交换第k行和第t行方程的函数，选主元的后续工作*/
void change(double a[n][n+1],int k,int t)
{
 double x;
 for(int j=k;j<=n;j++)
 {
  x=a[k][j];
  a[k][j]=a[t][j];
  a[t][j]=x;
 }
}

/*下面是高斯消去法的核心函数，其中调用求主元、交换方程函数*/
void gauss(double a[n][n+1])
{
 for(int k=0;k<n-1;k++)
 {
  int t=Max(a,k);
  if(t!=k)
   change(a,k,t);
  for(int i=k+1;i<n;i++)
  {
   a[i][k]=a[i][k]/a[k][k];
   for(int j=k+1;j<=n;j++)
   {
    a[i][j]=a[i][j]-a[i][k]*a[k][j];
   }
  }
 }
}

/*下面是解上三角方程组函数*/
void LYB(double a[n][n+1])
{
  for(int i=n-1;i>=0;i--)
  {
    for(int j=i+1;j<=n-1;j++)
      a[i][n]-=a[i][j]*a[j][n];
    a[i][n]/=a[i][i];
  }
}

/*主函数*/
void main()
{
 double A[n][n+1]={{1,-1,1,-4,2},{5,-4,3,12,4},
{2,1,1,11,3},{2,-1,7,-1,0}};
 gauss(A);
 LYB(A);
 for(int i=0;i<n;i++)
  cout<<"x["<<i<<"]="<<A[i][n]<<endl;
 cout<<endl;
}