[数值算法]求解线性方程组的高斯消元法

[数值算法]求解线性方程组的高斯消元法

高斯消元法在理论上还是很好理解的,但是由于在矩阵规模变大时,算法的可靠性极差,因此,它也是一个理论价值大于实用价值的算法,但同时也是后面求解行列式算法的基础.

这是一个高度模块化的算法,在实现高斯消元法之前,让我们做一些准备工作.

首先是一个求解上三角阵的回代算法,如下,主要思想就是通过前面求得的xn来计算下面的xn-1,还是很好理解的.如下:

void upTriangleMethod(Type** matrixArr,Type* bList,Type* xAnsList,int len)

{

       int i=0,j=0;

       double e=0.0001;/*precise controller*/

       Type tmpSum;

      

       assertF(matrixArr!=NULL,"in upAngleMethod,matrixArr is NULL/n");

       assertF(bList!=NULL,"in upAngleMethod,bList is NULL/n");

       assertF(xAnsList!=NULL,"in upAngleMethod,xAnsList is NULL/n");

      

for(i=0;i<len;i++)

assertF(fabs(matrixArr[i][i])>=e,"in upAngleMethod ,matrixArr[i][i] has value closed to 0/n");

      

/*Core part of the up triangle method*/

       xAnsList[len-1]=bList[len-1]/matrixArr[len-1][len-1];/*Get the last answer*/

       /*data initialiate*/

       i=len-2;

      

       while(i>=0)

       {

              tmpSum=0;

              for(j=i+1;j<len;j++)

                     tmpSum+=matrixArr[i][j]*xAnsList[j];

             

              xAnsList[i]=(bList[i]-tmpSum)/matrixArr[i][i];

              i--;

       }

       /*End of up triangle method*/

}

接下来要做的,就是准备好计算行列式所要用的三种基本操作:交换两行,a行数据自乘一个数据,a行乘以一个数据后加至B行.

代码如下:

/*Ulity matrix transformation functions*/

void swapTwoRow(Type** matrixArr,int row1Index,int row2Index,int colNums)

{

       int j;/*the iterator number*/

      

       /*assertion*/

       assertF(matrixArr!=NULL,"in swapTwoRow,*matrixArr is NULL/n");

      

       /*main swap process*/

       for(j=0;j<colNums;j++)

              swap(&matrixArr[row1Index][j],&matrixArr[row2Index][j]);

}

 

void valueByToRow(Type** matrixArr,Type valueBy,int rowIndex,int colNums)

{

       int j;/*the iterator number*/

      

       /*assertion*/

       assertF(matrixArr!=NULL,"in  valueByToRow,*matrixArr is NULL/n");

      

       for(j=0;j<colNums;j++)

              matrixArr[rowIndex][j]*=valueBy;

      

}

 

void valueByRowAPlusToRowB(Type** matrixArr,Type valueBy,int rowIndexFrom,

int rowIndexTo,int colNums)

{

       int j;/*the iterator number*/

       assertF(matrixArr!=NULL,"in  valueByRowAPlusToRowB,*matrixArr is NULL/n");

      

       for(j=0;j<colNums;j++)

              matrixArr[rowIndexTo][j]+=matrixArr[rowIndexFrom][j]*valueBy;

      

}

有了以上这些准备,再实现我们的高斯消元法就简单多了,再实现之前,还是先用算法语言描述一下,其实就是把我们手工计算时的操作映射到计算机上的一个过程.

输入:

系数矩阵a[n][n]以及右端列向量b[n]

主过程:

将a[n][n]和b[n]拼成增广矩阵a|b  这里用matrixA_b表示

For i<-1 to n-1

       For j<-i+1 to n

              m=matrixA_b [j][i]/ matrixA_b [i][i]

              matrixA_b第j行中的每个元素-m*matrixA_b第i行中的每个对应元素

(相当于上面的函数valueByRowAPlusToRowB)

              Next j

       Next i

       至些, matrixA_b [n][n]已化成上三角矩阵.

       再对它施加一次上三角矩阵的回代算法,即得结果.

       输出:

       线性方程组的解: x1-xn.

       下面给出我写的高斯消元法在C下的实现,为了防止对原数据的更改,这里在函数内部用了局部的指针变量,使得代码看上去长了些.

/*Guess Algorithm

This version is implemented:

by EmilMatthew 05/8/15

You can use these code freely , but I don’t assure them could totally fit your needs.

*/

void gaussMethod(Type** matrixArr,Type* bList,Type* xAnsList,int len)

