由于求解三解方程较易，所以，考虑将系数矩阵A分解成两个三角矩阵的乘积，
即: A=LU的形式
其中，L为下三解矩阵，U为上三解矩阵，则线性方程组:
Ax=b
可改写为 LUx=b
令 Ux=y
得 Ly=b
然后，用前代方法求Ly=b得列矩阵y,再用回代方法求Ux=b，得到的列矩阵x即为所求的线性方程组的解.
分解的细节请参考相应的数值算法书籍，这里就不做过细的解释了.
首先是一个对下三解矩阵的前代方法，很简单的:
void downTriangleMethod(Type** matrixArr,Type* bList,Type* xAnsList,int len)
{
int i=0,j=0;
double e=0.0001;/*precise controller*/
Type tmpSum;
Type* tmpBList;
/*---assertion---*/
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");
/*mem apply*/
tmpBList=(Type*)malloc(sizeof(Type)*len);
for(i=0;i<len;i++)
tmpBList[i]=bList[i];
/*copy the data*/
for(i=0;i<len;i++)
assertF(fabs(matrixArr[i][i])>=e,"in upAngleMethod ,matrixArr[i][i] has value closed to 0/n");
/*data initialiate*/
tmpBList[0]=tmpBList[0]/matrixArr[0][0];/*Get the first answer*/
i=1;
/*Core part of the up triangle method*/
while(i<len)
{
tmpSum=0;
for(j=0;j<i;j++)
tmpSum+=matrixArr[i][j]*tmpBList[j];
tmpBList[i]=(tmpBList[i]-tmpSum)/matrixArr[i][i];
i++;
}
/*copy to answer*/
for(i=0;i<len;i++)
xAnsList[i]=tmpBList[i];
/*End of down triangle method*/
free(tmpBList);
}
/*The core part of the LUPation Method.*/
void LUPationMethod(Type** matrixArr,Type* bList,Type* xAnsList,int len)
{
/*Core think of this algorithm:
matrix_L*yAnsList=bList;
matrix_U*xAnsList=yAnsList;
*/
Type** matrix_L,** matrix_U;/*The core two matrix object of the LU method.*/
Type* yAnsList;
int i,j,k;/*iterator num*/
int i2,j2;
Type sum;/*the temp tween data*/
/*input assertion*/
assertF(matrixArr!=NULL,"in LUPationMethod,matrixArr is NULL/n");
assertF(bList!=NULL,"in LUPationMethod,bList is NULL/n");
assertF(xAnsList!=NULL,"in LUPationMethod,xAnsList is NULL/n");
/*memory apply*/
matrix_L=(Type**)malloc(sizeof(Type*)*len);
matrix_U=(Type**)malloc(sizeof(Type*)*len);
for(i=0;i<len;i++)
{
matrix_L[i]=(Type*)malloc(sizeof(Type)*len);
matrix_U[i]=(Type*)malloc(sizeof(Type)*len);
}
yAnsList=(Type*)malloc(sizeof(Type)*len);
/*==end of memory apply==*/
/*---assertion after apply the mem---*/
assertF(matrix_L!=NULL,"in LUPationMethod,matrix_L is NULL/n");
assertF(matrix_U!=NULL,"in LUPationMethod,matrix_U is NULL/n");
assertF(yAnsList!=NULL,"in LUPationMethod,yAnsList is NULL/n");
/*---end of assertion---*/
// printf("arr mem applyed/n");
/*----data initialization----*/
for(i=0;i<len;i++)
for(j=0;j<len;j++)
{
matrix_L[i][j]=0;
matrix_U[i][j]=0;
}
for(i=0;i<len;i++)
matrix_L[i][i]=1;
for(i=0;i<len;i++)
{
/*matrixU make*/
for(j=i;j<len;j++)
{
sum=0;
for(k=0;k<i;k++)
sum+=matrix_L[i][k]*matrix_U[k][j];
matrix_U[i][j]=matrixArr[i][j]-sum;
}
/*matrixL make*/
if(i<len-1)
{
j2=i;
for(i2=j2+1;i2<len;i2++)
{
sum=0;
for(k=0;k<j2;k++)
sum+=matrix_L[i2][k]*matrix_U[k][j2];
matrix_L[i2][j2]=(matrixArr[i2][j2]-sum)/matrix_U[j2][j2];
}
}
}
printf("matrix L:/n");
show2DArrFloat(matrix_L,len,len);
printf("matrix U:/n");
show2DArrFloat(matrix_U,len,len);
downTriangleMethod(matrix_L,bList,yAnsList,len);
upTriangleMethod(matrix_U,yAnsList,xAnsList,len);
/*memory free*/
for(i=0;i<len;i++)
{
free(matrix_L[i]);
free(matrix_U[i]);
}
free(matrix_L);
free(matrix_U);
free(yAnsList);
}
/*The matrix’s LUPation Method
This version is implemnted by EmilMatthew
You can use these code freely , but I'm not assure them could totally fit your needs.
*/
/*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*************/
LUPationMethod(matrixArr,bList,xAnsList,len);
printf("after the LUPationMethod:the ans x rows is:/n");
fprintf(outputFile,"after the LUPationMethod: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 upTriangleMethod: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;