[数值算法]求解线性方程组的高斯消元法
高斯消元法在理论上还是很好理解的,但是由于在矩阵规模变大时,算法的可靠性极差,因此,它也是一个理论价值大于实用价值的算法,但同时也是后面求解行列式算法的基础.
这是一个高度模块化的算法,在实现高斯消元法之前,让我们做一些准备工作.
首先是一个求解上三角阵的回代算法,如下,主要思想就是通过前面求得的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]
主过程:
For i<-1 to n-1
For j<-i+1 to n
m=a[j][i]/a[i][i]
a第j行中的每个元素-m*a第i行中的每个对应元素
(相当于上面的函数valueByRowAPlusToRowB)
Next j
Next i
至些,a[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 guassMethod(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 guassMethod,matrixArr is NULL\n");
assertF(bList!=NULL,"in guassMethod,bList is NULL\n");
assertF(xAnsList!=NULL,"in guassMethod,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 guassMethod,matrixA_b is NULL\n");
assertF(tmpMatrixArr!=NULL,"in guassMethod,tmpMatrixArr is NULL\n");
assertF(tmpBList!=NULL,"in guassMethod,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 guass 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 guass 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*************/
guassMethod(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 upTriangleMethod: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 upTriangleMethod:the ans x rows is:(from x0 to xn-1)
15.500000,37.500000,-23.000000;