{

       Type** matrixA_b;

       Type** tmpMatrixArr,* tmpBList;

       Type m;/*the tween data*/

       int i,j;/*iterator num*/

      

       /*input assertion*/

       assertF(matrixArr!=NULL,"in gaussMethod,matrixArr is NULL/n");

       assertF(bList!=NULL,"in gaussMethod,bList is NULL/n");

       assertF(xAnsList!=NULL,"in gaussMethod,xAnsList is NULL/n");    

      

       /*memory apply*/

       matrixA_b=(Type**)malloc(sizeof(Type*)*len);

              for(i=0;i<len;i++)/*the col's size should be row's size+1*/

                     matrixA_b[i]=(Type*)malloc(sizeof(Type)*(len+1));

      

       tmpMatrixArr=(Type**)malloc(sizeof(Type*)*len);

              for(i=0;i<len;i++)

                     tmpMatrixArr[i]=(Type*)malloc(sizeof(Type)*len);

      

       tmpBList=(Type*)malloc(sizeof(Type)*len);

      

       assertF(matrixA_b!=NULL,"in gaussMethod,matrixA_b is NULL/n");

       assertF(tmpMatrixArr!=NULL,"in gaussMethod,tmpMatrixArr is NULL/n");

       assertF(tmpBList!=NULL,"in gaussMethod,tmpBList is NULL/n");

      

             

       /*matrixA_b make*/

       for(i=0;i<len;i++)

       {

              for(j=0;j<len;j++)

                            matrixA_b[i][j]=matrixArr[i][j];

              matrixA_b[i][j]=bList[i];

       }

      

       /*Core part of gauss matrix method*/

              for(i=0;i<len-1;i++)

                     for(j=i+1;j<len;j++)

                     {

                            m=matrixA_b[j][i]/matrixA_b[i][i];

                            valueByRowAPlusToRowB(matrixA_b,-m,i,j,len+1);                

                     }

              /*The matrixA_b has become a triangle matrix*/

              /*split the matrixA_b back to matrixA and bList*/

             

              for(i=0;i<len;i++)

              {

                     for(j=0;j<len;j++)

                            tmpMatrixArr[i][j]=matrixA_b[i][j];

                     tmpBList[i]=matrixA_b[i][j];

              }

      

              /*implement the up triangle algorithm*/

              /*debug code*/

show2DArrFloat(tmpMatrixArr,len,len);

              showArrListFloat(tmpBList,0,len);

              /*end of debug code*/

       upTriangleMethod(tmpMatrixArr,tmpBList,xAnsList,len);

 

       /*End of gauss method*/      

      

       /*memory free*/

      

       for(i=0;i<len;i++)

              {

                     free(matrixA_b[i]);

                     free(tmpMatrixArr[i]);

              }

       free(matrixA_b);

       free(tmpMatrixArr);

       free(tmpBList);

}

 

/*Test Program*/

/*upTiangleMethod  Algorithm test program*/

#include "Global.h"

#include "Ulti.h"

#include "MyAssert.h"

#include "Matrix.h"

#include <time.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

 

 

char *inFileName="inputData.txt";

/*

       input data specification

       len,

       a00,a01,...,a0n-1,b0;

       .....

       an-10,an-11,...,an-1n-1,bn-1;

*/

 

 

char *outFileName="outputData.txt";

#define DEBUG 1

 

void main(int argc,char* argv[])

{

       FILE *inputFile;/*input file*/

       FILE *outputFile;/*output file*/

 

       double startTime,endTime,tweenTime;/*time callopsed info*/

      

       /*The read in data*/

       int len;

       Type** matrixArr;

       Type* bList,* xAnsList;

 

      

       int i,j;/*iterator index*/

      

       /*input file open*/

       if(argc>1)strcpy(inFileName,argv[1]);

       assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");

       printf("input file open success/n");

      

       /*outpout file open*/

       if(argc>2)strcpy(outFileName,argv[2]);

       assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");

       printf("output file open success/n");    

      

       fscanf(inputFile,"%d,",&len);

       printf("len is:%d/n",len);

      

       /*Memory apply*/

       matrixArr=(Type**)malloc(sizeof(Type*)*len);

              for(i=0;i<len;i++)

                     matrixArr[i]=(Type*)malloc(sizeof(Type)*len);

      

       bList=(Type*)malloc(sizeof(Type)*len);

       xAnsList=(Type*)malloc(sizeof(Type)*len);

             

       /*Read info data*/

       for(i=0;i<len;i++)

       {

              for(j=0;j<len;j++)

                     fscanf(inputFile,"%f,",&matrixArr[i][j]);

              fscanf(inputFile,"%f;",&bList[i]);

       }

      

      

       /*Check the input data*/

       showArrListFloat(bList,0,len);

       show2DArrFloat(matrixArr,len,len);

      

#if  DEBUG

       printf("/n*******start of test program******/n");

       printf("now is runnig,please wait.../n");

       startTime=(double)clock()/(double)CLOCKS_PER_SEC;

       /******************Core program code*************/

              gaussMethod(matrixArr,bList,xAnsList,len);

              printf("after the upTriangleMethod:the ans x rows is:/n");

              fprintf(outputFile,"after the upTriangleMethod:the ans x rows is:(from x0 to xn-1)/r/n");

              showArrListFloat(xAnsList,0,len);

              outputListArrFloat(xAnsList,0,len,outputFile);

       /******************End of Core program**********/

       endTime=(double)clock()/(double)CLOCKS_PER_SEC;

       tweenTime=endTime-startTime;/*Get the time collapsed*/

       /*Time collapsed output*/

       printf("the collapsed time in this algorithm implement is:%f/n",tweenTime);

       fprintf(outputFile,"the collapsed time in this algorithm implement is:%f/r/n",tweenTime); 

       printf("/n*******end of test program******/n");

#endif

       for(i=0;i<len;i++)

              free(matrixArr[i]);

       free(matrixArr);

      

       free(xAnsList);

       free(bList);

      

       printf("program end successfully,/n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer./n");

       getchar();/*Screen Delay Control*/

       return;

}

/*测试结果*/

       测试一:

输入数据:

       3,//矩阵长度

2,2,3,3;//a11-a13,b1

4,7,7,1;

-2,4,5,-7;//a31-a33,b3.

 

输出:

after the gauss algorithm:the ans x rows is:(from x0 to xn-1)

2.000000,-2.000000,1.000000;

 

测试2:

输入数据:

3,

2,0,1,8;

0,4,6,12;

1,1,1,30;

输出:

after the gauss algorithm:the ans x rows is:(from x0 to xn-1)

15.500000,37.500000,-23.000000;