数值作业：Guass全选主元消去法之C语言代码

比眉伴天荒 2022-06-17 06:39 201阅读 0赞

全选主元的基本思想:当变换到第K步时，从系数矩阵的右下角(n-k+1)阶子阵中选取绝对值最大的元素，然后通过行变换与列变换将它交换到主元素a\[k\]\[k\]的位置上．行变换是不影响最后求解的结果的，但是列变换，相应的两个未知数的次序也被交换了．因此，全选主元的高斯消去法，必须记住进行的一切列交换，一边最后的还原.下面是代码：

/********************************************
        > File Name: Allselect.c
        > Author:chendiyang
        > School:WUST_CST_1501班
        > Myblog:www.chendsir.com
        > Mail:1441353519@qq.com 
        > Created Time: 2017年04月19日 星期三 13时14分40秒
     *****************************************************************
    #include <stdlib.h>
    #include <math.h>
    #include <stdio.h>
    #define MAX 20
    
    
    
    int main()
    {
        int n;       //矩阵的维数　
        static double A[MAX][MAX],b[MAX],x[MAX];//有A*x=b,分别代表不同系数矩阵
        int *js,l,k,i,j,is,p,q;//用于记忆行列交换信息的需要的临时变量，其中指针js用于开辟动态空间
        double d,t;
        printf("\n 输入方程组的维数:");//输入AX=b的维数
        scanf("%d",&n);
        if(n>MAX)
        {
            return 1;
         } 
         if(n<=0)
         {
            return 1;
         }
         //输入A矩阵的值
         printf("\n请输入A矩阵的值:"); 
         for(i=0;i<n;i++)
           for(j=0;j<n;j++) 
             scanf("%lf",&A[i][j]);
               //输入b矩阵 
          printf("\n请输入系数B矩阵的值:");      
          for(i=0;i<n;i++)
          scanf("%lf",&b[i]);
          double *a=(double *)A;
        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++)  //当然是绝对值最大一个，这里用到了数学库绝对值函数fabs()
            {
                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("原方程的系数矩阵奇异，无法进行高斯消去法");
                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("原方程的系数矩阵奇异，无法进行高斯消去法");
            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);          //释放动态空间
    
    
    
        if(l!=0)
        for(i=0;i<=3;i++)
        printf("x(%d)=%e\n",i,b[i]);
    
    }

运行结果：  
　　　　![这里写图片描述][SouthEast]　　　　　　　  
　　　　很好，我们看到结果显示：原方程的系数矩阵奇异，无法进行高斯消去，喔，这个就表明，这个TMD运行有问题，啊啊啊，惨，敲这么久你给我看这个，实在是难以静下心去调试，这个有问题的代码就先放这吧，至少原理是正确的．衰．．

[SouthEast]: /images/20220617/9e7344b899d44964b90e8e7ffb796785.png