[数值算法]高斯消元法改进版---列主消元法
高斯消元在求解一些系数矩阵中含有极小数的情况下，会产生巨大的舍入误差，导致算法失效。一个简单而有效的改进方法是每次在进行将当前列中元素的消成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;
