高斯消元法简单实现

这里只是简单对原理的实现，以一个简单的已输入实例为例

#include<stdio.h>

#define row 3
#define colomn 3

void swap(double *a,double *b);

int main()
{
    int i,j,m,k,ii,jj;
    double a[row][colomn]={{2,3,1},{4,2,3},{7,1,-1}};
    double b[colomn]={4,17,1};
    int index;
    char str[]="matrix.txt";

    for(ii=0;ii<row;ii++)
    {
        for(jj=0;jj<colomn;jj++)
        {
            printf("%5.2f\t",a[ii][jj]);
        }
        printf("%5.2f\n",b[ii]);
    }
    printf("\n*******************************\n");
    for(i=0;i<row;i++)
    {    
        index=i;
        while(a[index][i]==0)
        {
            index++;
        }
        if(index!=i)
        {
            for(k=0;k<colomn;k++)
            {
                swap(*(a+i)+k,*(a+index)+k);
            }
            swap(b+index,b+i);
            //swap row i and index
        }
        for(m=i;m<row;m++)
        {
            if(m==i)
            {
                b[m]=b[m]/a[i][i];
            }
            else
            {
                /////////***********//////
                b[m]=b[m]-b[i]*a[m][i]/a[i][i];
            }
            for(j=colomn-1;j>=0;j--)
            {
                if(m==i)
                {
                    a[m][j]=a[m][j]/a[i][i];
                }
                else
                {
                    /////////***************//
                    a[m][j]=a[m][j]-a[i][j]*a[m][i];
                }
            }
        }
        //print every step
        for(ii=0;ii<row;ii++)
        {
            for(jj=0;jj<colomn;jj++)
            {
                printf("%5.2f\t",a[ii][jj]);
            }
            printf("%5.2f\n",b[ii]);
        }
        printf("\n*******************************\n");
    }
    index=0;
    while(a[index][index]!=0&&index<row)
        index++;
    for(i=index-1;i>0;i--)
    {
        for(j=i-1;j>=0;j--)
        {
            //can't change turn
            b[j]=b[j]-a[j][i]*b[i];
            a[j][i]=a[j][i]-a[i][i]*a[j][i];
        }
        //print a 
        for(ii=0;ii<row;ii++)
        {
            for(jj=0;jj<colomn;jj++)
            {
                printf("%5.2f\t",a[ii][jj]);
            }
            printf("%5.2f\n",b[ii]);
        }
        printf("\n*******************************\n");
    }
    printf("the answer is :\n");
    for(i=0;i<colomn;i++)
    {
        printf("\t%6.2f\n",b[i]);
    }
    return 0;
}

void swap(double *a,double *b)
{
    double temp;
    temp=*a;
    *a=*b;
    *b=temp;
}

以下使用了c文件的读取来完成，避免了数据输入的繁琐，程序直接从txt文件中读取。

#include<stdio.h>
#include<malloc.h>

void swap(double *a,double *b);

int main()
{
    int i,j,m,k,ii,jj;
    int row,colomn;
    double *a,*b;
    int index;
    char str[]="matrix.txt";

    FILE *fr;
    if(fr=fopen(str,"r"))
    {
        fscanf(fr,"%d %d",&row,&colomn);
        a=(double *)malloc(sizeof(double)*row*colomn);
        for(i=0;i<row;i++)
        {
            for(j=0;j<colomn;j++)
                fscanf(fr,"%lf",a+i*row+j);
        }
        b=(double *)malloc(sizeof(double)*row);
        for(i=0;i<row;i++)
        {
            fscanf(fr,"%lf",b+i);
        }
        fclose(fr);
    }
    else
    {
        printf("read error !");
    }
    //print a and b
    for(ii=0;ii<row;ii++)
    {
        for(jj=0;jj<colomn;jj++)
        {
            printf("%5.2f\t",*(a+ii*row+jj));
        }
        printf("%5.2f\n",b[ii]);
    }
    printf("\n*******************************\n");
    for(i=0;i<row;i++)
    {    
        index=i;
        while(*(a+index*row+i)==0)
        {
            index++;
        }
        if(index!=i)
        {
            for(k=0;k<colomn;k++)
            {
                swap(a+i*row+k,a+index*row+k);
            }
            swap(b+index,b+i);
            //swap row i and index
        }
        for(m=i;m<row;m++)
        {
            if(m==i)
            {
                b[m]=b[m]/a[i*row+i];
            }
            else
            {
                /////////***********//////
                b[m]=b[m]-b[i]*a[m*row+i]/a[i*row+i];
            }
            for(j=colomn-1;j>=0;j--)
            {
                if(m==i)
                {
                    a[m*row+j]=a[m*row+j]/a[i*row+i];
                }
                else
                {
                    /////////***************//
                    a[m*row+j]=a[m*row+j]-a[i*row+j]*a[m*row+i];
                }
            }
        }
    }
    index=0;
    while(a[index*row+index]!=0&&index<row)
        index++;
    for(i=index-1;i>0;i--)
    {
        for(j=i-1;j>=0;j--)
        {
            //can't change turn
            b[j]=b[j]-a[j*row+i]*b[i];
            a[j*row+i]=a[j*row+i]-a[i*row+i]*a[j*row+i];
        }
        //printf a 
        for(ii=0;ii<row;ii++)
        {
            for(jj=0;jj<colomn;jj++)
            {
                printf("%5.2f\t",*(a+ii*row+jj));
            }
            printf("%5.2f\n",b[ii]);
        }
        printf("\n*******************************\n");
    }
    
    printf("the answer is :\n");
    for(i=0;i<colomn;i++)
    {
        printf("\t%6.2f\n",b[i]);
    }
    free(a);
    free(b);
    return 0;
}

void swap(double *a,double *b)
{
    double temp;
    temp=*a;
    *a=*b;
    *b=temp;
}