gaussian.c

高斯消去法的c语言实现

//Using Gaussian elimination method and Gaussian elimination method with row scaled method to solve the following tri-diagonal system
// 
//for n=10 and 100. 
//
//解决方案：由于使用高斯消去发的前提是，在运算过程中A矩阵的主对角线元素必须是非零的，但该矩阵不满足此条件无法进行高斯消去，在实验中根据高斯法编写代码进行运算过程中系统报错，因为计算中出现初以零运算。以下为高斯行尺度法的C语言实现。
 
#include <stdio.h>
#include <math.h>
#define N 10

/*swap function is used to exchange value between x and y.*/

swap(double *p1,double *p2)
  {
   double temp;
   temp=*p1;
   *p1=*p2;
   *p2=temp;
  }

print_column_N(double array[N][N])
  {
   int i,j;
   for(i=0;i<N;i++)
     {
      for(j=0;j<N;j++)
        printf("%5f ",array[i][j]);
      printf("\n");
     }
  getchar();
  }

print_column_1(double array[N])
  {
   int i;
   for(i=0;i<N;i++)
     {
      printf("%5f ",array[i]);
      printf("\n");
     }
  getchar();
  }

main()
{
 int i,j,m,max;
 double *pointer_1,*pointer_2;
 double k;
 double s[N];
 double a_s[N];
 double l[N][N];
 double u[N][N];
 double a[N][N];
 double b[N];
 double z[N];
 double x[N];

/*initial a[N][N],b[N] ,l[N][N],u[N][N] and row scaled s[N].*/

     for(i=0;i<N;i++)
     {
      for(j=0;j<N;j++)
      {
       a[i][j]=0;
      }
      b[i]=0;
     }

     for(i=0;i<N-1;i++)
      {
	a[i][i]=3;
	a[i][i+1]=1;
	a[i+1][i]=9;
      }
      a[N-1][N-1]=3;
      b[0]=4;
      b[N-1]=12;
    for(i=1;i<N-1;i++)
      {
       b[i]=13;
      }

 for(i=0;i<N;i++)
 {
   for(j=0;j<N;j++)
   {
    u[i][j]=0;
    if(i==j)
      l[i][j]=1;
    else
      l[i][j]=0;
   }
 }


for(i=0;i<N;i++)
{
 s[i]=fabs(a[i][0]);
 for(j=1;j<N;j++)
 {
  if(fabs(a[i][j])>fabs(s[i]))
    s[i]=fabs(a[i][j]);
 }
}

/*examine a[N][N],b[N].*/

print_column_N(a);

print_column_1(b);

print_column_1(s);

/*gaussian_scaled arithmetic core.*/

 for(j=0;j<N-1;j++)
 {
  for(m=j;m<N;m++)
  {
   a_s[m]=fabs(a[m][j])/s[m];
  }

  max=j;
  for(m=j+1;m<N;m++)
  {   if(a_s[max]<a_s[m])
   max=m;
  }

 if(max!=m)
 {
  for(m=0;m<N;m++)
  {
   pointer_1=&a[j][m];
   pointer_2=&a[max][m];
   swap(pointer_1,pointer_2);
  }
  pointer_1=&s[j];
  pointer_2=&s[max];
  swap(pointer_1,pointer_2);
  pointer_1=&b[j];
  pointer_2=&b[max];
  swap(pointer_1,pointer_2);
 }

  for(i=j+1;i<N;i++)
  {
   k=a[i][j]/a[j][j];
   a[i][j]=k;

   for(m=j+1;m<N;m++)
   {
    a[i][m]=a[i][m]-k*a[j][m];
   }
  }
 }

/*generate l[N][N] and u[N][N].*/

 for(i=0;i<N;i++)
 {
  for(j=0;j<N;j++)
  {
   if(i>j)
     l[i][j]=a[i][j];
   else
     u[i][j]=a[i][j];
   }
 }

/*print l[N][N] and u[N][N].*/

 printf("\n");
 printf("L:\n");
 print_column_N(l);

 printf("\n");
 printf("U:\n");
 print_column_N(u);

/*calculate z[N] and print.*/

 for(i=0;i<N;i++)
 {
  if(i==0)
  z[i]=b[i];
  else
  {
   z[i]=b[i];
   for(j=0;j<i;j++)
   {
    z[i]=z[i]-l[i][j]*z[j];
   }
  }
 }

 printf("\n");
 printf("Z:\n");
 print_column_1(z);

/*calculate x[N] and print.*/


 for(i=N-1;i>=0;i--)
 {
  if(i==N-1)
  x[i]=z[i]/u[i][i];
  else
  {
   x[i]=z[i];
   for(j=N-1;j>i;j--)
   {
    x[i]=x[i]-x[j]*u[i][j];
   }
   x[i]=x[i]/u[i][i];
  }
 }

 printf("\n");
 printf("X:\n");
 print_column_1(x);

}