ex11.c

采用全选主元的高斯消去法解线性方程组

  #include "stdlib.h"
  #include "math.h"
  #include "stdio.h"

  //file <RGAUSS.C>全选主元高斯消去法

  int rgauss(n,a,b)
					//函数返回整型标志值。若返回0，表示原方程组的系数矩阵奇异；不为0，则表示正常返回
  int n;			//方程组的阶数
  double a[],b[];	//a[]存放方程组的系数矩阵，返回时被破坏；
					//b[]存放方程组右端常数向量，返回方程组的解向量；
  { 
    int *js,l,k,i,j,is,p,q;
    double d,t;
    js=malloc(n*sizeof(int));	//开辟用于记忆列交换信息的动态空间
    l=1;						//置非奇异标志
    for (k=0;k<=n-2;k++)		
      { d=0.0;
        for (i=k;i<=n-1;i++)	//全选主元
          for (j=k;j<=n-1;j++)
            { t=fabs(a[i*n+j]);
              if (t>d) { d=t; js[k]=j; is=i;}	//记忆行、列交换信息
            }
        if (d+1.0==1.0) l=0;	//置奇异标志
        else
          { if (js[k]!=k)
              for (i=0;i<=n-1;i++)	//列交换
                { p=i*n+k; q=i*n+js[k];
                  t=a[p]; a[p]=a[q]; a[q]=t;
                }
            if (is!=k)
              { for (j=k;j<=n-1;j++)	//行交换
                  { p=k*n+j; q=is*n+j;
                    t=a[p]; a[p]=a[q]; a[q]=t;
                  }
                t=b[k]; b[k]=b[is]; b[is]=t;
              }
          }
        if (l==0)	//奇异返回
          { free(js); printf("fail\n");
            return(0);
          }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)	  //归一化
          { p=k*n+j; a[p]=a[p]/d;}
        b[k]=b[k]/d;
        for (i=k+1;i<=n-1;i++)	  //消元
          { for (j=k+1;j<=n-1;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];
          }
      }
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)	//奇异返回
      { free(js); printf("fail\n");
        return(0);
      }
    b[n-1]=b[n-1]/d;	//计算解向量的最后一个分量
    for (i=n-2;i>=0;i--)	//回代
      { t=0.0;
        for (j=i+1;j<=n-1;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);	//返回正常标志
  }


  main()
  { 

	int i;
    static double a[4][4]=
           { {1,3,2,13},
             {7,2,1,-2},
             {9,15,3,-2},
             {-2,-2,11,5} };

     static double c[4][4]=
           { {1,3,2,13},
             {7,2,1,-2},
             {9,15,3,-2},
             {-2,-2,11,5} };
 
	static double b[4]={9,6,11,-2};
	static double d[4]={0,4,7,-1};

	if (rgauss(4,a,b)!=0)
		for (i=0;i<=3;i++)
			printf("x(%d1)=%e\n",i+1,b[i]);
	printf("\n");

	if (rgauss(4,c,d)!=0)
		for (i=0;i<=3;i++)
			printf("x(%d2)=%e\n",i+1,d[i]);
  }