[数值算法]高斯消元法改进版---列主消元法

[数值算法]高斯消元法改进版---列主消元法

高斯消元在求解一些系数矩阵中含有极小数的情况下，会产生巨大的舍入误差，导致算法失效。一个简单而有效的改进方法是每次在进行将当前列中元素的消成0的运算时，选择当前列j对应的行(j to n)中绝对值较大的一个数作为主元行，这样，误差就会减小很多。具体详见算法描述:

一些必备的实用矩阵变换函数已写在高斯消元法中了，这里不再赘述了，参：

http://blog.csdn.net/EmilMatthew/archive/2005/08/15/454999.aspx

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

输入:

系数矩阵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

              得到确定当前i列中 matrixA_b [i][i]至matrixA_b [n][i]中绝对值最大的元素的行标，赋值给tmpColIndex（指标当前列中的(i行至n行中)最大值）

              If  tmpColIndex!=i then

                            Swap a中matrixA_b中第tmpColIndex行和第i行

End if

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 Col Main Algorithm

This version is implemented:

by EmilMatthew 05/8/16

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

*/

/*Aided function maxRowIndex,get the absolute max element’s row index in now colIndx*/

int maxRowIndex(Type** matrixArr,int colIndex,int len)

{

       int j;/*the iterator number*/

       int maxIndex;/*the tween min data*/

      

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

      

       maxIndex=colIndex;

      

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

              if(fabs(matrixArr[maxIndex][colIndex])<fabs(matrixArr[j][colIndex]))maxIndex=j;

      

       return maxIndex;

}

 

 

/*main gauss col main method*/

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

{

       Type** matrixA_b;

       Type** tmpMatrixArr,* tmpBList;

       Type m;/*the tween data*/

       int i,j;/*iterator num*/

       int tmpColIndex;/*the tween index data for main col*/

      

       /*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-col main version*/

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

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

                     {

                            tmpColIndex=maxRowIndex(matrixA_b,i,len);

                           

       if(tmpColIndex!=i)

       swapTwoRow(matrixA_b,i,tmpColIndex,len+1);/*still here,be attentionf:len+1!*/

                                  

                            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*/

              show2DArrFloat(tmpMatrixArr,len,len);

              showArrListFloat(tmpBList,0,len);

      

       upTriangleMethod(tmpMatrixArr,tmpBList,xAnsList,len);

 

 

       /*End of gauss method-col main version*/  

      

       /*memory free*/

      

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

              {

                     free(matrixA_b[i]);

                     free(tmpMatrixArr[i]);

              }

       free(matrixA_b);

       free(tmpMatrixArr);

       free(tmpBList);

}

 

 

/*Test program*/

/*Gauss column main 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*************/

              gaussMethodColMain(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;

}

 

 

测试结果：

test-1:

//input

3,

2,0,1,8;

0,4,6,12;

1,1,1,30;

//output

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

15.500000,37.500000,-23.000000;

the collapsed time in this algorithm implement is:0.000000

 

 

test-2:

//input

2,

0.0001,1,1;

1,1,3;

//output

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

2.000200,0.999800;

the collapsed time in this algorithm implement is:0.000000

 

 

test-3:

//input

3,

2,2,3,3;

4,7,7,1;

-2,4,5,-7;

//output

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

2.000000,-2.000000,1.000000;