[数值算法]列主元三角分解法

By EmilMatthew 05/9/13

列主元三角分角法是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,就是避免小数作为分母项.

有一些麻烦是的是要对原矩阵的增广矩阵部分进行LU分解,而且每次要对当前未处不景行作求列主行的下标,然后做变换.

其它的注意点和LU分解法其本一致,而最后求出U矩阵后,可以直接用回代法求出解x,而不必先求y再用Ux=y来求解向量x了.

其实,这里根本没有必要用到数学上推导时用的L和U两个三角矩阵,只要用一个增广矩阵就可以解决战斗了,不过我实现时还是用了老方法,以后有待改进.

/*

LUColMainMethod, coded by EmilMathew 05/9/13, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.

*/

void LUColMainMethod(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*/

/*LU

Col Main Spec Data*/

int maxRowIndex;

Type** matrixA_B;

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

twoDArrMemApply(&matrix_U,len,len+1);

for(i=0;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

{

for(j=0;j

{

matrix_L[i][j]=0;

matrix_U[i][j]=0;

}

matrix_U[i][j]=0;

}

for(i=0;i

matrix_L[i][i]=1;

/*matrix A_B prepare*/

twoDArrMemApply(&matrixA_B,len,len+1);

matrixCopy(matrixArr,matrixA_B,len,len);

for(i=0;i

matrixA_B[i][len]=bList[i];

printf("show copied matrix:/n");

show2DArrFloat(matrixA_B,len,len+1);

for(i=0;i

{

maxRowIndex=getMaxRowIndexLU(matrixA_B,i,matrix_L,matrix_U,len);

if(i!=maxRowIndex)

{

swapTwoRow(matrixA_B,i,maxRowIndex,len+1);

swapTwoRow(matrix_L,i,maxRowIndex,len);

swapTwoRow(matrix_U,i,maxRowIndex,len);

}

/*matrixU make*/

for(j=i;j

{

sum=0;

for(k=0;k

sum+=matrix_L[i][k]*matrix_U[k][j];

matrix_U[i][j]=matrixA_B[i][j]-sum;

}

matrix_U[i][j]=matrixA_B[i][j]-sumArr1_IKByArr2_KJ(matrix_L,matrix_U,i,j,0,i-1);

matrixA_B[i][j]=matrix_U[i][j];

/*matrixL make*/

if(i

{

j2=i;

for(i2=j2+1;i2

{

sum=0;

for(k=0;k

sum+=matrix_L[i2][k]*matrix_U[k][j2];

matrix_L[i2][j2]=(matrixA_B[i2][j2]-sum)/matrix_U[j2][j2];

}

}

}

//copyBack

for(i=0;i

bList[i]=matrixA_B[i][len];

printf("show bList/n");

showArrListFloat(bList,0,len);

/*Adjusting of matrix_L*/

for(i=0;i

{

for(j=i;j

{

if(i==j)

matrix_L[i][i]=1;

else

matrix_L[i][j]=0;

}

}

printf("matrix L/n");

show2DArrFloat(matrix_L,len,len);

printf("matrix U/n");

show2DArrFloat(matrix_U,len,len);

for(i=len-1;i>=0;i--)

{

xAnsList[i]=(bList[i]-sumArr_JKByList_K(matrix_U,xAnsList,i,i+1,len-1))/matrix_U[i][i];

}

/*

downTriangleMethod(matrix_L,bList,yAnsList,len);

upTriangleMethod(matrix_U,yAnsList,xAnsList,len);

*/

/*memory free*/

for(i=0;i

{

free(matrix_L[i]);

free(matrix_U[i]);

}

free(matrix_L);

free(matrix_U);

free(yAnsList);

free(matrixA_B);

}

//辅助函数,求列主元的行位置.

int getMaxRowIndexLU(Type** inMatrixArr,int currentI,Type** matrix_L,Type** matrix_U,int size)

{

int i,maxRowIndex;

Type tmpSMax,tmpCompare;

assertF(inMatrixArr!=NULL,"in getMaxRowIndexLU,inMatrixArr is NULL/n");

assertF(matrix_L!=NULL,"in getMaxRowIndexLU,matrix_L is NULL/n");

assertF(matrix_U!=NULL,"in getMaxRowIndexLU,matrix_U is NULL/n");

tmpSMax=0;

//maxRowIndex=currentI;

for(i=currentI;i

{

tmpCompare=(float)fabs(inMatrixArr[i][currentI]-sumArr1_IKByArr2_KJ(matrix_L,matrix_U,i,currentI,0,currentI-1));

if(tmpCompare>tmpSMax)

{

tmpSMax=tmpCompare;

maxRowIndex=i;

}

}

return maxRowIndex;

}

我所写的程序的源码下载:

//如果上面这个链接无法响应下载(有可能是被网站给屏蔽掉了)，则可使用下载工具(如迅雷等)下载

测试结果:

test1:

input:

3;

4,-2,-4,10;

-2,17,10,3;

-4,10,9,-7;

output:

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

2.000001.00000-1.00000

test2:

input:

3;

12,-3,3,15;

18,-3,1,15;

-1,2,1,6;

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

1.000002.000003.00000